Revisiting Horizontal Diffusion of Perturbations over Terrain for NCEP RSM

1 Department of Computer Science and Information Engineering, Nan Jeon Institute of Technology, Tainan, Taiwan, ROC 2 Chung Cheng Institute of Technology, National Defense University, Taoyen, Taiwan, ROC 3 Environmental Modeling Center, National Centers for Environmental Prediction, Washington, DC, USA * Corresponding author address: Dr. Hann-Ming Henry Juang, Environmental Modeling Center, National Centers for Environmental Prediction, Washington, DC, USA; E-mail: Henry.Juang@noaa.gov ̇The perturbation in the National Centers for Environmental Prediction (NCEP) Regional Spectral Model (RSM) is defined by the difference between the RSM and its outer model, such as the NCEP global model, on the same given coordinate surfaces. This perturbation is dominated by the difference between terrains of the RSM and the global model if model terrain exists over the domain. After being transformed into spectral space, this difference in terrains is characterized by two distinct groups of coefficients i.e., long-wave and short-wave groups. The long-wave group is considerably larger in amplitude than the short-wave group, and considered to be the cause of error in horizontal diffusion of perturbations on sigma surfaces. A higher order of diffusion and diffusion without long waves on sigma surfaces with perturbations are formulated to investigate the possibility of using diffusion on sigma surfaces. In a typical case, erroneous circulation occurred when fourth-order diffusion on a sigma surface was used; however, reasonable circulation was simulated with second-order diffusion on the pressure surface. Nevertheless, diffusion on a pressure surface is a non-economical method in spectral modeling. The results of sixthand eighth-order horizontal diffusions of perturbations and no-diffusion on the group of long waves for all orders (second, fourth, sixth, and eighth) on sigma surfaces did not show erroneTerr. Atmos. Ocean. Sci., Vol. 18, No. 1, March 2007 68 ous circulations. Thus, diffusions greater than and equal to sixth-order diffusions and no-diffusion on a long-wave group can be implemented into NCEP RSM as an economical computation as it leads to an absence of erroneous circulation comparable to diffusion on pressure surfaces. (

ous circulations.Thus, diffusions greater than and equal to sixth-order diffusions and no-diffusion on a long-wave group can be implemented into NCEP RSM as an economical computation as it leads to an absence of erroneous circulation comparable to diffusion on pressure surfaces.

INTRODUCTION
Some numerical schemes can result in smoothed prognostic fields, such as an interpolation in a semi-Lagrangian scheme (Robert 1993) or enstrophy conserving finite difference (Arakawa 1966); however, most discretized numerical models produce computational noise, which is generated through short wave growth due to energy cascade and/or aliasing of nonlinear integration.This kind of noise could grow into computational instability (Phillips 1959), so that integration cannot be continued.There are two common methods to control this kind of noise, one is subscale parameterization (Tag et al. 1979;Leith 1968), and the other is horizontal diffusion (Shapiro 1970;Cullen 1976;von Storch 1978).
Horizontal diffusion on terrain-following coordinates produces erroneous vertical mixing (Pielke 1984) when such diffusion is on coordinate surfaces over mountains.Solutions to this include applying horizontal diffusion on constant-height or quasi-horizontal surfaces (e.g., Juang et al. 2005;Zängl 2002) and applying perturbation diffusion to terrain-following coordinates.As mentioned in Juang et al. (2005, hereafter J05), having horizontal diffusion on pressure surfaces in spectral models is complicated.Horizontal diffusion on pressure surfaces from variables on sigma surfaces requires coordinate transformation or vertical interpolation.Furthermore, in spectral computation there is more than one term in coordinate transformation, and most of the corrected terms are nonlinear, requiring arduous computation in grid-point space.Thus it is not pragmatic to use higher-order diffusion because of a geometric increase in nonlinear terms.
The problem of horizontal diffusion shows clearly in temperature and humidity fields because they are usually vertically stratified, either normal (decreasing with height) or inverted (increasing with height).As shown in Fig. 1 of J05, under normal stratification, diffusion on a sigma surface will warm and moisten the atmosphere over mountaintops yet cool and dry it over slopes and/or foothills.Most horizontal diffusion schemes use perturbations of temperature and/or humidity with respect to certain reference fields such as NCEP GFS (Global Forecast System, Kanamitsu et al. 1991) and CWB (Central Weather Bureau) GFS (Liou et al. 1997).The global horizontal mean field of temperature is used for NCEP GFS to do perturbation diffusion on temperature, and time-mean 3-dimensional temperature is used for horizontal diffusion in CWB GFS.Thus, it is a common problem in the spectral model to use perturbation diffusion.Nevertheless, J05 demonstrated erroneous circulations with fourth-order horizontal diffusion on sigma surfaces with perturbations in temperature and humidity.As a consequence of this study it became apparent that further investigation of the characteristics of perturbation is required to determine why horizontal diffusion of perturbations on sigma surfaces failed.This research continues that of J05 and revisits horizontal diffusion of perturbations on sigma surfaces over mountains.Section 2 illustrates perturbation diffusion on sigma surfaces, perturbation characteristics, and the high-order form of diffusion on sigma surfaces.Section 3 describes the typical case selected.Section 4 describes the results from horizontal diffusion on perturbations using different orders and from horizontal diffusion without a long-wave group.Section 5 is the summary and conclusion.

HORIZONTAL DIFFUSION OF PERTURBATION ON SIGMA SURFACES
In this section, we briefly describe the model, investigate the characteristics of perturbations, paraphrase the design of horizontal diffusion of perturbations on a sigma surface, and illustrate different orders of diffusion with and without long waves.

Model Description and the Definition of Perturbation
The National Centers for Environmental Prediction (NCEP) Regional Spectral Model (RSM) was used in this study.A description of the model can be found in Juang and Kanamitsu (1994) for the hydrostatic version and in Juang (2000) for the non-hydrostatic version.The perturbation method in spectral computation, described in detail in Juang and Kanamitsu (1994) and Juang et al. (1997), provides a spectral truncation as well as filtering to preserve largescale features; this approach is beneficial in long-term integrations (Juang and Hong 2001).
A definition of a perturbation is a deviation between the inner regional model and the outer global model on the same sigma surface.The sigma surface in the RSM and in the outer model may not be at the same physical height because terrain heights in the two models may differ.In this case, the perturbation is a mathematical, not a physical, value.Figure 1 of J05 shows that a different model terrain will provide perturbations over all sigma surfaces.In other words, model perturbation is related to the difference in terrain between the NCEP RSM and the outer model.If horizontal diffusion is applied to this perturbation on the sigma surface in spectral space, perturbations preexisting due to the difference in model terrains may introduce a systematic error that will grow with time integration.Because the main goal of this study is to revisit diffusion on sigma surfaces in spectral space with perturbations, we had to identify the characteristics of this kind of perturbation in spectral space.

Characteristics of Perturbation
In cases of model-terrain proximity between the NCEP RSM and its outer, coarse-resolution model -for example, the NCEP Global Spectral Model (GSM), the perturbation will initially be negligible and the noise generated after nonlinear integration may be smoothed out by horizontal diffusion on sigma coordinate surfaces with few diffusion errors.Although using a perturbation field for horizontal diffusion on sigma surfaces might result in fewer errors than using the full variable, the error might be non-negligible near steep mountain areas when the difference in terrain between the regional and outer models is sufficiently large due to increasing regional resolution.
Since, RSM uses rectangular truncation, we can plot all wave coefficients with respect to x-direction wave numbers and y-direction wave numbers, or we can plot coefficients with respect to total wave numbers, by combining x-direction and y-direction wave numbers.This is similar to the more familiar triangular truncation for global modeling.In triangular truncation, we have zonal wave numbers, meridian wave numbers, and total wave numbers.As terrain perturbation exists (J05) and it is temporally independent, if the problem is present, it should become clear during entire integration.Thus, let's examine terrain perturbation with plots as mentioned.Figure 1 shows the logarithm of the perturbation in x and y waves with respect to two-dimensional wave spaces in relation to the difference in terrain between inner regional and outer global models for four different domains.Domain (a) covers South China from 112 to 118.5°E and from 25.5 to 31°N with 35 × 33 wave numbers; domain (b) covers Taiwan from 118 to 124°E and from 21 to 26.5°N with 35 × 44 wave numbers; domain (c) covers the southwestern US and Mexico from 99 to 117.5°W and from 20 to 37.5°N with 57 × 57 wave numbers; and domain (d) covers Hawaii from 158.8 to 157.4°W and from 21 to 22.1°N with 79 × 77 wave numbers.There is a high gradient of contours in each panel of Fig. 1.These contours along with frame (as in x and y axies) form a rectangular box.Within the box, the wave numbers are smaller (long waves), and the magnitude of the coefficients is larger than those out of the box, which has all short waves.The boundary that separates the long and short waves shows a wave number roughly at 0.45 times the total wave number in x direction, and at 0.55 times total wave number in y direction for (a), (b), and (c), but not (d).The boundary in Fig. 1d may be roughly 0.5 for both directions.Note that, the feature shown here is under the existence of terrain perturbation, and is regardless of whether terrain perturbation has or has not been subject to spectral smoothing by Lanczos' method.The terrain data is NCAR NAVY at 10 min by 10 min resolution.Furthermore, we checked spectral coefficient plots for surface pressure perturbation and all temperature and humidity.Scale separation exists on surface pressure perturbations, but is not so clear in temperature and humidity.However, coefficients of temperature and humidity have larger magnitude for longer waves and smaller magnitude for short waves, indicating scale dependentce in amplitude though there is no clear line of separation between the two groups.
The waves in Fig. 1 separate into two groups in x and y space; these can be plotted to combine the two groups of waves in x and y space into one, as mentioned and shown in Fig. 2 (domains as in Fig. 1).No matter which domain is used, the perturbation shows two groups of magnitudes in Fig. 2. The gap in Fig. 2 resembles the high-gradient boundary feature in Fig. 1.The group of long waves dominates the perturbation and the group of short waves contributes little or negligible amplitudes of perturbation.This finding led us to investigate horizontal diffusion in terms of no or negligible diffusion in the group of long waves, because the longwave group may contribute systematic errors due to a larger amplitude in spectral space.However, the same definition of the separation of long and short waves should be used as that described for Fig. 1.In Fig. 2, for example, wave number 20 in panel (a) has three coefficients in the large-magnitude group with its remaining coefficients being in the small-magnitude group.It is not clear whether total wave number 20 belongs in the long-wave group or shortwave group.For rectangular waves, we are using wave numbers in x and y space to group long waves and short waves with respect to their magnitudes.

The General Form of Horizontal Diffusion on Sigma Surfaces and Diffusion Without Long Waves
The NCEP RSM/Mesoscale Spectral Model (MSM) applies fourth-order horizontal diffusion on sigma surfaces to the perturbation in spectral space (Juang and Kanamitsu 1994).The perturbation equation with horizontal diffusion can be paraphrased here as: where A can be any prognostic variable and κ is a constant diffusion coefficient; A R is a regional variable, A G is a global variable, and A' is a perturbation which should be A R -A G as defined.The first term on the right-hand side (RHS), ∂ ∂ A t R * / , is the total forcing computed on the regional grid without horizontal diffusion.The second term on the RHS, ∂ ∂ A t G / , is the total tendency (including outer-grid diffusion by the base field) from the outer model interpolated on the regional grid.The third term on the RHS is the horizontal diffusion of perturbations on sigma surfaces, assuming that the horizontal diffusion of the base field from the outer model was done within the second term on the RHS (more detail in Juang and Kanamitsu 1994).
The perturbation horizontal diffusion on coordinate surfaces in Eq. ( 1) is a linear term and can be computed in spectral space.We followed the procedure discussed in Sardeshmukh and Hoskins (1984) to rewrite the last term in Eq. ( 1) in spectral form for each wave as: where m and n are the wave numbers in x and y directions, respectively, M and N are the maximal wave numbers in x and y directions, respectively, and h is the order of diffusion.

The [ ']
A m n represents the variable amplitude in spectral space at wave (m, n) and τ σ is e-folding time as a weighting coefficient to control the effect of the diffusion so that κ π τ σ = 1 / [( x / ) ] h/2 ∆ .For numerical stability, an implicit time scheme is used as: and it can be solved as: This form is unconditionally stable for any value of τ σ and can have a strong diffusion, while τ σ is as small as two time step.For diffusion with the group of short waves only, a different order of diffusion on sigma surfaces can be used with a modification of µ = 0 when the wave is less than or equal to 0.45 times the maximal wave number in the x direction and less than or equal to 0.55 times the maximal wave number in the y direction.Figure 3 shows the contributions of the (a) second-, (b) fourth-, (c) sixth-, and (d) eighth-order diffusions for all waves; Fig. 4 shows the (a) second, (b) fourth, (c) sixth, and (d) eighth orders without diffusion for long waves.Those marked by a gray X are coefficients of terrain perturbation, those marked by a gray + are diffusion coefficients, and those marked by a black dot are diffusion coefficient (gray +) multiplied by perturbation coefficient (gray X).The higher the order of diffusion for all waves, the less the diffusion influences the group of long waves.Thus, a higher order of diffusion may have less bias error over mountainous areas.In the case of no diffusion on long waves (from wave 0 to wave 15 in Fig. 4), there may be less or no bias error as well.Nevertheless, the results from these different orders of diffusion with and without diffusion on long waves will be examined in Section 4.

CASES DESCRIPTION AND EXPERIMENTAL DESIGNS
Because J05 showed that stable conditions over mountainous areas have the best evidence of erroneous land-sea breezes from horizontal diffusion, we selected a stable-atmosphere case to investigate the bias effect of horizontal diffusion on sigma surfaces.Though we ran several other cases, this one revealed the clearest evidence of the impacts of perturbation diffusion.
The area selected was over the Oahu Island of Hawaii (see Fig. 5 for modeled horizontal domain); the initial conditions were from NCEP/NCAR reanalysis (Kalnay et al. 1996) on 0000 UTC 23 May 2002(1400HST 22 May 2002).The boundary conditions were from the reanalysis as well (at six-hour intervals).Reanalysis has horizontal resolution of T62 and 28 vertical layers.Three RSM domains were designed at resolutions of 60 km, 10 km, and 1 km of horizontal grid-spacing with time steps of 240, 60, and 5 seconds, respectively.The domain of 60-km resolution was over the area from 174 to 144°W and from 6 to 35°N; the 10-km resolution domain was over the area from 161.5 to 153.5°W and from 18 to 25°N; and the 1-km resolution domain was over the area from 158.8 to 157.4°W and from 21 to 22.1°N.Synoptic conditions for this case were shown in J05, which indicated stable conditions with a synoptic high.Note that, though elevation of the mountains at Oahu is not very high, the effect of improper horizontal diffusion producing erroneous circulation is apparent in this selected case.
The results from p-DIFF, a second-order horizontal diffusion of full fields on pressure surfaces (see J05), were treated as a standard solution.Results shown in Section 4 are from only the innermost domain with nesting into outer domains, which uses p-DIFF as the horizontal diffusion.The innermost domain uses either p-DIFF or different diffusions of perturbations on sigma surfaces.Figure 5 shows (a) the lifted index from model integrations, and (b) a 10-m wind (10 m above the model terrain resembles the height of a surface observation station) with p-DIFF.This figure indicates a stable atmosphere as shown in observations over the mountainous areas, and the p-DIFF shows the corrected circulation with a sea breeze during the daytime.

RESULTS
Results from only the innermost domain with different diffusions are shown here.Landsea breeze is the major phenomenon used to verify the results.A 10-m wind indicates model land-sea breezes.

Fig. 4.
As in Fig. 3, but without diffusion on the group of long waves.

Comparison Among Different Order of Perturbation Diffusions
Figure 6 shows the 10-m wind after six-hour integrations from (a) second-order, (b) fourthorder, (c) sixth-order, and (d) eighth-order diffusions on sigma surfaces.Again, the erroneous land breezes over the east coast and ocean areas are shown in the second-and fourth-order diffusions as expected.During the daytime with a stable atmosphere, the sea breeze should be over the shore and along the coastal areas.Results from the sixth order and eighth orders show the same sea-breeze circulation as from the p-DIFF.
Figure 7 shows the vertical cross section of virtual temperature from experiments of (a) p-DIFF, (b) fourth-order diffusion, and (c) sixth-order diffusion.Fourth-order diffusion produced a wave-like field over the mountains and showed warming on top of the mountains and cooling over mountain slopes as indicated in Fig. 1 of J05.Results from p-DIFF and sixthorder diffusion showed no wave-like distribution over mountain areas, which indicates that a higher order of horizontal diffusion (sixth-order) with perturbations produced no erroneous warming or cooling fields.

Sensitivity Test with No Diffusion on a Group of Long Waves
Based on the assumption described in Section 2, horizontal diffusion on sigma surfaces can be used if the diffusion effect over the group of long waves is negligible and will not accumulate and produce an erroneous circulation over mountainous areas.The question now is whether there is no erroneous circulation if only the group of short waves is diffused.Figure 8 shows a 10-m wind after a six-hour integration from the (a) second-, (b) fourth-, (c) sixth-, and (d) eighth-order diffusions of perturbations on sigma surfaces with no diffusion on a group of long waves in the Hawaii case.These results may be compared with p-DIFF in Fig. 5 and with the diffusion of all waves in Fig. 6.There is no erroneous circulation of land breeze in any results in Fig. 8, even for the second-and fourth-order diffusions.Even though there may be some noise in the second order, the results support our concerns of grouping perturbations in spectral space, and it proves that the erroneous circulation on horizontal perturbation diffusion on sigma surfaces is due to the influence of diffusion on the group of long waves.Fig. 8.As in Fig. 6, except for no diffusion on the group of long waves.

CONCLUSIONS
Even with perturbations, horizontal diffusion on sigma surfaces produce erroneous circulations, thus, horizontal diffusion on quasi-horizontal pressure surfaces was proposed and implemented in the NCEP RSM (J05).Nevertheless, implementation of horizontal diffusion on pressure surfaces to spectral models contains time-consuming spectral transformations and complicated terms of coordinate transformations, which limit an extension to higher orders of diffusion.Thus, horizontal diffusion with perturbations on sigma surfaces was revisited in this paper and it was revealed that erroneous circulation of perturbation diffusion on sigma surface is due to horizontal diffusion on the long-waves of terrain perturbation, which we prove by diffusion without a long-wave group for all orders of horizontal diffusion.
Characteristics of terrain perturbations were investigated to determine if it is possible to use horizontal diffusion with perturbations on sigma surfaces.We found that, for any given domain containing terrain, terrain-perturbation amplitudes are distributed into two groups: larger perturbations into a group of long waves and small perturbations into a group of short waves.This two-group distribution may result mathematically from spectral discretization of original terrain data.Thus, no matter what domain with mountains is given, the perturbation distribution in wave spaces shows two groups of spectral coefficients.Furthermore, the higher the order of horizontal diffusion, the less diffusion appears in the group of long waves.Thus, we examined high-order horizontal diffusion with perturbations on sigma surfaces, and zero diffusion coefficients within the group of long waves.
Exploration of different orders of horizontal diffusion of perturbations on sigma surfaces indicated that sixth-order diffusion is good enough to smooth noise without causing the group of long waves to produce erroneous circulation; and second-and fourth-orders produce erroneous circulations.Thus, sixth-order diffusion can be used as optimal horizontal diffusion with perturbations and the NCEP RSM can take advantage of spectral computation of horizontal diffusion.
The concept of no diffusion for the group of long waves produced positive results for all orders of diffusion, including the second and fourth orders.This indicates that the concept of revisiting horizontal diffusion on sigma surfaces through the characteristics of terrain-perturbation distributions in spectral space is acceptable.This method of horizontal diffusion can be used in any order, especially second-order diffusion to allow strong diffusion on short waves.This wave-selected diffusion (diffusion only on short waves) has been used in NCEP GFS as well.Based on theoretical and numerical findings of Leith (1968 and1971), NCEP GFS has diffusion on short waves starting at a wave number equal to 0.55 times total wave number.The purpose of this wave-selected diffusion is to avoid bias in diffusion on long waves, but to allow strong diffusion on short waves.However, such wave-selected diffusion can only be used on global and regional spectral models.Furthermore, we note that our finding is a numerical method based on real data, not based on any physical meaning or mathematical proof, thus, the factors (for example, 0.45 in x-direction and 0.55 in y-direction) may change in different domains as shown in Fig. 1d, even though it is possible to have no gap if the entire domain has no terrain perturbation between the global and regional model, such as over oceans or flat plains in the entire regional domain.However, the existence of a gap in terrain perturba-tion provides us a starting point to separate long waves and short waves especially for horizontal diffusion, or possibly for other wave-related numerical techniques.Although, we have shown only one dry stable case, conclusions from other moist unstable cases should be the same, because this scale selected diffusion relates only to terrain perturbation in spectral space, not to the stability of the atmosphere.

Fig. 1 .
Fig. 1.Spectral coefficients of terrain perturbation with respect to two-dimensional wave numbers for the domains of: (a) South China, (b) Taiwan, (c) south-east US and Mexico, and (d) Hawaii.

Fig. 2 .
Fig. 2. Spectral coefficients of terrain perturbation with respect to combined wave numbers for the domains of: (a) South China, (b) Taiwan, (c) south-east US and Mexico, and (d) Hawaii.

Fig. 3 .
Fig. 3. Spectral coefficients of terrain perturbations (symbol x), diffusion coefficients (symbol +), and the multiplication of both (symbol period) with respect to wave numbers are plotted for (a) second-, (b) fourth-, (c) sixth-, and (d) eighth-order horizontal diffusion on terrain surfaces.

Fig. 5 .
Fig. 5. (a) Lifted-index of model integration and (b) 10-m wind after 6 hours of integration by p-DIFF.Shaded areas are model terrain at interval of 300 m.