## Abstract

Although a terrain-following vertical coordinate is well suited for the application of surface boundary conditions, it is well known that the influences of the terrain on the coordinate surfaces can contribute to increase numerical errors, particularly over steep topography. To reduce these errors, a hybrid sigma–pressure coordinate is formulated in the Weather Research and Forecasting (WRF) Model, and its effects are illustrated for both an idealized test case and a real-data forecast for upper-level turbulence. The idealized test case confirms that with the basic sigma coordinate, significant upper-level disturbances can be produced due to numerical errors that arise as the advection of strong horizontal flow is computed along coordinate surfaces that are perturbed by smaller-scale terrain influences. With the hybrid coordinate, this artificial noise is largely eliminated as the mid- and upper-level coordinate surfaces correspond much more closely to constant pressure surfaces. In real-data simulations for upper-level turbulence forecasting, the WRF Model using the basic sigma coordinate tends to overpredict the strength of upper-air turbulence over mountainous regions because of numerical errors arising as a strong upper-level jet is advected along irregular coordinate surfaces. With the hybrid coordinate, these errors are reduced, resulting in an improved forecast of upper-level turbulence. Analysis of kinetic energy spectra for these simulations confirms that artificial amplitudes in the smaller scales at upper levels that arise with the basic sigma coordinate are effectively removed when the hybrid coordinate is used.

## 1. Introduction

Topography is one of the main sources for atmospheric gravity waves, and simulating their behavior in atmospheric modeling systems is crucial for numerical weather prediction. With the availability of increasing computational resources, atmospheric numerical models are progressing to higher resolutions and are investigating the influences of small-scale terrain on gravity wave dynamics (e.g., Kirshbaum et al. 2007; Sheridan and Vosper 2012). Historically, terrain-following vertical coordinates have been widely used in atmospheric models to represent flow over topography, although it is well known that the terrain influences on the coordinate surfaces can contribute to increased errors in the model numerics, particularly in computing the horizontal pressure gradient terms (Mahrer 1984; Dempsey and Davis 1998). The numerical treatment of real topography is challenging because the terrain may contain significant structure at the smallest scales that are represented in the model, scales for which the model numerics are least accurate. Park et al. (2016, hereafter P16) assessed the influence of terrain smoothing on upper-level turbulence forecasting and documented that the overprediction of upper-air turbulence can be significantly reduced by filtering the high-frequency terrain structure.

To mitigate errors associated with a terrain-following coordinate, there have been a number of studies examining alternative numerical approaches, including alternate numerics to increase the accuracy in calculating the horizontal pressure gradient (cf. Mahrer 1984) and cut-cell techniques to remove terrain influences from the coordinate surfaces (cf. Adcroft et al. 1997; Steppeler et al. 2002). A hybrid terrain-following vertical coordinate provides another approach in reducing numerical errors by more rapidly removing the influences of the terrain on the coordinate surfaces with increasing height, such that they transition to constant height (or pressure) surfaces at mid- and upper levels in the model domain (Simmons and Burridge 1981; Bubnová et al. 1995; Schär et al. 2002; Klemp 2011).

The Advanced Research WRF (WRF-ARW) Model solves the nonhydrostatic equations employing a terrain-following hydrostatic pressure (sigma) vertical coordinate, often characterized as a mass coordinate (Laprise 1992; Klemp et al. 2007; Skamarock et al. 2008). In this study, we formulate a hybrid coordinate for the WRF-ARW following the design in Park et al. (2013) and document its impact in both an idealized and a real-data simulation. For the real-data case, we investigate the role of the hybrid coordinate in the upper-level turbulence forecast analyzed by P16 in assessing the effects of terrain smoothing. Since upper-level turbulence forecasting is quite sensitive to topography and high wind shear at upper levels, this case provides a good test of the dependence of the coordinate formulation on the accuracy of the model numerics over terrain.

Our paper is organized as follows. In section 2, formulation of the hybrid coordinate in the WRF-ARW is introduced, and an idealized test case is presented in section 3 to illustrate the performance of the hybrid coordinate. Section 4 provides an analysis of a real-data forecast for air turbulence and demonstrates that the hybrid vertical coordinate reduces artificial noise (and associated overprediction of turbulence) at mid- and upper levels. The kinetic energy spectra for these simulations provide further confirmation of the beneficial influence of the hybrid coordinate. Section 5 summarizes the results of this study.

