## Abstract

Mesoscale atmospheric models are increasingly used for high-resolution (<3 km) simulations to better resolve smaller-scale flow details. Increased resolution is achieved using mesh refinement via grid nesting, a procedure where multiple computational domains are integrated either concurrently or in series. A constraint in the concurrent nesting framework offered by the Weather Research and Forecasting (WRF) Model is that mesh refinement is restricted to the horizontal dimensions. This limitation prevents control of the grid aspect ratio, leading to numerical errors due to poor grid quality and preventing grid optimization. Herein, a procedure permitting vertical nesting for one-way concurrent simulation is developed and validated through idealized cases. The benefits of vertical nesting are demonstrated using both mesoscale and large-eddy simulations (LES). Mesoscale simulations of the Terrain-Induced Rotor Experiment (T-REX) show that vertical grid nesting can alleviate numerical errors due to large aspect ratios on coarse grids, while allowing for higher vertical resolution on fine grids. Furthermore, the coarsening of the parent domain does not result in a significant loss of accuracy on the nested domain. LES of neutral boundary layer flow shows that, by permitting optimal grid aspect ratios on both parent and nested domains, use of vertical nesting yields improved agreement with the theoretical logarithmic velocity profile on both domains. Vertical grid nesting in WRF opens the path forward for multiscale simulations, allowing more accurate simulations spanning a wider range of scales than previously possible.

## 1. Introduction

Advances in high-performance computing have made multiscale atmospheric simulations possible, covering scales ranging from global to large-eddy simulations (LES). The Weather Research and Forecasting (WRF) Model represents one such multiscale simulation framework including efficient parallel computing routines, a suite of physical process parameterizations appropriate for a broad range of scales, and dynamic downscaling capabilities enabled via concurrent grid nesting, a procedure in which multiple computational domains of increasing resolution are integrated simultaneously. With grid nesting, information from the coarse “parent” domain is interpolated and provided as lateral boundary conditions to the fine “child” domain (referred to as one-way nesting). Two-way nesting additionally aggregates information from the fine domain, and feeds the information back onto the coarse domain. While WRF has one- and two-way concurrent grid nesting capabilities, resolution may only be increased in the horizontal dimension, with all domains using a common vertical grid.

It is well known in computational fluid mechanics that grid quality affects the accuracy of numerical solutions (Lee and Tsuei 1992; You et al. 2006), meaning that the lack of flexible gridding in WRF constrains the model’s ability to produce high-quality multiscale simulations. When assessing grid quality, properties such as aspect ratio, orthogonality of coordinate surfaces, and cell volume are commonly considered. Mesoscale models generally use terrain-following coordinates where the vertical grid lines are aligned with the gravity vector, with large aspect ratios (Δ*x*/Δ*z*) near the surface. As multiscale simulations are increasingly used, the high vertical resolution desired on the finest domain yields extremely large aspect ratios on coarser domains. Additionally, a high degree of nonorthogonality (skewness) is introduced in the vicinity of steep terrain slopes. Both grid skewness and aspect ratio have been shown to contribute to numerical errors in atmospheric simulations using the WRF Model (Chow et al. 2008, Lundquist et al. 2010; Mirocha et al. 2013); therefore, it is desirable to reduce these numerical errors by optimizing the grid independently for each domain.

Terrain-following coordinates used by most mesoscale models can be a large source of error for simulations over complex terrain. With terrain-following coordinates, a coordinate transformation is introduced that maps a domain with an irregular lower boundary onto a rectangular grid. While this simplifies the application of lower boundary conditions, it also introduces additional numerical errors (Gal-Chen and Somerville 1975). Inaccuracies from the coordinate transformation are present in each spatially discretized term of the governing equations, and arise from truncation errors due to the coordinate transformation (including grid stretching and skewness) as well as the numerical calculations of the additional metric terms. Numerical errors have been noted in the calculation of the horizontal pressure gradient (Janjić 1977, 1989; Klemp 2011; Zängl 2012), diffusion (Zängl 2003; Zängl et al. 2004), and horizontal advection (Schär et al. 2002) terms in the presence of sloping coordinate surfaces and steep topography.

Mahrer (1984) notes that numerical errors in the calculation of horizontal derivatives occur when large grid aspect ratios are used with sloping coordinates, and states that the minimum vertical grid spacing must be larger than the elevation difference over a grid cell. Mahrer explains that when this condition is not satisfied, accurate horizontal derivatives cannot be calculated based on a traditional numerical stencil using adjacent nodes, and suggests alternative (larger) stencils for the calculation, which would increase accuracy. Most mesoscale models, including WRF, do not use the alternative stencil. Figure 1 illustrates the skewness of computational cells as a function of terrain slope and grid aspect ratio. The gray area on this plot indicates the parameter space which violates Mahrer’s condition that the vertical grid spacing be larger than the elevation change over the horizontal span of the cell. For illustrative purposes, aspect ratios of ½ to 2 are used in this figure; however, a mesoscale model would normally use much larger aspect ratios. For example, a horizontal resolution of 3 km and a vertical resolution of 30 m yields an aspect ratio of 100, and a terrain slope of just 0.6° would violate the condition.

Grid aspect ratio has additionally been shown to affect the accuracy of LES results over flat terrain. Mirocha et al. (2010, 2013) conducted LES of a neutral boundary layer using WRF, and found that the near-surface grid aspect ratio impacted WRF’s ability to accurately reproduce the theoretical logarithmic velocity profile. Nested simulations were constrained by the lack of vertical grid nesting, so that the optimal aspect ratio could be used on either the parent or the child domain, leading to errors on either the parent or the nested domain. Mirocha et al. (2013) hypothesized that vertical grid nesting would improve WRF’s agreement with similarity theory by allowing the optimal grid aspect ratio to be used on each domain.

Many authors have investigated WRF as a multiscale modeling tool (Talbot et al. 2012; Munoz-Esparza et al. 2014; Marjanovic et al. 2014; Mirocha et al. 2014), which could span from meso- to LES scales; however, this is cumbersome without vertical grid nesting. For example, Lundquist et al. (2012) used the WRF Model with 2-m horizontal and 1-m near-surface vertical resolution for LES over Oklahoma City, Oklahoma. Vertical resolution on the order of 1 m generally requires the use of idealized lateral boundary conditions, as it is impractical to use the high-resolution vertical grid on all nested domains. Only with vertical nesting does it become practical to perform multiscale simulations with high-resolution nested domains receiving realistic lateral boundary conditions through nesting from mesoscale simulations and meteorological reanalysis data.