## 2. Formulation of the hybrid coordinate in the WRF mass coordinate

We employ a terrain-following hybrid sigma–pressure coordinate in the nonhydrostatic WRF-ARW Model using a formulation similar to Park et al. (2013), who introduced the hybrid, mass vertical coordinate in the hydrostatic version of the Model for the Prediction across Scales (MPAS). In this hybrid coordinate system, the vertical coordinate *η* is defined in terms of the hydrostatic pressure of dry air :

where is a reference sea level pressure, and are the hydrostatic surface pressure and the top-level pressure for dry air, respectively. Here, defines the relative weighting of the terrain-following and pure dry hydrostatic pressure coordinate such that

Note that for in (2), the coordinate formulation is the same as original mass coordinate in WRF, which was described by Klemp et al. (2007). To smoothly transition from a sigma coordinate near the surface to a pressure coordinate at upper levels, we define using a third-order polynomial

subject to the boundary conditions

such that

Here, is a specified value of *η* at which it becomes a pure pressure coordinate, and thus for all . Figure 1 displays the profiles for the basic terrain-following (BTF) sigma coordinate and for the hybrid terrain-following (HTF) coordinate for several values of the parameter . Plotted as a function of *η* (Fig. 1a), these profiles depict the form of the polynomial defined in (3)–(5). However, plotting as a function of height (Fig. 1b) provides a better physical sense of the transition toward a pressure coordinate with increasing height. For example, with a model domain having a depth of 30 km, for , the vertical coordinate becomes a pure pressure coordinate at an altitude of about 12 km.

With the vertical coordinate definition (1), the coordinate metric term has the formulation

where is proportional to the mass in each vertical column. Defining contravariant flux-form variables as , the continuity equation is then expressed as

where . Vertically integrating from the surface to the top of the domain yields

Stepping (8) forward in time for , the coordinate metric is then recovered from (6). The vertical wind component is also derived from vertical integration of (6) level by level:

In implementing the hybrid mass coordinate in the WRF-ARW, there is little change in the model equations and numerics. Essentially, the only change is that vertically integrating the horizontal divergence in (8) now yields at the new time level (instead of ), and is then recovered from (6). Thus, the only change is that becomes a 3D function instead of a 2D function as in the BTF formulation. The 3D array need not be stored, as it can be readily reconstructed wherever it is used from (6), which simply scales the 2D array by 1D functions of *η*. A detailed description of the WRF-ARW Model equations and numerics is contained in Skamarock et al. (2008), and an updated version of their tech note (version 4) will soon be released that includes description of the hybrid-coordinate implementation.

With the hybrid coordinate, as the terrain influence is removed more rapidly from the coordinate surfaces, the coordinate surfaces become more compressed above the higher terrain. In the extreme, this could cause adjacent coordinate surfaces to intersect, which violates the monotonicity requirement of the coordinate. This situation will be avoided provided at all levels above the location of the minimum surface pressure (i.e., the highest terrain). Using (3)–(6), this monotonicity requirement is equivalent to

where represents the minimum surface pressure. The maximum value of is readily obtained by setting the *η* derivative of (3) equal to zero and yields the requirement

where the coefficients defined in (5) are functions of only. The smallest allowable surface pressure for which the hybrid coordinate does not violate monotonicity [applying the = sign in (11)], designated as is plotted in Fig. 2 as a function of and . However, in order to maintain solution accuracy and stability, it is important to configure the hybrid coordinate such that the anticipated minimum surface pressure is significantly greater than In WRF-ARW version 4, the hybrid coordinate is employed with a default value of

## 3. An idealized test case

We explore an idealized 2D test case for upper-level flow over terrain to evaluate performance for the hybrid coordinates as described in the previous section. Following Schär et al. (2002), the terrain shape is specified as

with , , and m. The wind profile is

and a two-layer stability profile is defined with a Brunt–Väisälä frequency *N* as

The top of the model domain is set at approximately 20 km. Two different vertical coordinates, BTF and HTF, are tested, and their level distributions are shown in Fig. 3. For the HTF coordinate as represented in (1), the profile is defined using (3)–(5) with . The influence of the topography on the coordinate surfaces is still significant at upper levels in the domain for the BTF coordinate (Fig. 3a), whereas with the HTF coordinate with (Fig. 3b), vertical levels become flat at about 11 km. The simulations are conducted on a grid with m and with vertical increments in *η* that correspond to a nearly constant vertical mesh spacing of about 500 m. The model equations are integrated forward in time using split-explicit numerics as described by Klemp et al. (2007) and Skamarock et al. (2008) using a 12-s time step and a 2-s acoustic step.

Because there is no wind at or below mountain-top level in this test case, the upper-level flow should remain undisturbed as it passes over the mountain. Perturbations of vertical velocity at 5 h are depicted in Fig. 4 and exhibit significant differences. For the BTF simulation in Fig. 4a, stationary waves are excited at upper levels in the domain with maximum vertical velocities reaching about 1 m s^{−1}. In contrast, with the HTF coordinate, the influence of topography is largely removed from the coordinate surfaces in the upper half of the domain, and the perturbations are much weaker, with the maximum vertical velocity less than 0.07 m s^{−1}. This significant reduction of the spurious perturbations is to be expected because for the HTF case, the coordinate surfaces are much more closely aligned with the direction of the flow (see Schär et al. 2002; Shaw and Weller 2016).

## 4. Real-data test case

### a. A case on 2 November 2015

To explore the influences of the vertical hybrid coordinate in NWP applications, we have conducted a retrospective simulation for a case over the western United States on 2 November 2015, previously studied by P16. Here, the WRF-ARW is configured as close as possible to the operational NOAA Rapid Refresh (RAP) model, which predicted a huge false alarm of the light intensity of upper-level turbulence that was not supported by numerous null (smooth) observations.