A stand-alone program called ndown processes WRF simulation output and writes initial and boundary condition files to allow for grid nesting of domains integrated in series. This program allows for both increased horizontal and vertical resolution by including interpolation routines for all dimensions, and provides linear interpolation between output times as is done in WRF when forcing from an external file. The vertical interpolation routines in ndown use a high-order monotonic interpolant, and were developed in Moustaoui et al. (2009), where use of vertical nesting improved agreement with observations for simulations near Lompoc, California. Mahalov and Moustaoui (2009), Moustaoui et al. (2010), and Mahalov et al. (2011) coupled WRF to a microscale model using a procedure similar to that in ndown, and found increased vertical resolution was critical to resolving small-scale processes associated with regions of high shear and strong stratification created by mountain lee waves in the Terrain-Induced Rotor Experiment (T-REX). Vertical grid refinement using ndown within the WRF Model was used in Shaffer et al. (2015), and showed an improvement in WRF’s ability to predict the near-surface structure of temperature inversions and low-level jetlike features through comparisons to ground-based observations of temperature and wind speed.

Our vertical nesting method follows from the work of Moustaoui et al. (2009); it utilizes the vertical interpolation method included in ndown, but has been integrated into the WRF framework for concurrent simulation, which has several advantages. First, when ndown is used, a zero gradient boundary condition is applied to vertical velocity, rather than passing vertical velocity between domains. Second, lateral boundary condition updates are limited to the frequency of output on the parent domain, with linear interpolation applied between output times. Both of these limitations make the use of ndown inappropriate for high-resolution simulations and simulations over complex terrain, where high-frequency information and vertical velocity are critical. This is especially true in LES, where it is necessary to pass information on the time scale of turbulent fluctuations between domains. Mahalov and Moustaoui (2009) investigated the treatment of vertical velocity at lateral nested boundaries in their microscale model; they showed reduced errors when passing vertical velocity, as opposed to applying a zero gradient boundary condition, and also noted that a large savings in computational resources would be achieved if the frequent data input/output required was eliminated. Michioka and Chow (2008) found that more frequent lateral updates improved predictions of passive scalar concentration fields when compared to observations of a tracer gas at Mount Tsukuba, Japan. Our method addresses these limitations by passing boundary data to nested domains, including vertical velocity, at each time step, while avoiding computational overhead from writing frequent output files, and facilitating the use of multiple nests. An additional benefit is that while ndown requires vertical resolution be increased with integer refinement (i.e., every *N*th grid level is aligned), our method permits levels to be specified independently on each domain, allowing additional control of the grid.

This work details the implementation, validation, and use of our new one-way vertical grid nesting capability, which enhances WRF’s ability as a multiscale solver. Work presented here uses a modified version of WRFv3.6.1; however, vertical nesting is planned for inclusion in the WRFv3.8.0 release. Details of the WRF solver and the vertical nesting implementation are given in section 2. Vertical nesting is then used in three types of simulations: idealized, mesoscale, and large eddy. The idealized simulations in section 3 validate our implementation of vertical nesting and allow for quantification of errors associated with the nesting procedure. Mesoscale simulations with vertical nesting are presented in section 4 for the T-REX field campaign and comparisons are made to observations. These simulations demonstrate the use of vertical grid nesting in mesoscale mode with a full suite of atmospheric physics. In section 5, LES of a neutral boundary layer is performed to investigate the effects of vertical nesting and grid aspect ratio on WRF’s ability to reproduce the theoretical logarithmic velocity profile. Finally, conclusions and further work are discussed in section 6.

## 2. Implementation of vertical nesting in the WRF solver

### a. Governing equations and treatment at lateral boundaries

WRF is a conservative finite-difference model that solves the nonhydrostatic compressible Navier–Stokes equations (Skamarock et al. 2008). Model equations are cast in an isobaric terrain-following coordinate *η*, which is defined as , where is the dry hydrostatic pressure and is the dry column mass of fluid per unit area. In WRF, the moist Euler equations are transformed into the isobaric terrain-following coordinate, while additional terms such as Coriolis, diffusion, and parameterized physics (represented by *F*) are computed in physical space. A velocity , defined as the contravariant velocity of the vertical coordinate, is introduced in the coordinate transformation, necessitating the solution of an additional equation. Perturbation variables are introduced to reduce numerical errors, and are defined as the deviation from a time-invariant hydrostatically balanced reference state. Pressure *p*, inverse density *α*, geopotential *ϕ* = *gz*, and dry column mass are cast as mean and perturbation values as , where *φ* represents a generic variable and the overbar indicates the hydrostatic base state. The subscripts *d* and *m* represent dry and moist variable states. The transformed equations for conservation of mass and momentum are given in Eqs. (1a)–(1d):

In the above equations the velocity vector is , includes only the horizontal velocity components *u* and *υ*, and operates on coordinate surfaces. Furthermore, a conservation equation, Eq. (2), is solved for each additional scalar quantity, such as potential temperature *θ*, water vapor , and passive scalars:

In addition to solving the prognostic equations above, diagnostic relationships are used to solve for thermodynamic variables. The hydrostatic relationship is used to diagnose perturbations to the dry inverse density , which in the transformed coordinate is given as

Pressure is diagnosed from the equation of state below, where is the ratio of heat capacities of dry air , is a reference pressure, and is the universal gas constant:

Several options for lateral boundary conditions exist in the WRF Model, of which the “specified” and “nested” options are relevant here. Specified boundary conditions are used when boundary conditions are being supplied by an external forcing file, such as from a forecast, analysis, or the program ndown. With specified boundary conditions, variables are temporally interpolated between times in the external forcing file, which are often an hour or more apart. Nested boundary conditions are used on the child domain when multiple domains are integrated simultaneously, and are the option for which we have developed vertical grid nesting. In this case, variables are passed from the parent domain to the child domain, and are temporally interpolated over the parent time step, which is generally on the time scale of seconds. Specified boundary conditions apply to , *θ*, , , and . Nested boundary conditions apply to each variable for which a prognostic equation is solved [Eqs. (1) and (2)], which additionally includes vertical velocity and other scalars ranging from microphysical variables (i.e., moisture constituents such as ice) to passive tracers. Specified and nested boundary conditions utilize a specified zone, where the variable values are directly imposed, and a “relaxation” zone, which uses a forcing term to nudge the solution on the child domain toward the boundary condition value. The width of these zones is run-time configurable; here, the default width of one point for the specified zone and four additional points for the relaxation zone is used.

### b. Vertical interpolation algorithm

The interpolation algorithm used here follows the implementation in ndown (Moustaoui et al. 2009), but has been integrated into the WRF framework for one-way concurrent nesting. The interpolation employs an intermediate vertical coordinate based on log-pressure height, which yields more accurate results than direct use of the *η* coordinate, according to Moustaoui et al. (2009). A cubic monotonic polynomial described in Steffen (1990) is used, in which a piecewise polynomial is constructed that satisfies the variable value and derivative at given interpolation points (Hermite-type interpolation). One known problem of cubic polynomials is producing overshoot–undershoot in interpolation regions between given values. Steffen (1990) solves this problem by adjusting the derivative at interpolation points, if the resulting polynomial will produce local extrema. This method is appropriate for arbitrary datasets, and behaves monotonically on each data interval, by not permitting minima–maxima to occur between known data points. Fields defined at half levels require extrapolation to the model top and surface of the coarse domain, as the new fine grid levels may lie in this region. This is accomplished using the standard Lagrange polynomial, which WRF uses elsewhere for the purpose of extrapolating variables on half levels to the model top and surface.

The interpolation in ndown requires integer refinement, which requires that the number of vertical levels follow Eq. (5), where and are the number of levels on the refined and coarse grids, and is the integer refinement factor:

Integer refinement can be used in the vertical nesting implementation developed here; however, additional capability is developed here to use independent *η* levels on each domain. Independent vertical grid levels may be set by specifying *η* levels for each domain in the namelist, or using default *η* levels as calculated by WRF. When more than two domains are used, the user must choose either integer refinement or independent *η* levels, and cannot use a mix of the two methods. It is possible, however, to vertically refine some domains, while other domains are not vertically refined. Test cases presented here use two domains, and only basic testing has been completed for setups with three or more domains.

### c. Integration of vertical interpolation into the WRF framework for one-way concurrent nesting

When variables are interpolated from a parent to a child domain, data from the parent grid is first passed to an intermediate grid, as part of the parallel communications procedure. The intermediate grid has the vertical resolution of the child, but the horizontal resolution of the parent domain. The vertical interpolation procedure is called for each column of data on the intermediate grid, and then the data are horizontally interpolated onto the nest using the existing routines in WRF.

Our first implementation of vertical nesting followed the method used for concurrent simulations with horizontal nesting, by interpolating the base state and perturbation values from the parent to the child domain at instantiation of the nest, and then passing values at lateral boundaries for all of the prognostic equations during integration. This procedure left the reference state for the child domain out of hydrostatic balance, causing errors, as was also noted in Moustaoui et al. (2009), where reference fields were rebalanced to minimize the transients introduced by these errors. Moustaoui et al. (2009) additionally recalculated from the coordinate definition, which is followed here. Therefore, as part of the vertical grid nesting implementation, rebalancing routines are added to the nest instantiation procedure, and during integration, values for the prognostic variable [Eq. (1d)] are recalculated at lateral boundaries, rather than passed directly from the parent to the child domain. The additional routines added for initialization and integration recalculate the reference state, as well as on the vertically nested grid using the set of Eqs. (6a)–(6d), which are also the equations used at initialization in the WRF Model [see Skamarock et al. (2008) for details]:

### d. Activation of vertical nesting

Vertical grid nesting is planned for inclusion in WRF v3.8.0 for real cases and the ideal LES case. A namelist variable, vert_refine_method, activates vertical nesting and determines which vertical refinement method to use. The default value of 0 indicates no vertical nesting, 1 is used for integer refinement, and 2 is used for arbitrary *η* levels. If option 2 is selected, the user may specify *η* levels in the namelist for each domain using eta_levels. If *η* levels are not specified with option 2, default WRF levels will be calculated. A sample namelist for three domains is shown in Table 1, where vertical levels are set arbitrarily for each domain. The variable e_vert should be set independently for each domain, and feedback (i.e., two-way nesting) must be turned off. If setting arbitrary *η* levels, the variable eta_levels must be a concatenated vector of values ranging from 1 to 0 for each domain. If integer refinement is used, then e_vert should be set so that the number of levels results in the desired nesting ratio according to Eq. (5), and vert_refine_method should be set to 1.

### e. Atmospheric physics

The vertical nesting code was tested with a variety of atmospheric physics parameterizations. Our vertical nesting modifications were found to be compatible with most of the parameterizations without modification. Parameterizations that have been successfully used with vertical nesting without modification include the Yonsei University (YSU) and Mellor–Yamada–Nakanishi–Niino Level 2.5 (MYNN2) planetary boundary layer, thermal diffusion and Noah land surface, WRF single-moment 3- and 5-class microphysics, Dudhia and Goddard shortwave radiation, slope-dependent radiation, topographic shading, Kain–Fritsch cumulus, and Smagorinsky LES schemes [see Skamarock et al. (2008) for references]. Initially, the Rapid Radiative Transfer Model (RRTM), Community Atmosphere Model (CAM), and RRTM for general circulation models (RRTMG) longwave radiation schemes all failed, as they were written with the expectation that concurrently run domains use a common number of vertical grid levels. This limitation is removed from the RRTM scheme by recalculating the number of levels for the scheme (i.e., by passing the number of WRF levels for each domain to RRTM) and adding required supplemental levels for integration to 0 mb. Also, calculation of the number of levels for RRTM is done at each call to the scheme, rather than once at initialization.

## 3. Idealized simulations

### a. Model setup