As shown in the 500-hPa analysis in Fig. 5, a large-scale synoptic trough extended from the eastern Pacific through British Columbia at 0000 UTC 2 November 2015. Large gradients of geopotential heights persisted along the southern portion of planetary waves, and strong westerly winds were present in the western part of the Rockies. During the next 24 h, the wind direction changed to southwesterly as the synoptic trough approached the western part of Washington (not shown). For this case, we have run the WRF-ARW forecast over the United States with the same configuration as used for the P16 simulation initialized at 0000 UTC 2 November 2015, but including the hybrid vertical coordinate as specified in (1) with given by (3)–(5) for . A single domain with a 13-km grid covering the contiguous United States (CONUS) using 758 × 567 mesh points was used. The initial fields were interpolated from GFS data, and the same physics parameterizations were used as in the operational RAP (see http://rapidrefresh.noaa.gov). In varying the smoothness of the terrain, P16 found that the filtering of topography had a significant influence on mitigating the false alarms of light intensity turbulence over the western mountainous United States. In this study, we extend the tests from P16 to determine the impact of including the hybrid vertical coordinate. The list of simulations is summarized in Table 1. RAP refers to simulations with the WRF-ARW similar to the operational NOAA RAP, and RAP-HYBRID is the same configuration except for inclusion of the hybrid vertical coordinate. As introduced in P16 for simulations with terrain smoothing, SMTH1 uses one pass of a fourth-order terrain filter (the default topography shape in the WRF-ARW), and SMTH2 and SMTH4 have increased smoothing of the topography with two and four passes of the terrain filter, respectively. Each of the different amounts of terrain smoothing is also tested with the hybrid coordinate as indicated in Table 1 to investigate the role of the vertical coordinate in upper-level turbulence and mountain-wave prediction.

### b. Turbulence forecast

To illustrate the influences of terrain smoothing and the hybrid coordinate on the shape of the coordinate surfaces, vertical cross sections of these surfaces are displayed in Fig. 6 along the line AB in Fig. 5. Clearly, the coordinate surfaces with the unfiltered terrain exhibit much more small structure (Fig. 6a) than with the heavily filtered terrain (Fig. 6c). However, regardless of the terrain smoothing, the terrain influence on upper-level coordinate surfaces remains significant. The upper level influences of terrain are essentially eliminated with the hybrid vertical coordinate, as shown in Figs. 6b and 6d for the RAP-HYBRID and SMTH4-HYBRID cases, respectively. Terrain features influencing the coordinate surfaces are rapidly reduced and become almost negligible at levels above 11 km, even with high-frequency topography in RAP-HYBRID.

The 18-h forecast for the calculated eddy dissipation rate (EDR; m^{2/3} s^{−1}) turbulence and wind speed at 35 000 ft (10 668 m) are displayed in Fig. 7. Here, EDR is derived from the absolute vertical velocity divided by the Richardson number (), and more detailed information can be found in Sharman and Pearson (2017) and Kim et al. (2018).

Pilot reports (PIREPs) are also shown in Fig. 7. Here, flight observations of null turbulence (EDR < 0.05) are displayed with black dots and asterisks, and light turbulence reports (0.05 < EDR < 0.1) are marked with green asterisks. In the RAP configuration in Fig. 7a, large EDR is predicted over eastern Oregon, northern Idaho, and western Montana associated with cyclonic wind shear in the vicinity of the jet stream. However, with smoother topography, weaker turbulence is predicted in this region as shown in Figs. 7c and 7e, even though simulated jets are almost indistinguishable in the left panels in Fig. 7. Based on pilot reports, light EDR is expected over northeastern Oregon, and thus comparing the left panels in Fig. 7, the SMTH4 results in Fig. 7c provide the closest match. P16 showed that the sensitivity of the mountain-wave structure depended on the terrain features, but they did not explain the reason for the strong sensitivity of turbulence in the northern part of the upper-level jet. To further investigate this behavior, additional experiments with the hybrid vertical coordinate are shown in the right-hand panels in Fig. 7. In the HTF simulations, a large reduction of the turbulence in the jet stream region is evident, as shown in Figs. 7b, 7d, and 7f, and this reduction is largely independent of any terrain smoothing. The significant impact of the HTF coordinate in the absence of terrain smoothing (Figs. 7b vs 7a) is intriguing, as the reduced turbulence occurs in a region where terrain influences are relatively weak. This behavior raises questions as to the influence of the hybrid coordinate on the structure of mountain waves forced by the terrain and the role of the upper-level jet in promoting small-scale disturbances.

To address these questions, we evaluate two different vertical cross sections indicated in Fig. 7: one over the high terrain of the Rocky Mountains (Fig. 8), taken along line AB in Fig. 7e, and one along the region in which RAP predicted strong turbulence (Fig. 9), taken along line CD in Fig. 7f. The left-hand panels in Fig. 8 were shown in P16 and demonstrate the strong sensitivity of mountain waves to the shape of topography. We analyzed this same vertical cross section for simulations including the hybrid vertical coordinate, which are displayed in the right-hand panels of Fig. 8. While the amount of terrain smoothing has a significant influence on the mountain waves at smaller scales, the inclusion of the hybrid coordinate has relatively little influence on the mountain-wave structure (cf. the left and right panels for each terrain profile). There are noticeable small-scale waves in both the RAP and RAP-HYBRID for the unfiltered terrain (Figs. 8a,b), and the maximum/minimum magnitudes of vertical velocities (see brackets in each plot) depend much more on the terrain shape than the vertical coordinate. However, it does not mean that HTF is ineffective for upper-level turbulent forecasting in mountainous regions. As shown in the EDR analysis of Fig. 7b, with HTF, light turbulence over Nevada and Utah is reduced from Fig. 7a and is comparable to that in Figs. 7c and 7e using filtered topography.

Along the vertical cross section in Fig. 9, taken along line CD in Fig. 7f, there is no steep terrain, and the shape of the topography is quite similar in the experiments with differing amounts of terrain filtering. However, with the BTF coordinate, there are noticeable differences in the vertical velocity perturbations in the vicinity of the tropopause where the jet stream level winds are strong. In Fig. 9a, the phase of these disturbances do not tilt upwind with height (which would reflect vertical momentum transport), and their horizontal scale is less than 5*dx* (65 km) in RAP. These unphysical waves appear to result from numerical errors in computing the advection of strong winds near the troposphere along coordinate surfaces that retain small-scale structure due to the underlying terrain. As seen in Fig. 9, either terrain smoothing (Figs. 9c,e) or use of the hybrid coordinate (Fig. 9b) can be effective in reducing these errors. For the BTF coordinate, the minimum and maximum vertical velocities (indicated in brackets) become small as the topographic smoothing is increased, while for the HTF coordinate, the amplitudes are small regardless of the terrain smoothing.

### c. Kinetic energy spectra

The kinetic energy spectra for all of the simulations discussed above are summarized in Fig. 10. As described by Skamarock (2004) and Errico (1985), we calculated the one-dimensional spectra, averaging over the 12–24-h integrations at 1-h intervals beginning at 1200 UTC 2 November. The one-dimensional Fourier transforms are calculated along isobaric levels, and results at 500, 300, and 100 hPa are depicted in the top, middle, and bottom panels of Fig. 10, respectively. The left-hand panels in Fig. 10 show spectra from the BTF simulations, and the right-hand panels are for the HTF simulations. Here, we focus on the higher wavenumbers in the simulations, corresponding to length scales less than 500 km. The overall behavior of these spectra is similar to those previously analyzed using WRF (e.g., Skamarock 2004; Waite and Snyder 2013). However, for the simulations presented here, some noticeable differences in the high-wavenumber tails of the spectra are apparent in Fig. 10. The BTF simulation with no terrain smoothing (RAP) shows an unphysical buildup of energy at small scales less than 4*dx* at all isobaric levels. As Skamarock (2004) indicated, energy buildup in small scales can be caused by insufficient diffusion near the grid scale. However, in this case, the artificial energy at small scales appears to be caused by the influences of small-scale terrain structure on the coordinate surfaces.

With filtered topography (SMTH1 and SMTH4), the energy buildup at small scales is effectively removed in both the BTF and HTF simulations, as shown in Fig. 10. For the hybrid coordinate with no terrain smoothing (RAP-HYBRID), while this small-scale energy buildup is apparent at 500 hPa (Fig. 10b), it is not present at higher levels (Figs. 10d,f), where the terrain influences on the coordinate surfaces have been largely removed. This behavior suggests that the artificially large energy at the smallest scales in the absence of terrain smoothing is due to small-scale terrain influences on the coordinate surfaces rather than the direct forcing by small-scale features in the unfiltered terrain. This interpretation is also consistent with the results shown in Fig. 9, where small-scale perturbations are apparently due to numerical errors arising as strong horizontal winds are advected along coordinate surfaces that are perturbed by small-scale terrain influences.

### d. Sensitivity test for

The HTF simulations discussed above were obtained with , representing the level at which the hybrid coordinate transitions to a pure pressure coordinate. For larger (smaller) values of , the level where coordinate surfaces become flat is located at a lower (higher) altitude. With unfiltered topography in RAP, we also conducted simulations with and , and their results are displayed in Fig. 11. In comparison with Fig. 7, for the air turbulence analysis in the top panels of Fig. 11, the main differences are in the region of moderate-intensity forecast turbulence in the vicinity of the line CD. As one might expect, for the larger (smaller) value of , smaller (larger) air turbulent indexes are calculated because terrain influences on the coordinate surfaces are removed at lower (higher) altitudes. This behavior is also apparent in the bottom panels of Fig. 11 in vertical cross sections along line CD. With in Fig. 11c, the magnitudes of artificial vertical velocity perturbations are slightly larger than those in the control simulation with RAP-HYBRID in Fig. 9b, while the perturbations with in Fig. 11d are somewhat smaller. This sensitivity to the value of is not particularly large, and other aspects of the numerical forests should also be evaluated in considering the optimal choice for this parameter.

## 5. Summary

We have formulated a hybrid terrain-following sigma-pressure vertical coordinate for the nonhydrostatic WRF-ARW Model and evaluated it for an idealized test case and for a real-data case for an upper-air turbulence forecast. This hybrid formulation allows the terrain influences on the coordinate surfaces to be removed more rapidly with height than in the basic terrain-following sigma coordinate. The formulation is designed to transition more aggressively toward a pressure coordinate with increasing altitude, to limit the collapse of the vertical grid increments just above high terrain, and to smoothly transition to a pure pressure coordinate at a user-specified level. For the idealized test case, we demonstrate that in BTF coordinate simulations, significant artificial perturbations can arise as the advection of upper-level winds is computed along coordinate surfaces that are perturbed by residual influences of the terrain in the sigma coordinate formulation. This artificial behavior largely disappears with the HTF coordinate, which removes the terrain influences from the coordinate surfaces more rapidly with increasing altitude [as also documented by Schär et al. (2002) and Shaw and Weller (2016)].

In considering the impact of the hybrid coordinate in real-date simulations, we do not attempt here to address the influences on overall forecast accuracy. Rather, we evaluate a specific case for an upper-air turbulence forecast that was previously studied by P16 in assessing the influences of terrain smoothing on the resulting turbulence forecast. As P16 demonstrated, this terrain smoothing can remove the artificial amplification of upper-level disturbances in WRF-ARW simulations that may arise when small scales in the terrain are not filtered. We extend the analyses in that study by demonstrating that the overprediction of upper-level turbulence can also be reduced by using the HTF coordinate, even in the absence of any terrain smoothing. In the EDR turbulence analysis, with RAP configurations, the EDR is overestimated in regions where the upper-level jet is strong and exhibits strong horizontal wind shear. The disturbances that produce these large EDR values in BTF simulations appear to result from numerical errors in representing strong advection along irregular coordinate surfaces that are present due to the influence of small-scale terrain features. With the HTF coordinate, these errors are much reduced as the coordinate surfaces are nearly constant pressure surfaces in the vicinity of the jet. These simulations suggest that either terrain smoothing or a hybrid vertical coordinate may be effective in reducing the overestimation of upper-level waves that contribute to false alarm predictions of air turbulence.

We also conducted simulations (not discussed above) to see if larger internal diffusion is effective in reducing these artificial disturbances that arise with the BTF coordinate. For one test, we switched the calculation of horizontal advection from fifth to third order (increasing the implicit diffusion), and in another, we turned on an explicit sixth-order horizontal filter [as configured and implemented by Knievel et al. (2007)]. Although these filters should increase the dissipation at small scales, they did not noticeably reduce the disturbances contributing to the overprediction of turbulence. The ineffectiveness of these filters may be due in part to the fact that they are computed along coordinate surfaces instead of along constant height or pressure surfaces and therefore may be adversely affected by small-scale terrain influences on the coordinate surfaces.

In evaluating the model kinetic energy, we find that with the BTF coordinate and no terrain smoothing (the RAP case) unphysical energy accumulation is significant at high wavenumbers at all levels. With filtered topography, both the SMTH1 and SMTH4 cases exhibit no artificial accumulation of energy at high wavenumbers. In RAP-HYBRID, since the HTF coordinate transitions more rapidly to constant pressure coordinate surfaces, unphysical buildup of energy is not observed at the upper levels. The reduction in the noise that contributes to the overprediction of turbulence is somewhat sensitive to how rapidly the hybrid coordinate transitions to a pressure coordinate, as determined by the parameter . However, the selection of an optimal value for will also depend on its impact on other important forecast parameters and is beyond the scope of this investigation.

## Acknowledgments

(SHP) This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MEST; 2018R1A2B6008078) and (in part) by the Yonsei University Future-leading Research Initiative of 2018-22-0021. (JBK) Funding for the NCAR portion of this research was provided through support from the National Science Foundation under Cooperative Support Agreement AGS-0856145 and from the NOAA NCAR IDIQ T0006, distributed through the Developmental Testbed Center. NCAR is sponsored by the National Science Foundation.

## REFERENCES

*12th Conf. on Numerical Weather Prediction*, Phoenix, AZ, Amer. Meteor. Soc., 236–239.

## Footnotes

© 2019 American Meteorological Society.

This article is licensed under a Creative Commons Attribution 4.0 license (http://creativecommons.org/licenses/by/4.0/).