Idealized nested simulations were carried out over flat terrain to validate the newly implemented vertical nesting method, and quantify errors associated with nesting. The two test cases included here use a two domain setup with periodic lateral boundary conditions on the parent domain, and nested boundary conditions on the child domain. Common soundings of potential temperature and moisture are used for initialization of both cases, however, velocity differs. For all cases, initial conditions are horizontally homogeneous. The first test case is initialized as quiescent and the simulation includes no forcing terms, so that velocities should remain zero and errors are easier to observe. We will refer to this case as the “quiescent case.” The second test case is initialized with a constant velocity profile of *U* = 10 m s^{−1} and includes forcing through a horizontal pressure gradient and surface drag. This case allows for validation with the additional complexity of flow through the lateral boundaries at the parent–child interface, and is referred to as the “forced case.” These cases were chosen to validate our modified nest instantiation and treatment of geopotential, assess the method’s ability to accurately interpolate vertical profiles, and examine the effects of surface forcing where the height of the first grid point, and therefore the reference height used in the similarity theory applied at the first grid level is discontinuous across the nest interface.

Three grid nesting scenarios that differ only in the vertical grid are presented for each case, shown in Fig. 2, which we label: “vertical coarse,” “vertical nested,” and “vertical fine.” Grid parameters are summarized in Table 2. All three grid setups use the same horizontal grid, with Δ*x* = 99 m on the parent domain (d01) and Δ*x* = 33 m on the child domain (d02). Vertical coarse and vertical fine are control cases which do not use vertical grid nesting, and have 40 and 118 vertical levels, respectively. Vertical nested has coarse vertical resolution on d01 (40 levels), and fine vertical resolution on d02 (118 levels). Integer refinement is used in the vertical, so that collocated points exist on the coarse and fine grids, enabling direct comparisons between simulations. The number of vertical levels was chosen to satisfy Eq. (5) with , , and refinement ratio . The coarse grid has m for the first full grid level, and the fine grid has m. Superscript 1 denotes the first grid point above the surface here and throughout. The top of the domain is 4000 m, to allow for a realistic scale for the vertical profiles of potential temperature and moisture, assuming a boundary layer height of ~1000 m.

Test cases are initialized with the potential temperature (*θ*) and specific humidity () profiles shown in Fig. 3 at *t* = 0. Potential temperature is neutral to 1000 m, and then stable above with a lapse rate of 10 K km^{−1}. Similarly, specific humidity is held constant at 5 g kg^{−1} for the first 1000 m, then decreases linearly over 500 m to a value of 0.05 g kg^{−1}, and is held constant at that value to the domain top. These soundings of temperature and moisture were chosen, in part, for the profiles they yielded in geopotential [*ϕ*, determined by Eq. (6d)], which is additionally shown in Fig. 3. In contrast, using a neutral temperature profile in a dry atmosphere leads to zero perturbation geopotential, and, therefore , did not test our treatment of this variable. In tests performed with a neutral, dry atmosphere, errors were so small that they were indistinguishable from errors in simulations with horizontal nesting alone. The cases shown here challenge the vertical interpolant as they include sharp vertical gradients, which are inherently represented less accurately on the coarse vertical grid than the fine, creating a discontinuity across the nest interface. The aim of these tests is to demonstrate that our vertical nesting implementation is accurate and errors introduced at the interface are sufficiently small (i.e., are of the same magnitude as errors introduced by horizontal interpolation alone).

Simulations are run for a total of 48 h with a time step of 1 s on the parent and 1/3 s on the nested domain. This lengthy simulation run time is chosen to rigorously quantify the error growth within the nested domains. All simulations are carried out with a constant eddy viscosity of m^{2} s^{−1}, which contributes to the highly idealized nature of this simulation, and avoids the added complication of allowing the eddy viscosity to change discontinuously across the nest interface. A Rayleigh damping layer is used in the top 500 m of the domain, which damps toward the initial sounding. The surface roughness coefficient is specified as m, and a standard parameterization of the surface stresses , is used, according to the following relation:

where drag coefficient , is the horizontal wind speed, and is a component of the horizontal wind vector.

### b. Results and discussion

Profile comparisons for *θ*, , *U*, and from the quiescent case result in the same conclusions as profile comparisons for the forced case. We, therefore, present and discuss only profiles from the forced case here, so that the development of the forced *U* profile can be examined. Profiles of *θ*, , *U*, and are shown at initialization and after 24 and 48 h of simulation at the center of the nested domain for the forced test case in Fig. 3. Differences between vertical coarse, vertical nested, and vertical fine are nearly imperceptible, indicating that the vertical nesting implementation is correct. With further investigation, however, the effects of vertical nesting become evident. Figure 4 shows *θ*, , *U*, and *p* in the specified zone of d02, along with the collocated profile from d01 for all three sets of simulations at initialization and after 1 h. These profiles focus on regions of strong gradients to show how the solution differs between vertical coarse and vertical fine, and also how the interpolator represents the vertical profile on vertical nested. It can clearly be seen in profiles of *θ* and that the Hermite-type polynomial matches the value and slope at given points, and behaves monotonically on interpolation intervals. Profiles of *U* show that variables defined on half levels are extrapolated onto lower points on the fine grid. Note that while *θ*, , and *U* are interpolated from d01 to d02 at initialization and during integration, pressure is determined from a diagnostic relationship, and, therefore, does not match exactly at collocated points. Also for the vertically nested case, geopotential is diagnosed from the interpolated variables on d02, rather than being interpolated from the coarse domain as in the cases without vertical nesting. We have not included profiles of geopotential in Fig. 4 because differences between solutions on the three vertical grids are nearly imperceptible given the mild gradient of this variable.

Next, errors arising from nesting are quantified for the quiescent case, for which in a perfect simulation with no numerical errors, velocities would remain zero. In this case, *U* and *V* remain identically zero, and vertical velocity component *W* is small throughout d01. While the solution on d01 remains constant, errors arising from both horizontal and vertical nesting appear on d02. Because of the different grid resolution on d01 and d02, the solutions evolve in slightly different ways. This creates small horizontal gradients between the specified zone where the solution is given by d01, and the interior of d02 where the governing equations are solved on the finer grid. Flows result from these differences, and are quantified in Fig. 5, which shows the maximum value of each velocity component within the nested domain as a function of time. Errors due to horizontal grid nesting are seen in all three simulations; additional errors induced by the vertical nesting procedure are included in the maximum error for vertical nested. Peak maximum error values of 0.008 m s^{−1} in the horizontal components of velocity occur in vertical nested 30 min after initialization, and in *W* (0.011 m s^{−1}) 1 h after initialization. The velocity fields subsequently adjust to the vertical nest configuration, likely as the sharp vertical gradients become smoother, but small errors (~ m s^{−1}) continue to be observed over the course of 48 h. Although errors in *U* and *V* in vertical nested are roughly twice those of the simulations without vertical nesting, these errors are still small in magnitude, and do not grow in time. Vertical nesting errors are approximately the same order of magnitude of those introduced by horizontal nesting.

Contours of velocity components for the quiescent case are shown in Fig. 6 at 1 h after initialization, corresponding to the peak errors in *W* observed in Fig. 5. The maximum errors indicated in the time series in Fig. 5 do not appear in Fig. 6, because these values occurred near the corners of the domain, in the relaxation zone where nudging to boundary values from two intersecting lateral boundaries is effectively added together. Small horizontal gradients between the specified zone (along the lateral boundaries) and the interior of d02 cause weak flows, which are strongest at the edges of the relaxation zone, near the height with the greatest change in vertical gradients (~1000 m), evident in *U* and *W* (Figs. 6b,h). These errors appear just after initialization and move downward to the surface over the course of the simulation. Contours of *V* (Figs. 6d–f) show errors an order of magnitude smaller than those in *U* and *W*, and largest in the simulation with the coarse vertical grid. Contours of *V* in a *y–z* slice would appear with errors of the same magnitude as those shown here for *U* in the *x–z* slice. Horizontal “stripes” can be seen in *U* above 1500 m in all simulations. This numerical error is a result of the “ideal” initialization procedure in WRF, which calculates inverse density and pressure via iteration of a set of coupled equations. Following iteration, inverse density is visibly not converged between the vertical intervals of the input sounding. This lack of convergence yields a jagged “shark tooth” pattern that affects the model solution throughout integration.

Figure 7 shows velocity contours as appearing in Fig. 6, but at the end of the simulation period (48 h). At this time, errors are closer in value between the three grid configurations, and remain small. The largest errors in the vertically nested case are seen near the surface, just inside of the relaxation zone. Errors at this location are also present in the cases without vertical nesting, but to a lesser degree. Here, again, the solution is equilibrating between the specified zone, and the domain interior where the prognostic equations are solved.

The forced case has the added complication of flow through the lateral boundaries, and the application of a surface drag model, resulting in errors near the surface at the inflow and outflow boundaries on d02. As seen in Fig. 3, the differences between the three setups are almost imperceptible at the center of the domain, and this is also the case throughout the domain for the duration of the simulation. Figure 8 shows an *x* transect of *U* and *W* at the center of d02 at the end of 48 h. These values are taken from the first collocated point above the surface (for *U*, *k* = 1 for vertical coarse and *k* = 2 for vertical nested and vertical fine, for *W*, *k* = 2 for vertical coarse and *k* = 4 for vertical nested and vertical fine). The vertical coarse *U* value is slower than vertical fine, which also results in a slower velocity at the first grid point of vertical nested after extrapolation down from the coarse grid. The *U* in vertical nested enters the domain with the value of the coarse solution, and though it increases slightly due to the finer vertical discretization, it is forced back toward the coarse solution at the outflow boundary. Values of *W* are consistent with changes in *U* on vertical nested. The surface boundary condition for momentum [defined in Eq. (7)] experiences a discontinuity across the nest interface in vertical nested due to the change in reference height; however, effects of this discontinuity are not obvious with constant eddy viscosity. We will see in the discussion of section 5 that with the LES turbulence closure, in which the eddy viscosity depends on the scales being resolved, the effect of the discontinuity in on the velocity profile becomes pronounced.

Additional simulations not presented herein were carried out with different grids (e.g., domain top extended to 20 km, varying resolutions), physics parameterizations, initialization, and surface boundary condition options (e.g., surface heating, different initial temperature, humidity, velocity profiles). With the results of the simulations presented in this section, along with additional tests, we conclude that the interpolation procedure is working properly and that errors due to vertical interpolation are within an acceptable range.

## 4. Mesoscale simulations

The T-REX field campaign took place in Owens Valley, California, a region of complex, mountainous terrain, between the Sierra Nevada and Inyo mountain ranges, shown in Fig. 9. The intensive observation periods (IOPs) of T-REX focused on observing mountain waves and rotors under strong synoptic forcing, while enhanced observation periods (EOPs) focused on valley flows and boundary layer development under weak synoptic forcing. We selected T-REX for our case study because a large observational dataset, including high-resolution soundings, is available for model validation. Additionally, the complexity of the topography coupled with the strength of the synoptic winds during the IOPs, allow us to investigate model errors over complex terrain and the effects of vertical grid nesting.

### a. Model setup

A set of mesoscale simulations are carried out corresponding to IOP6 (24–26 March 2006), as well as EOP4 and EOP5 (28–30 April 2006) of T-REX. These simulations use a two domain nested setup for four different vertical grid configurations, with details summarized in Table 3. For all configurations, horizontal resolution is Δ*x* = 3 km on d01 and Δ*x* = 1 km on d02, and the model top is approximately 20 km above mean sea level (MSL). We employ the same terminology as in the previous section, where vertical coarse uses 40 vertical levels on both domains without vertical nesting, vertical nested uses 40 levels on the outer domain and 120 levels on the inner domain, and vertical fine uses 120 vertical levels on both domains (no vertical nesting). This relatively large increase in the number of vertical grid points was chosen to provide a greater challenge to the interpolation method, with additional motivation following the work of Mahalov and Moustaoui (2009), Moustaoui et al. (2010), and Mahalov et al. (2011), who reported improvement in comparisons of simulations to observations at even higher vertical grid resolutions. Additionally, *nz* = 120 allows a factor of 3 increase in resolution in both the horizontal and vertical directions for our simulations. Grid stretching is used in the vertical dimension, and when vertical nesting is used, the vertical levels are defined independently on each domain. Additional simulations are performed for EOP4/5 using the stand-alone program ndown with vertical refinement, so that comparisons may be made to simulations using concurrent vertical grid nesting. In the simulations using ndown, 40 vertical levels are used on the larger domain, and 118 vertical levels are used on the subsequently run finer domain. Because ndown is limited to integer refinement, 118 vertical levels had to be used, according to Eq. (5). History is output at 15-min intervals for all simulations, and these output files are processed to create the lateral forcing files for the ndown nested simulations.

Initial and boundary conditions are obtained from the North American Mesoscale Forecast System (NAM) model at 12-km horizontal grid resolution every 6 h. Simulations are run for an 11-h spinup period before the start of the sounding observation comparison period. The MYNN Level 2.5 planetary boundary layer (PBL) scheme is employed, with MM5 surface layer physics and the Noah land surface model. The RRTM and Dudhia radiation schemes are used, with slope-dependent radiation and topographic shading effects included, on all domains. We employ the WRF single-moment 5-class microphysics scheme. Rayleigh damping is applied in the top 5000 m of the domain, with a coefficient of 0.003 s^{−1}, damping toward the base state. References for these schemes can be found in Skamarock et al. (2008). A time step of 12 s is used on d01, and 3 s on d02. In this case, the time step requires greater refinement than is used in the spatial dimension, which is found to be common with vertical grid refinement, especially in the presence of complex terrain where the resolved terrain slope can become steeper on the higher-resolution grid.

### b. Results and discussion

Four nested simulations are performed for the duration of EOP4 and 5: vertical coarse, vertical nested, vertical fine, and a simulation using ndown. Figure 10 shows vertical (*x–z*) slices of *W* at the middle of d02 for all four simulations. In general, the three concurrent simulations appear similar, with the vertical nested solution more closely resembling the vertical fine than the vertical coarse solution due to the increased vertical resolution. Most noticeable are the general differences in structure and magnitude between the concurrent simulations and the simulation produced using ndown. The ndown solution shows lower vertical velocity magnitudes of the mountain waves coming off the Sierra crest, as well as locations where the flow direction is opposite that in our concurrent simulations, namely, the top of the mountain crest between 60 and 70 km on the *x* axis. At this location in the concurrent simulations, there is a small patch of positive *W* surrounded by negative or near-zero values, while the ndown simulation has a wide capping region of relatively strong positive *W* values at the top of the ridge. This is likely due to the different lateral boundary conditions in the ndown simulations, as they are the only significant difference from the concurrent simulations.

These simulations are compared to observational soundings, with data for the comparison extracted from the center of the slices shown in Fig. 10, at 50 km, the center of the domain, corresponding to the sounding release location at Independence, California. Horizontal position data were not available from the EOP soundings taken during T-REX. A total of 30 sounding observation profiles from T-REX during EOP4 and EOP5 are compared to simulated profiles at times within 15 min after the recorded sounding release time. In general comparisons of wind speed and direction, potential temperature, and specific humidity were as expected, based on our previous modeling experience with T-REX (Daniels et al. 2006, 2008; Schmidli et al. 2009). These sounding comparisons, however, indicate that the increase in vertical resolution seems to provide little advantage in accuracy for our setup. Figure 11 shows an example comparison from the 1804 UTC sounding on 29 April 2006, where all simulations give similar results, with the exception of wind direction below ridge crest height. There are exceptions where increased vertical resolution improved the comparisons to observations, such as where abrupt changes in wind direction near ridge crest height are partially captured by both the vertical fine and vertical nested simulations, and completely missed by vertical coarse (not shown). The ndown simulation frequently diverged from the concurrent simulations. Wind direction is the variable with the highest variability between simulations and observations as well as in comparing the sets of simulations to each other, which is not surprising given the complex terrain of the simulation domain. As observed in Fig. 11, vertical nested simulations produced results similar to vertical fine, despite receiving boundary conditions from a simulation with less vertical resolution.

Figure 12 shows height- and time-averaged bias and root-mean-squared error (RMSE), given by Eqs. (8) and (9), for 28 soundings from EOP4/5 (2 of the 30 original soundings had too few observation points to be included in the averages):

Here is the predicted variable from the simulation interpolated to the observed sounding heights and is the observed variable from the sounding observations. The variable *N* is the total number of observations recorded for each sounding, and *n* is an index referring to a particular time and height where the comparison between simulations and observations is made. Special consideration was given to wind direction to calculate bias as the smallest angle between two wind vectors. RMSE is the root-mean-square error. Bias and RMSE are calculated over each sounding and then averaged over all viable soundings (28) taken during the observation period to produce the values shown in Fig. 12. In general, all four simulations show similar values of both RMSE and bias in wind magnitude and direction, potential temperature, and specific humidity, and increased vertical resolution did not yield better comparisons with observations, though simulations using vertical nesting performed similarly to those on the fine grid.

This set of simulations demonstrates that our new method for concurrent vertical nesting is correctly implemented and able to perform predictions similarly to WRF without vertical nesting or WRF using ndown. Although previous researchers have shown improved agreement with increased vertical resolution, this was not seen in the T-REX simulations here. One advantage to using vertical nesting in these mesoscale simulations was to decrease the vertical resolution on outer grids without loss of accuracy on inner grids with finer horizontal resolution. This results in computational savings, but may also be used to reduce numerical errors associated with poor grid quality in simulations over complex terrain. For example, in our simulations of IOP6, which had stronger winds than in EOP4/5, d01 of the vertical fine simulation exhibited numerical errors not seen in the vertical coarse and vertical nested simulations. While vertical coarse and vertical nested run to completion for IOP6 and compare reasonably well with sounding observations, the vertical fine simulation for IOP6 “blows up,” as shown in Fig. 13 where the *U*s and *W*s are overlaid with vertical slices of *V* contours and spurious velocity values are apparent. Decreasing the time step, even to 1 s, did not reliably allow the simulation to run to completion; however, *decreasing* the vertical resolution did allow the simulation to run with results similar to the sounding observations. Chow et al. 2008 and Lundquist et al. (2010) noted in their investigations of errors arising from terrain-following coordinates in WRF that although errors due to the terrain-following coordinate began directly above steep topography, where grid skewness is at a maximum, these errors continued to grow as they advected downstream of the steep topography. Their results may point to why we see large unphysical velocity values downstream of the Sierra ridgeline. Given that the vertical nested and vertical coarse simulations (both with smaller aspect ratios on the parent grid than vertical fine) run to completion while vertical fine becomes numerically unstable, our IOP6 simulation results could be an indication that the ability to control grid aspect ratio with vertical grid nesting, and thus improve grid quality over steep terrain, may not only allow greater numerical stability, but could also save computational resources by allowing fewer grid points and a larger time step.

Additional vertically nested mesoscale simulations, while not described in detail herein, were performed for different grids, locations, and time periods, including the January 2000 snow storm and Hurricane Katrina test cases from the WRF tutorial, and standard testing for inclusion in the WRF distribution, including simulations over five summer and five winter days. All simulations showed expected model behavior, and good agreement was found between vertically nested and nonvertically nested solutions on the inner domain(s) when the number of grid points was the same (equivalent to our vertical nested to vertical fine comparisons).

## 5. Large-eddy simulations

### a. Model setup

LES of neutral boundary layer flow are performed with the goal of examining how grid aspect ratio affects model accuracy in achieving the theoretical solution of a mean logarithmic velocity profile in the lowest region of the boundary layer. Mirocha et al. (2010) and Mirocha et al. (2013) demonstrated that the accuracy of LES in the logarithmic layer was dependent on grid aspect ratio, and found that using the WRF Model, an aspect ratio in the range of 2–4 gave the most accurate results. Following our convention in the previous sections, three nested simulations are performed here, and are referred to as vertical coarse, vertical nested, and vertical fine in reference to the vertical grid resolution. The model setup is summarized in Table 4, where domain size, resolution, and aspect ratio follow Mirocha et al. (2013). Vertical coarse has 46 vertical levels on both domains, vertical nested has 46 vertical levels on the outer domain, and 67 vertical levels on the inner domain, while vertical fine has 67 vertical levels on both domains. All simulations are horizontally nested, using the same horizontal grid, which has Δ*x* = 33 m on the parent domain and Δ*x* = 11 m on the child domain. Two additional stand-alone (i.e., no nesting) simulations with periodic boundary conditions are performed on grids at the same resolution as d02 in the nested setup as control simulations for comparison with the nested simulations. These stand-alone simulations will be referred to as “SA.” For all simulations, a time step of 0.25 s was used on the parent domain, and 1/12 s on the child domain (one-third of the parent time step).

In the nested setup, d01 uses periodic lateral boundary conditions, while d02 uses nested boundary conditions. Vertical grid spacing is prescribed similarly to Mirocha et al. (2013), with 5% grid stretching up to the model top of 1400 m. Vertical grid spacing at the first level is m on the coarse grid, and m on the fine vertical grid. Rayleigh damping is applied in the top 300 m of each domain, with a coefficient of 0.003 s^{−1}, damping toward the initial sounding. Standard parameterization of the surface stresses , *i* = 1, 2 is used, as defined in Eq. (7). The simulations are forced with a geostrophic wind of m s^{−1} and have a Coriolis parameter *f* = 0.0001 s^{−1} (≈45°N), as in many LES investigations (e.g., Andren et al. 1994; Chow et al. 2005; Mirocha et al. 2013). Simulations are run for 24 h of spinup to allow damping of inertial oscillations. Time averaging in postprocessing was performed over 4 h following the spinup period for each simulation. The simulations presented here use the standard Smagorinsky turbulence closure model included in WRF. The Smagorinsky model is known to have difficulty in accurately reproducing the log-law (e.g., Andren et al. 1994; Porté-Agel et al. 2000; Chow et al. 2005) compared to more sophisticated models; however, its simplicity and known sensitivity to grid aspect ratio make it a desirable choice for this study.

### b. Results and discussion

Figure 14 shows instantaneous contours of *U* velocity after 28 h of simulation (i.e., at the end of the averaging time) for d01 and d02 of the nested simulations, along with the stand-alone simulations. The effect of grid nesting can be seen in d02 (Figs. 14c–e), where large-scale features from d01 pass into d02 and persist to the outlet, though small-scale features do appear within these larger features in the second half of d02 (after about 2000 m), also observed by Mirocha et al. (2013). In contrast, the SA simulations (Figs. 14f,g) display small-scale features throughout the horizontal and vertical extents of the domain. The snapshots in Fig. 14 show qualitatively that the SA simulations exhibit more intricate structure than any of the nested simulations; all of the nested cases exhibit development of some smaller-scale features in addition to the main structures passed from the parent domain, and vertical nesting produces similar results to using both vertical and horizontal nesting.

Spatially and temporally averaged surface stresses are shown in Fig. 15 along the *x* direction for all simulations. Here it can be seen that SA simulations and d01 with common vertical grids have different values, despite both having periodic boundary conditions, due to the different horizontal resolutions (Δ*x* = 33 on d01 and (Δ*x* = 11 for the SA simulations). For d02 of vertical coarse and vertical fine, takes the value from the outer domain at the inlet, dips, recovers, and then dips and recovers again just before the outlet. This behavior is consistent with results in Mirocha et al. (2013), where it was shown that smaller scales responsible for the vertical transport of momentum take time to develop as the flow enters and moves across the nested domain, so that downward transport of momentum is reduced in this developing “transition zone,” which creates a velocity deficit near the surface of the nested domain, and thus reduced surface stresses. Vertical nested follows this pattern of development along the streamwise direction; however, is much higher at the inlet and outlet than on the parent domain. This is because the velocity used to calculate on d02 at the inflow and outflow boundaries must be extrapolated from the first point above the surface on the coarse vertical grid of d01, and this extrapolated value ends up being larger (~1 m s^{−1}) than the value on d02 of vertical fine. The SA simulations represent an ideal value for d02 of the nested simulations, which the nested simulations should be “adjusting” to. The spike at the outlet of d02 for the nested simulations is due to readjustment leading up to and through the relaxation zone to the conditions on d01.

Figure 16 shows a comparison of each of the three nested simulations on their respective inner (d02) and outer (d01) domains, along with the SA simulations performed on the d02 grids. Figures 16a,c,e show mean velocity profiles along the streamwise direction, while Figs. 16b,d,f show the same profiles normalized by versus normalized height. On the inner domain, spatial averages are performed at individual *I* indices, and averaged over (30 < *j* < 210) to capture the developing turbulent flow as it progresses across the domain. On the outer domain, flow is fully developed throughout with periodic boundary conditions, so spatial averages are performed over the entire domain. The time average is performed as a postprocessing step over 4 h following the 24-h spinup, using data saved at 1-min history intervals.

Figures 16a,c,e show that the mean horizontal wind velocity on the inner domain is heavily influenced by the lateral boundary conditions coming from the outer domain, to a distance past 20% of the length of the inner domain (*i* = 80). Profiles from the inner and outer domains are nearly identical in this region near the inflow boundary. In all cases, as the fully developed flow passes *i* = 320 (within 15% of the outflow boundary), the mean wind speed profiles closely match the SA simulations up to approximately 20 m above the surface. Above this height, the flow adjusts to match the lateral boundary conditions (mean wind speed profile of the outer domain). Horizontal nesting alone accounts for significant divergence of the inner-domain solution from the outer-domain solution and from the SA solution (Figs. 16a,e). At just below 120 m (*k* = 23 and *k* = 11 on the fine and coarse vertical grids, respectively) all d01 and d02 vertical coarse and vertical fine profiles differ from their respective corresponding SA simulations by 0.1 m s^{−1}, while the vertical nested simulation differs from SA fine by ~0.2 m s^{−1}, approximately double the difference in the other two simulations; however, these differences decrease with height above 120 m in vertical nested and increase with height above 120 m in vertical coarse and vertical fine (not shown).

These observations can also be made in Figs. 16b,d,f where comparison is made to the log law. In this column of profiles, *H* is the boundary layer height of ~1000 m, the approximate height above which stresses are attenuated. According to Mirocha et al. (2013), the best comparisons with the log law, using a roughness length of 0.1 m, are achieved when the grid aspect ratio is AR ≈ 4. In their study, it was only possible to optimize AR on either the inner or the outer domain, while with vertical nesting, we are able to control AR on each domain independently. Thus, in Fig. 16b the outer domain (d01) is optimized with AR ≈ 4, while the inner domain (d02) has AR ≈ 1.3; d01 shows the optimal comparison with the log law, while the lines corresponding to vertical coarse d02 and SA coarse (also on d02) are farther from the theoretical log law line. We see that in Fig. 16f, d02 aligns with the log law, having AR ≈ 4, while d01 with AR ≈ 11.5 does not compare as well. Mirocha et al. (2013) hypothesized that if AR ≈ 4 could be achieved on both domains, then good agreement with the log law could also be achieved on both domains; our results shown in Fig. 16d demonstrate that this is true in our simulations. The trade-off is slightly slower equilibration within the nested domain due to inlet boundary conditions from a domain with a coarser vertical grid, and, therefore, fewer resolved small-scale structures.

## 6. Summary and conclusions

This paper describes our implementation of concurrent vertical grid nesting in WRF based on the program ndown created by Moustaoui et al. (2009). Vertical nesting is now fully integrated into the nesting framework of WRF so that lateral boundary condition updates for nested domains can take place at every parent grid time step. Thus vertically nested simulations can now be run concurrently.

Results of idealized simulations and rigorous error analysis validate the implementation and indicate that errors resulting from vertical nesting are bounded (i.e., reach small, constant values by 48 h of simulation time). Mesoscale simulations of three case study days from T-REX show no measurable change in accuracy when employing vertical nesting (cf. WRF without vertical nesting), based on comparisons to sounding observations. Large-eddy simulation results validate the hypothesis of Mirocha et al. (2013), that the ability to control grid aspect ratio through vertical nesting allows better agreement of the log law on both the outer and the nested domains.

Our results indicate that it may be possible to decrease the vertical grid resolution on the outer domain without sacrificing accuracy on the inner, nested domain. This means that grid aspect ratio adjustment is possible in regions of complex terrain, so that errors over steep slopes may be reduced, and overall simulation accuracy may be increased even when vertical grid resolution is decreased. This result could also mean computational savings for mesoscale weather forecasting and for cutting-edge research in multiscale modeling.

In the interest of providing useful guidance and opening the way forward for more modeling studies with vertical nesting, we have intentionally exercised the model in ways that were likely to expose problems and errors. Vertical nesting is one step toward a robust, general multiscale modeling framework, but there remains an array of possible improvements and open questions to be explored. Ongoing and future work, therefore, includes the following: exploring effects of using more sophisticated turbulence closure models with vertical nesting, improving the surface boundary condition by modifying the extrapolation function used near the surface with the goal of eliminating the discontinuity in surface stress across the lateral boundary, further exploring errors due to grid aspect ratio especially over complex terrain, implementation of vertical nesting with two-way nesting, and modifications to use vertical nesting with additional radiation models and physical parameterizations.

## Acknowledgments

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344, and was supported by the U.S. DOE Office of Energy Efficiency and Renewable Energy and the LLNL Laboratory Directed Research and Development program as Project 14-ERD-024.

## REFERENCES

*13th Conf. on Mountain Meteorology*, Whistler, BC, Amer. Meteor. Soc., 9A.5. [Available online at https://ams.confex.com/ams/13MontMet17AP/webprogram/Paper141221.html.]

*12th Conf. on Mountain Meteorology*, Santa Fe, NM, Amer. Meteor. Soc., 7.2. [Available online at https://ams.confex.com/ams/SantaFe2006/techprogram/paper_114757.htm.]

*18th Symp. on Boundary Layers and Turbulence*, Stockholm, Sweden, Amer. Meteor. Soc., 11A.1. [Available online at https://ams.confex.com/ams/18BLT/techprogram/paper_139961.htm.]

*14th Conf. on Mountain Meteorology*, Lake Tahoe, CA, Amer. Meteor. Soc., 10.1. [Available online at https://ams.confex.com/ams/14MountMet/techprogram/paper_173801.htm.]

*10th WRF Users’ Workshop*, Boulder, CO, National Center for Atmospheric Research, 7 pp. [Available online at http://www2.mmm.ucar.edu/wrf/users/workshops/WS2009/abstracts/6-03.pdf.]

_{3}and CO induced by mountain waves in the upper troposphere and lower stratosphere during terrain-induced rotor experiment

## Footnotes

This article is included in the Terrain-Induced Rotor Experiment (T-Rex) special collection.