Crossing Multiple Gray Zones in the Transition from Mesoscale to Microscale Simulation over Complex Terrain

This review paper explores the field of mesoscale to microscale modeling over complex terrain as it traverses multiple so-called gray zones. In an attempt to bridge the gap between previous large-scale and small-scale modeling efforts, atmospheric simulations are being run at an unprecedented range of resolutions. The gray zone is the range of grid resolutions where particular features are neither subgrid nor fully resolved, but rather are partially resolved. The definition of a gray zone depends strongly on the feature being represented and its relationship to the model resolution. This paper explores three gray zones relevant to simulations over complex terrain: turbulence, convection, and topography. Taken together, these may be referred to as the gray continuum. The focus is on horizontal grid resolutions from ∼10 km to ∼10 m. In each case, the challenges are presented together with recent progress in the literature. A common theme is to address cross-scale interaction and scale-awareness in parameterization schemes. How numerical models are designed to cross these gray zones is critical to complex terrain applications in numerical weather prediction, wind resource forecasting, and regional climate modeling, among others.


Introduction
Atmospheric models are becoming an indispensable source of knowledge about our planet's current and future state.In addition to predicting daily changes in weather, numerical models are used to drive air quality models, to provide wind resource assessments, and to predict long-term climate such as spatial and temporal changes in weather patterns.These forecasts play a critical role in everyday life, affecting economic, personal, and environmental decisions among others [1].
Most of these atmospheric predictions are made over complex terrain, where topography and land-use heterogeneity directly impact the evolution of the atmosphere.The influence of complex terrain can drive local circulations (e.g., thermally-driven winds [2]) or regional scale flow (e.g., upstream blocking from mountain ranges [3]).Lehner and Rotach [4] reviewed many of these mountain transport and exchange processes, including gravity waves and moist convection.
For numerical models, complex terrain poses a number of challenges.First, accurately representing turbulent processes is essential to models at all scales, to properly dissipate kinetic energy and represent mixing in the atmosphere.The influence of complex terrain on turbulence parameterizations must also be considered.Second, the parameterization of convective processes is required at coarser resolutions, and is impacted by orographic triggering of convection over complex terrain.This is difficult to capture because local moisture and temperature fields are strongly modulated by mountain circulations.Finally, the accurate representation of topography is critical, yet is limited over steep terrain by choices in the underlying numerical grid representation.
In all of these areas, there are unique challenges posed by the transition from mesoscale to microscale modeling regimes and the gray zones encountered in the process of moving across resolutions.The gray zone here refers to the challenge of transitioning from parameterization of physical processes to explicit resolution of these processes as the grid is refined [5,6].(In the context of turbulence modeling, the gray zone is often called the terra incognita, a phrase introduced by Wyngaard [6].)There is of course strong motivation for atmospheric models to move to higher and higher grid resolutions to more accurately represent physical flow processes.In particular, the representation of topography improves dramatically at higher resolutions.Similarly, more detailed land-use forcing can improve the representation of vertical exchange processes.What is unclear, however, is how to properly transition from the mesoscale, where much of these processes are parameterized, to the microscale, where a large fraction are resolved.
For example, at 1 km resolution, the topography representation over very complex terrain like the European Alps begins to improve greatly, so that valley-mountain circulations and moist convection can be adequately resolved.At 100 m resolution, local circulations in narrow valleys and moisture flux gradients between slopes and the valley can be captured [7].At 10 m resolution, the wake of a wind turbine can be resolved [8], but the resolved topography becomes quite steep, meaning that new numerical techniques, such as alternative grids, are required [9,10].There are then decisions to be made about the transition between different numerical and gridding techniques.Similarly, at 1 km horizontal resolution, deep convection can be resolved, meaning that convective parameterizations are no longer needed [11].However, as we move to finer scales, we find that the convective cell size formed still depends on the choice of turbulence closure and may be unphysically large [12,13].This also indicates a need to study in detail the transition from coarse to fine scales.
Experience has shown that, as the resolution is refined into the gray zone, there can be violations of some underlying assumptions used in many physical parameterization schemes, for example, turbulence equilibrium, isotropic eddies, surface similarity theories, and local down-gradient flux closure assumptions [6].In addition there are other basic simplifications which come into question, e.g., the use of 1D radiation and 1D soil column models, static vegetation models, and the choice of anelastic or hydrostatic approximations to the governing equations.At the same time, many new research studies both in weather and regional climate modeling are performing and interpreting results at 1 km resolution e.g., [14], which is in the middle of at least one gray zone Bryan et al. [15].Are the simulation results accurate?What should be done for example about turbulence modeling at these intermediate scales?The challenge becomes to properly traverse the gray zones and apply current physical parameterizations at scales other than those originally intended, or to develop new methods suitable for the gray zone.
Currently, there is a tremendous thrust to refine horizontal resolution of operational and research atmospheric models (see, e.g., Figure 6 from [4]).On the one hand, this trend is driven by the recognition that higher resolution will refine the representation of crucial atmospheric processes including convection, topography and turbulence.On the other hand, it is motivated by recent advances in computational power, which includes not only conventional CPU-based hardware architectures, but increasingly heterogeneous many-core architectures with accelerators and graphics processing units (GPUs).Such hardware has already been used in recent applications of numerical weather prediction models [16], idealized atmospheric models [17,18], and regional climate models [19].
Studies suggest that within a decade there will be kilometer-resolution global coupled climate models that will further enhance the interest in the gray zone [20].
This review paper focuses on the gray zone as it applies to simulations over complex terrain.Our goal is to address issues that arise as we move across scales to finer resolution, specifically: representation of turbulence, convection, and topography.There are certainly other processes that face their own gray zone transitions, e.g., orographic gravity waves [21], which are not discussed here.It can also be argued that even very flat terrain can have complexity which directly affects the flow, e.g., the CASES-99 experiment in Kansas (see [22]), hence the gray zone discussions herein are sometimes general and hence broadly applicable beyond mountainous regions.The model resolutions of interest here (roughly 10 km to 10 m) are relevant for a large range of models, from global models that are now being run at 1 km resolution to 10 m simulations explicitly representing flow in urban areas or wind turbines and their wakes.Here, we specifically consider this range of scales from the perspective of mesoscale and microscale atmospheric modeling.

Scales of Interest and Governing Equations
The governing equations for all atmospheric scales of motion are given by the Navier-Stokes equations combined with the continuity equation, to solve for the velocity field, u i , and pressure, p, along with transport equations for potential temperature, θ, and microphysics variables.At the mesoscale, Reynolds averaging is used to separate mean and fluctuating quantities in the Reynolds-averaged Navier-Stokes (RANS) approach.At the microscale, large-eddy simulation (LES) is often used, where a spatial filter separates large energy-containing eddies from smaller, unresolved motions.In both RANS and LES frameworks, the Navier-Stokes equations are filtered, with a time or space filter, respectively, and additional equations are used to represent the effects of unresolved turbulent fluxes.We illustrate this here for the momentum equation, where the effect of the filter is represented using an overbar: where ρ is density, and Coriolis and viscous terms have been neglected for brevity.For RANS, the effect of turbulent motions is contained in the Reynolds stress term, τ ij = u i u j , where u i is the fluctuating velocity field.For LES, τ ij = u i u j − u i u j and is defined as the subfilter-scale stress.In many mesoscale models, only the vertical components of the turbulent transport are accounted for, i.e., only the terms τ 1,3 and τ 2,3 are retained.This approximation is sometimes referred to as 1D turbulence.Zhong and Chow [23] discussed the emerging need for 3D PBL schemes over complex terrain, as also evidenced by the recent work of B. Goger [24].Mesoscale models typically use eddy-viscosity turbulence closures to represent all the turbulent motions.The grid spacing ∆ is larger than the length scales L of interest.L is typically on the order of the boundary layer depth (1 km), so with ∆ ∼ 10 km for typical mesoscale models, none of the turbulence is resolved.The success of mesoscale RANS models has relied on the premise that atmospheric turbulence can be ensemble averaged, meaning that grid cells must be large enough to contain a large statistical sample of eddies [6], as illustrated in the schematic in Figure 1a.
In contrast, large-eddy simulation (LES) is used at grid resolutions of about 10-100 m, with the goal of resolving most of the energy-containing eddies, as shown in Figure 1b.In this case, ∆ << L and most of the turbulence is resolved.Traditionally, LES models have been applied in idealized studies with small domain sizes because of computational constraints.
Advances in computational speed and availability are pushing mesoscale and LES modeling regimes to new limits.Atmospheric mesoscale models are now used across a vast range of spatial scales, with horizontal grid resolutions ranging from ∼ 50 km to < 1 km (e.g., [7,[25][26][27][28][29]).Grid refinement is used to span this large resolution range, and simulations with resolution at intermediate scales enter the terra incognita or gray zone [6,30].In the gray zone, ∆ ∼ L, and the turbulence is only partially resolved, as shown in Figure 1c.Furthermore, the assumptions for turbulence closure become problematic; with partially-resolved eddy motions, it is unclear how to best represent turbulent motions.At kilometer-scale resolutions, for example, some convective processes become partially resolved (see Section 5).Similarly, turbulence in the atmospheric boundary layer becomes partially resolved at intermediate resolutions (see Section 4).For convective boundary layer (CBL) conditions, this might happen at about 500 m resolution.For stable boundary layer (SBL) conditions, this might happen at 50 m resolution [22].Thus, the gray zone can be defined by particular flow characteristics and can vary over diurnal time and space scales.Reconciling the different approaches to turbulence modeling may smooth the path towards simulations that apply at coarse, fine, and even intermediate resolutions, within the terra incognita.

Transitioning across Scales: Grid Refinement
Numerical weather prediction has long relied on grid refinement to focus on the domain of interest at higher resolution, while reducing computational requirements compared to uniformly refining the grid throughout the domain.Refinement in regional scale models is typically accomplished through grid nesting, i.e., a series of telescoping grids of increasing resolution using an integer refinement factor (see, e.g., [23]).Global scale models often use gradual grid refinement over a region of interest (see, e.g., [31]).Zhong and Chow [23] provided a comprehensive list of a dozen mesoscale models (see their Table 10.1), which all use grid nesting.Given this ubiquitous grid nesting framework, we focus here on this traditional nesting process.Another related issue is the treatment of surface land-use data and associated land-surface models at different resolutions, as well as selection of appropriate top and bottom boundary conditions (these are not discussed here, but the interested reader can refer to Zhong and Chow [23]).Challenges at grid refinement zones are common to both global models with regional refinement and to nested grid strategies.
To predict weather at local or regional scales with a grid nesting framework, multiple nests are used to downscale information from the global scale, through the mesoscale to the microscale, and, optionally, information from the fine-scale domain is up-scaled onto coarser domains.For regional climate modeling, this process is termed dynamical downscaling, which uses a global climate model to provide boundary conditions for a regional climate model, thus providing physical mechanisms for deriving a regional prediction of future climate [32,33].Integer nesting ratios are required and ratio values of 1:3 or 1:5 are typical, with odd numbers leading to slightly more accurate interpolation.Several groups have recently experimented with much larger ratios, e.g., 1:11 or 1:15 [22,34,35].In this way, it is possible to transition from the mesoscale to the microscale, as shown in the schematic in Figure 2. Changing the grid nest ratio provides an opportunity to explore the grid nesting configuration as an avenue towards optimizing the mesoscale to microscale transition across the gray zone.
Depending upon the application, one-way or two-way nesting approaches may be employed.In the former case information is communicated only from the outer to the inner nest, and the simulations may be considered independently.In the latter case, information is allowed to propagate from the inner to the outer nest, thereby enhancing consistency between the two nests, though the benefits of two-way nesting are still being studied [36,37].Relaxation lateral-boundary conditions are employed for nested domains (e.g., [38]).Relaxation for instance serves to facilitate a smooth transition from the solution on the coarse domain to the finer-scale solution on the nested domain.Over complex terrain, additional issues exist with grid nesting regarding placement of grid boundaries over steep topography [39].As the resolution becomes finer, these nested simulations transition from the mesoscale into the terra incognita, and eventually to the LES scale.It is in this nesting process where it becomes unclear how to represent subgrid motions at different resolutions.At some point, the resolution becomes fine enough that the turbulence is explicitly resolved and does not need to be parameterized.However, most of the commonly used closures do not take into account any grid dependence and thus incorrectly continue to use the same parameterizations at all scales, resulting in so-called "double counting" of the turbulence effects.

Challenges in the Gray Zone
Comparing dynamics of the dry convective atmospheric boundary layer to the classic Rayleigh-Bénard thermal instability theory, Zhou et al. [13] laid out the difficulties of the gray zone in terms of flow dynamics.The gray zone, as discussed above, refers to the range of model grid spacing (∆) that is comparable to the length scale of interest in the flow (L).In the convective boundary layer (CBL) gray zone, erroneous flow predictions arise in terms of the onset of convection, the horizontal length scale of the convective cells, and higher order statistics such as variance and skewness.Figure 3 shows contours of vertical velocity in the convective boundary layer as simulated at LES scales at 25 m resolution, and, in the gray zone, at 1200 m resolution.The CBL is simulated here over idealized flat terrain using periodic boundary conditions.Heating at the surface is prescribed and drives the evolution of the boundary layer.The scale of the convective cells in the two simulations is strikingly different; with LES, the convective cell is about the same size as the boundary layer depth, z i , whereas with the gray zone, the convective cell is 6-8z i .Figure 3 also shows the LES results filtered at 1200 m, which can be considered the filtered true solution, showing a very different spatial pattern than what the gray zone simulation produces.Another striking feature of the gray zone simulation is the grid-dependence of the onset of resolved convection-convection begins earlier in the afternoon at 800 m than at 1000 m, than 1200 m, etc., as can be seen in Figure 14 of Zhou et al. [13], which shows resolved velocity variances as a function of time.Applications of model results in the gray zone could lead to erroneous conclusions, e.g., about scalar transport or moist convection, clouds, and precipitation patterns, since these rely on the structure of the CBL convection.This model error in the gray zone is rooted in the model's inability to resolve the wavelength of the most unstable convective eddies, because of grid resolution limitations.It also has to do with the turbulence closure, and we can see that this plays a role too in [40], where the TKE-1.5 and Smagorinsky closures show differences of 2 h in the onset of convection at gray-zone resolutions.Earlier onset of moist convection in higher-resolved simulations was also noted by Langhans et al. [41].A turbulence scheme that does not take into account the horizontal grid spacing is not likely to address the CBL gray zone problem, as discussed in Section 4 below.
To complicate matters further, errors on gray zone grids can be propagated down to finer nested grids.Following the discussion by Zhou et al. [13], we consider a 400-m horizontal resolution simulation one-way nested within a gray zone simulation of 1200 m resolution.This is compared to results from a stand-alone 400-m run with periodic lateral boundaries.Figure 4 shows horizontal contours of vertical velocity, where large square cells are visible near the south and east inflow boundaries of the nested run (Figure 4b).The size of the structures slowly decreases as the flow progresses from right to left across the domain, eventually matching the size of the structures in the 400-m stand-alone domain (Figure 4a).In cases where the boundary layer wind speed is strong, the preferred convection patterns are rolls instead of cells [42].Grid-dependent convection rolls formed on a parent gray zone grid often extend all the way across its inner nest (see, e.g., [43,44]).One approach to address these gray zone challenges is to invent new turbulence closures that are "scale-aware", and examples are given in the following sections.This, combined with techniques to enable a more efficient transition across nest boundaries can enable more accurate representation of turbulent structures in grids nested within the gray zone.

Considerations at Lateral Boundaries
Higher resolution enables the representation of increasingly smaller dynamical flow structures.Depending upon the resolution and weather situation considered, this may concern different mesoscale or microscale phenomena, such as the representation of narrow cold fronts and their associated rainbands, the generation of thunderstorms, or the formation of turbulent boundary-layer eddies.Each of these phenomena has its own characteristic life cycle, and the generation of these structures thus requires time and space for their formation.For instance, deep convective cells have a life time of 3-6 h [45], and they will thus typically travel 50-100 km before reaching maturity.The nesting strategy needs to account for these temporal and spatial scales.Several studies have systematically tested different approaches, addressing-for each of the meshes considered-domain size, resolution and up-dating frequency at the lateral boundaries.For instance, Brisson et al. [46] systematically conducted and validated a series of simulations with a target resolution of 2.8 km, driven by reanalysis data of 75 km resolution.They advocate a two-step nesting strategy with resolutions of 25 and 2.8 km, using boundary-updating frequencies of 6 and 1 h, respectively.The refinement is thus not only in space but also in time, which is critical as the synoptic-scale data may only be available every 6 h.Addressing the statistics of daily precipitation in a convective episode, Brisson et al. [46] also found that a minimum distance of about 150 km was needed between the lateral boundaries and the target region.
Nesting from the mesoscale to the microscale creates challenges in passing turbulent flow information across domains.At the mesoscale, turbulent motions are largely parameterized and the resolved flow field is relatively smooth.When this flow field is nested down to finer resolution, smaller scales are resolvable on the fine grid, yet the development of these structures is very slow.Indeed, over flat terrain, finer-scale turbulent structures can take more than half of the domain length to develop and are often still too large [47,48].This means that the benefits of higher resolution on the representation of the turbulent flow are often not fully realized, or require large computational domains to enable the build-up of the appropriate eddy statistics.Mirocha et al. [47] and Goodfriend et al. [48] showed that the choice of turbulence closure model greatly affects the transition of turbulent scales at a grid refinement interface.Dynamic models and mixed models which include scale-similarity terms performed the best.
A series of recent papers has examined the use of perturbations near the inflow domain boundaries on the nested grid.While perturbations near the Nyquist limit of the grid were found to quickly dissipate and hence be ineffective at triggering the appropriate turbulence, sinusoidal perturbations applied at 8∆x wavelengths are successful in enhancing turbulent length scales [43,49].This so-called cell perturbation method (CPM) applies the perturbations to the potential temperature field in three rows of 8∆x width around the lateral boundaries, extending vertically to a specified height within the PBL.Munoz-Esparza et al. [43] and Munoz-Esparza and Kosovic [50] tested the method for a range of large-scale forcing conditions to find appropriate parameters for applying the perturbations.
The CPM may prove highly useful in the gray zone, not just for transitioning from the mesoscale, where no turbulence is resolved, to an LES scale, but also when the domain has under-resolved turbulence as discussed in [13,49].Mazzaro et al. [51] investigated the impact of smooth lateral boundary conditions from parent mesoscale domains at gray zone resolutions on nested LES using a set of idealized simulations.They found that the LES nest performance improves with increasing mesoscale parent resolutions, and that the CPM helps to break down the large scale structures.Munoz-Esparza et al. [34] used CPM to reduce the fetch required to fill out the finer scales in a real case simulation with the CPM applied on the LES domain, nested within a mesoscale domain.Van Veen [52] considered a complex terrain case with a nested setup over the Perdigão valley, the site of the 2017 international field campaign in Portugal.Using the Weather and Research Forecasting (WRF) model, the authors found that, even over complex terrain, use of the perturbation method helps break down spurious elongated flow structures and enhances the turbulence spectra of the flow field at finer scales.Figure 5 shows snapshots of the flow field with and without the cell perturbation method on a nested domain at 150 m horizontal resolution (previously shown in Figure 2).The flow is from the southwest, and the turbulent structures are clearly modified when the CPM is used (Figure 5b) compared to when it is not (Figure 5c).These are compared with simulations from a much larger D03 domain (Figure 5a), which shows the extended fetch needed for similar turbulent structures to develop without CPM.Energy spectra taken along transects normal to the mean flow direction (Figure 6) indicate an increase in energy at the higher wavenumbers due to the CPM which agree well with results from the extended fetch of the large domain.This demonstrates that the CPM effectively decreases the fetch required to develop appropriate fine scale turbulence on the nested domain.2c.The diagonal lines are used for the w-spectra comparison in Figure 6.
Another transition at lateral boundaries is in the grid aspect ratio.Often the vertical resolution remains constant to facilitate interpolation at nest boundaries, yet the horizontal resolution changes by factors of 3-10 or so.This results in large changes in grid aspect ratio, which enable some improved representation of vertical gradients, but which can create errors that affect the representation of turbulent structures on the nested grids.Zhong and Chow [23] discussed cases where coordinate cells become highly skewed as resolution is refined and provide a rule of thumb for selecting grid aspect ratios.One way to alleviate issues with large grid aspect ratios is to change the refinement in the vertical on different nest levels.Daniels et al. [53] implemented vertical nesting capabilities into WRF so that the vertical grid could be adjusted on each nest level.They validated the vertical grid refinement feature for both mesoscale and large-eddy simulations, and showed improved model performance when vertical grids are purposefully chosen for each nest.Mirocha and Lundquist [54] examined the effects of vertical grid refinement on large-eddy simulations with two nested domains and found that although the vertical grid on the outer domain was coarsened, leading to coarser turbulent structures as the inflow to the finer domain, effects on the inner solution domain were minimal.

Alternatives Which "Skip" the Gray Zone
To avoid the challenges created by the spurious structures generated by grid nesting in the gray zone, one approach is to get rid of the convection that occurs on gray zone grids by nudging the gray zone flow towards a coarser mesoscale solution, as presented in [55].The mechanism behind this nudging is to increase horizontal eddy diffusivity, to prevent resolved convection from occurring on the gray zone grid.This effectively removes the erroneous flow features, but sacrifices the additional computational cost of running the model at a higher spatial resolution than that of a mesoscale grid.With a similar goal to counteract the emergence of under-resolved convection from the mesoscale domain, Munoz-Esparza et al. [34] increased the Smagorinsky constant c s on the highest resolution mesoscale domain.
As mentioned above, some global scale codes use adaptive grids which can locally be refined about the region of interest, as for example in Skamarock et al. [31].Note that, if the grid is locally refined by factors of 2 or 3, then these codes suffer from the same problems that nested grids show at the lateral boundaries.Goodfriend et al. [48] demonstrated the long transition required for appropriate turbulent scales to develop at the coarse-fine interface in an LES code using adaptive mesh refinement.
Another approach is to "skip over" the gray zone, as discussed by Zhou and Chow [22].Instead of using the typical 1:3 ratio, using 1:21 for example means that the model transitions from mesoscale to microscale in one jump with a very large domain size such that there is adequate time for turbulent length scales to develop on the finest grid.Over mildly complex terrain, Zhou and Chow [22] performed simulations at 6.4 km resolution down to 25 m in the horizontal, using nesting ratios of 10 or 5 with a total of four nest levels.To accommodate these large changes in resolution, the domain size was increased by nearly doubling the number of grid points in the horizontal directions.From a computational cost perspective, this means larger fine domains, but with the benefit of increased accuracy.Such an approach can be coupled with the cell perturbation method (CPM) discussed above for some computational savings.van Veen [52] used a nesting ratio of 1:15 between D02 and D03 with CPM, as shown previously in Figure 5, going from 2250 m to 150 m horizontal resolution.CPM can be used to greatly reduce the upwind fetch required and reduce the nested domain size when such a large nesting ratio is used to (partially) avoid the gray zone.

Representing Turbulence
Atmospheric motions are turbulent at almost all the scales relevant to weather and climate.This includes turbulence within the planetary boundary layer (PBL) and deep convection, although these processes are typically parameterized using different approaches [56].Below, we first briefly review the traditional turbulence closure schemes.The mass-flux approach used in cumulus convection parameterizations is discussed in the following Section 5.The challenges faced by traditional models in the gray zone are addressed.Then, we discuss recently developed turbulence parameterizations, some of which are intended to represent the subgrid-scale turbulence in both PBL flows and deep convection in terra incognita simulations.

Traditional Schemes and Challenges
Traditional turbulence closure schemes in large-eddy simulations (LES) and mesoscale models are based on the concept of eddy viscosity and diffusivity, although LES and mesoscale closures are designed with different assumptions that are appropriate for each regime.In the gray zone, some of these assumptions are violated, and, more importantly, many atmospheric processes become sensitive to the details of closure schemes associated with those assumptions.
In traditional LES codes, the subfilter-scale (SFS) turbulent stress (τ ij ) is often parameterized as where K m is the eddy viscosity, S ij = (∂u i /∂x j + ∂u j /∂x i )/2 is the strain rate tensor, and the overbar denotes a spatial filter.The SFS potential temperature flux (τ θj ), as an example for scalars, is written as where θ is potential temperature and K h is the eddy diffusivity.K h is often determined by dividing the eddy viscosity K m by a turbulent Prandtl number (Pr t ).The value of Pr t is often determined empirically and it typically ranges from 0.3 to 1.However, a recent study suggests Pr t can have much larger variability and K h should be independent of K m [57].In the turbulent kinetic energy (TKE) 1.5-order closure model [58], for example, K m = c m le 1/2 where c m is a constant, l is a length scale, and e is the SFS TKE, which is computed with a prognostic equation.
A key assumption of closure models in the form of Equations ( 2) and ( 3) is that SFS mixing processes are diffusive/dissipative, i.e., resolved scalar variances and kinetic energy are always transferred towards smaller unresolved scales using this closure.This assumption simplifies the representation of the closure model and often, although not always, gives satisfactory results at very fine resolutions.Backscatter of scalar variance and kinetic energy, however, can occur through countergradient fluxes at all scales, and it becomes particularly significant in the terra incognita.For example, when the grid spacing is at kilometer scales, the SFS thermal fluxes in deep convective clouds are shown to be dominated by countergradient fluxes [59].For the stratocumulus-capped boundary layer, where clouds and turbulence are strongly coupled through radiation, proper representation of the countergradient transport at subfilter scales is critical to the fidelity of the simulations [60].
Another assumption in eddy-viscosity closures is that SFS turbulence is isotropic, i.e., that the eddy diffusivity in LES is the same in the horizontal and vertical directions.In the gray zone, SFS eddies may be significantly influenced by mesoscale variability and therefore are anisotropic.For example, Bryan [61] found that to make their simulation of a hurricane match observational studies, it was necessary to use different eddy diffusivities for the horizontal and vertical directions.This is common practice in mesoscale simulations but not as much for large-eddy simulation.Green and Zhang [62] also demonstrated that the use of an anisotropic turbulence model at gray zone resolution is critical for producing the correct structure of a hurricane and for allowing coarser resolutions to be used without compromising simulation quality.The analysis of Kitamura [63] for a convective boundary layer (CBL) case further showed that, in the gray zone, remarkable anisotropy exists even when the grid aspect ratio is close to unity, thus the effect of anisotropy cannot be depicted by simple functions of grid spacing.
Turbulence schemes in mesoscale models are often referred to as PBL parameterizations.Traditional PBL schemes often have an eddy diffusivity based form such as Equations ( 2) and (3), although only the vertical fluxes are computed and the entrainment at the top of the PBL may be treated differently.Some of them rely on local gradients only like LES-type closures, and therefore face the same challenge in representing countergradient transport.Some other PBL schemes are able to produce countergradient fluxes by including a gradient adjustment term (γ), which was first proposed by Deardorff [64], as follows, Recall that in PBL schemes, the overbar is usually interpreted as a spatial average in the Reynolds-averaged Navier-Stokes (RANS) equations.PBL schemes that allow countergradient fluxes are referred as nonlocal schemes because they represent mixing over a depth larger than the scale of the vertical grid spacing.In contrast, local schemes rely on information at adjacent grid points, and thus consider local gradients only and do not produce countergradient mixing.
PBL schemes all assume horizontal homogeneity of the unresolved turbulence.Thus, the eddy diffusivity and viscosity mainly depend on vertical gradients of relevant variables, and they are responsible for one-dimensional mixing in the vertical only.In the terra incognita, when complex terrain is resolved, the assumption of horizontal homogeneity is not appropriate [23].The horizontal mixing in mesoscale models is usually calculated either using computational horizontal diffusion (e.g., [65]) or with a simple eddy-diffusivity model, with the horizontal eddy diffusivity being set as a constant or determined by a two-dimensional Smagorinsky model [66].Such simple treatment of the horizontal mixing process may be inappropriate in the terra incognita, because, at kilometer-scale resolutions, SGS horizontal mixing can significantly alter the characteristics of many atmospheric processes, such as distribution of pollutants [66], the organization of tropical convection [67], and the intensity of hurricanes [68].
A challenge for modeling in the gray zone is transitioning from use of a PBL scheme to a LES scheme, and thus one solution is use of a scale-aware turbulence closure.This is especially necessary for the simulations using a range of grid resolutions.Simulations using grid nesting or local grid refinement can involve not only the gray zone, but also traditional mesoscale-model and/or LES grid spacings.An ideal turbulence closure in these simulations should be able to reduce to simpler forms for mesoscale models or LES at sufficiently coarse or fine resolutions.A theoretical motivation for developing scale-aware schemes is the assumption that a model that can continuously transform between LES and PBL schemes is likely to work in the terra incognita as well [69].In addition, scale-aware parameterizations are desired because simulation results in the terra incognita may become resolution-dependent for some phenomena, such as the organization of convective cells in the CBL [13,55] and the size and strength of convective storms [70].
Lastly, for simulations involving moist processes, the parameterization of microphysics is often critical in determining the characteristics of clouds [71].Turbulence schemes can interact with microphysics and therefore modulate the strength of storms additionally in this indirect way [72].Thus, unifying the parameterization of turbulence and microphysics is another challenge in the terra incognita.

Recent Developments
Many recent efforts have been devoted to developing adequate turbulence closure schemes to meet the challenges faced in the gray zone, with different strategies and levels of complexity.Most of these new schemes have an origin in existing turbulence closures for the previous or current generation of LES or mesoscale models.Below, we first discuss the development of LES-type closures, and then review the new schemes originating from PBL schemes and other approaches.
Traditional LES turbulence closures have been modified to explicitly add scale-awareness to them.Kurowski and Teixeira [73] introduced a scale-adaptive formulation for the mixing length scale in a TKE closure, in which the mixing length increases with grid spacing until it saturates close to the boundary-layer length scale.Kitamura [74] also manipulated the length scale in a TKE model, while addressing the anisotropy of turbulence, for applying it in the terra incognita.A different approach was adopted by Bhattacharya and Stevens [69], who applied Reynolds averaging to the LES equations using a TKE closure, and obtained a system with two TKE equations, one corresponding to the turbulence at scales smaller than the vertical grid spacing, and the other corresponding to eddies residing in the scales between the vertical grid spacing and the boundary layer depth.
Another pathway of developing turbulence closure models for the terra incognita is to use simplified prognostic equations of SFS fluxes to improve upon traditional eddy-diffusivity closures [6].Such models have been demonstrated to be better than simple eddy-diffusivity models and can produce countergradient fluxes [6,75,76].However, integrating the prognostic equations for all SFS fluxes means a substantial increase in the computational cost [76].An alternative is to neglect the material derivatives of SFS fluxes in the full conservation equations and solve the resulting algebraic equations.This idea of an algebraic model can be traced back to the 1970s [77,78], and has been recently applied to the LES of atmospheric boundary layers by Enriquez [79], Shi et al. [80], and Shi et al. [81].Shi et al. [81] developed an implicit generalized linear algebraic subfilter-scale model (iGLASS) and applied it to stratocumulus-capped boundary layer simulations at traditional LES and at gray-zone resolutions.They found that iGLASS significantly outperformed traditional LES closures at both resolutions and maintained good vertical structure in the boundary layer.iGLASS results were similar to the dynamic reconstruction model (DRM) which also allows backscatter.
The DRM [82] is a more recent LES closure with different design, which has been applied to simulations of clouds in the terra incognita [60,83].It employs explicit filtering and reconstruction, recognizing the different effects of the explicit LES filter and grid discretization.DRM partitions the SFS motions into resolvable subfilter scales (RSFS) and subgrid scales (SGS).Its SFS flux of potential temperature, for example, is modeled as follows, where the tilde denotes the effect of discretization, the overline denotes the explicit filter, and star denotes reconstructed variables.The SGS component still has an eddy-diffusivity based form, but a dynamic eddy diffusivity model [84] is used and an anisotropic version has been developed [83].
One advantage of the dynamic model is that it allows the eddy diffusivity for scalars to be independent of the eddy viscosity.The latter traditionally determines the former through an empirical turbulent Prandtl number, and this simplification is found to cause significant error in representing scalar mixing for some dynamic regimes [57].The RSFS part can produce countergradient fluxes and is anisotropic, which has been found to be important in the terra incognita.The DRM was tested over Askervein Hill and was able to reproduce mean profiles and intermittent recirculation in the lee of the hill, which was not captured by traditional models [85].Shi et al. [60,83] examined the role of DRM in producing backscatter in shallow and deep clouds and found improved simulation results at gray zone resolutions.The mixed model tested by Moeng et al. [86], which is a simplified version of the zero-order DRM, was applied to simulations of deep convection in the terra incognita and also showed good performance [59].Simon et al. [40] found that in convective boundary layer simulations, the transition from LES to mesoscale does not always occur at the same resolution; instead, the transition depends on the closure model used.With the DRM, the gray zone shifted to coarser resolutions, meaning that the DRM could more accurately represent turbulent motions at coarser resolutions than traditional eddy-viscosity closures.Finally, Zhou and Chow [22] showed that, in the stable boundary layer gray zone over real terrain, the model results were highly sensitive to the choice of closure model; here again the DRM was able to better represent turbulence compared to standard models.Traditional PBL schemes have also been extended for simulations in the terra incognita recently.Zhou et al. [66] improved the representation of horizontal mixing for PBL schemes by relating horizontal eddy diffusivity to the characteristic length and velocity scales of the CBL.To make PBL schemes become scale adaptive, Ito et al. [87] modified various length scales in the traditional Mellor-Yamada model [88] by relating them to the horizontal grid spacing normalized by the CBL height, and improved other details, such as the parameterization of dissipation, of the PBL scheme.Shin and Hong [89] introduced explicit grid-spacing dependency functions to a nonlocal PBL scheme.Boutle et al. [90], who adopted a pragmatic approach, blended a 1D PBL scheme and a 3D LES scheme with the ratio of grid spacing to boundary layer depth as a key blending parameter.
A nontraditional PBL scheme that has been extended for the terra incognita is the eddy-diffusivity mass-flux (EDMF) model [91,92].Differing from traditional nonlocal PBL schemes, which employ a gradient-adjustment term, the EDMF scheme uses a mass flux term to represent countergradient transport, e.g., where M is a convective mass flux, and the subscript u refers to the strong updrafts.As discussed in the previous section, the concept of mass flux derives from cumulus convection parameterizations [93] and is supposed to represent coherent updraft-downdraft plumes in convection.The original EDMF formulation is not suitable for the terra incognita because many assumptions about mass flux, such as small fractional area of convective clouds and statistical quasi-equilibrium (QE), become invalid in the terra incognita [30,94].Neggers [95] developed a ED(MF) n model by applying a bin-macrophysics framework, in which a system of discretized size densities is used to represent convective plumes.
The benefit of such a bin-macrophysics approach is that scale-awareness is included in the basis of the closure models.Tan et al. [96] included explicit time dependence in an extended EDMF scheme, in which the prognostic plumes can interact with the environment and cover arbitrarily large area fractions of a grid cell.
The assumed probability density function (PDF) method has also been used to develop high-order closure schemes.Examples include the Cloud Layers Unified by Binormals (CLUBB) scheme [97] and the Simplified High-Order Closure by Bogenschutz and Krueger [98], both of which have been used in simulations at kilometer-scale resolutions [98,99].In these schemes, a double Gaussian functional form is assumed for the joint PDF of vertical velocity, temperature, and moisture content.Particularly for the terra incognita, they are better than the traditional mass-flux approach in that they do not need to assume the fractional area of convective clouds is small.In addition, because SGS variability is embodied in the PDFs, these schemes can drive microphysics schemes in the simulations involving clouds.
Stochastic methods are another sophisticated approach used for developing new turbulence closures for the terra incognita.The transition from a deterministic to a stochastic framework is partly motivated by the invalidation of the quasi-equilibrium assumption used in conventional convection parameterizations [ All of the above new turbulence closure models have shown some success in the terra incognita.Some of the new schemes focus on dry boundary-layer flows only, while others (DRM, EDMF-based models, and CLUBB) attempt to become unified closure schemes which are not only for the turbulence in boundary layer and shallow clouds, but also for the SFS mixing in deep convection.Some preliminary comparison studies of traditional turbulence schemes suggest that 3D closure models are better than 1D PBL schemes for simulating moist convection (e.g., [102][103][104]).Intercomparison studies including the recently developed schemes, especially those nontraditional approaches, will definitely be needed in the future to help illuminate the most suitable framework for turbulence closures in the terra incognita and in particular over complex terrain.

Convection in Mountainous Terrain
Continuing the theme of parameterizing unresolved motions, we move next to considering schemes used to represent convection over mountainous terrain.Kirshbaum et al. [105] gave a recent overview over moist orographic convection and Guichard and Couvreux [106] reviewed cloud-resolving models.Generally, the lifting of air masses in conjunction with mountains enables them to overcome convective inhibition and trigger the generation of new convective cells.The lifting can be of mechanical and/or thermal origin.Differential heating over sloping terrain drives thermal circulations that provide ascent for convective initiation.Moreover, these mountain circulations efficiently transport moisture from the adjacent plains into the mountainous region (alpine pumping [107]).While simulations using a grid spacing of 6-7 km are generally able to reproduce the mountain-plain circulation [108], the horizontal moisture transport may be overestimated even in simulations using a grid spacing of 1 km due to exaggerated thermal circulations.Mountain ranges can lead to the development of quasi-stationary convective bands (e.g., [109][110][111]) as well as propagating rainfall systems that are excited by sensible heating over elevated terrain.

Parameterization of Convection
In the following, an overview of the basics underlying the parameterization of convection is given, and some of the challenges of these schemes in the gray zone within mountainous terrain are assessed.A commonly used approach to parameterize convection is the so-called "bulk mass-flux approach" [112].The mass-flux approach as first pioneered by Arakawa and Schubert [93] approximates the vertical convective flux (of heat, moisture and momentum) w χ of a variable χ as: where χ can be temperature, the total water mixing ratio or horizontal winds, )w e and assuming that there is no mean ascent or descent (w ∼ 0), Equation (7) reduces to: where under the so-called "small area approximation" (i.e., w e ≈ w and χ e ≈ χ if σ u,d 1) the mass flux in the environment has been neglected (M e = 0).Thus, the total vertical flux is approximated by the flux within an updraft plume and a downdraft plume.
The tasks for a convection parameterization scheme thus consist of: 1. predicting the mass fluxes M u and M d 2. predicting the value of χ within the updraft (χ u ) and the downdraft (χ d )

determining where convection should occur
The mass flux is parameterized by an entraining plume, where entrainment dilutes the plume and detrainment represents the mass transfer into the environment: where we relate the fractional entrainment (ε) and detrainment (δ), and the mass entrainment (E) and detrainment (D), respectively.The crucial point lies in determining the mass flux at cloud base, M u,cb for the updraft, as this will scale the entire vertical profile of M u by solving Equation (9).Likewise, a start value for M d is required to solve Equation ( 9) for M d .This constitutes the so-called "closure" of the scheme.χ u and χ d are provided by the bulk cloud model, which describes how χ is modified during the ascent/descent in the updraft/downdraft.Thereby, the plume is assumed to be in steady state.Finally, the trigger function of the scheme decides when and where the convection scheme is activated.This could, e.g., be the thermodynamic state and large-scale ascent of air in the boundary layer, or the moisture convergence (e.g., [113]).Parameterization schemes for convection have a number of well-known problems in reproducing the observed convective precipitation (e.g., [105,114] and references therein).Over the Alps, models with parameterized convection show weak precipitation intensity and a low wet-day frequency compared to observations [115].In complex terrain, parameterization schemes have shortcomings in reproducing the distribution of precipitation with altitude correctly [11], which is especially true for convective precipitation [116].Going to higher resolution an improvement in the precipitation distribution is generally observed, which can partly be attributed to an improved representation of topography [117].Langhans et al. [118], however, noted that the discrepancy in the distribution of precipitation between cloud-resolving simulations and simulations using parameterized convection can mainly be attributed to the behavior of the convection parameterization rather than to the under-resolved small-scale topography.It has also been noted that parameterization schemes strongly overestimate precipitation on the windward side of the mountains, while underestimating precipitation in the lee (windward-lee effect) [119].Simulations at 1 km resolutions in the same study showed an improved spatial distribution of precipitation in comparison to the 7 km simulation using a parameterization scheme.Several aspects to explain this discrepancy were put forward, including the inaccurate representation of the flow at coarse resolution and the missing motion of convective cells and hydrometeor advection.Thereby, also the fall speed of the hydrometeors given by the microphysics scheme has an influence on the precipitation distribution.
Simply refining the grid spacing may not always yield improved simulation results.Convection schemes were originally developed for grid spacings in the order of 10-100 km.Going to finer and finer grid spacings, the convection parameterization scheme should ideally converge to what would be obtained from a cloud-resolving model (CRM).As the grid spacing approaches the size of convective clouds, several of the underlying assumptions of the bulk mass-flux parameterization break down.The requirement that σ u and σ d 1 is no longer valid.Furthermore, within a model column, convection is not uniquely determined by the large-scale state of the atmosphere but a number of possible realizations of subgrid-scale convection exist.In addition, convection excites gravity waves that propagate into neighboring grid cells.The horizontal advection of heat, moisture, or momentum from neighboring grid cells is typically neglected in parameterization schemes for convection.While this may be a fair simplification for classical mesoscale grids, it is not a good assumption for kilometer-scale models.This lack of horizontal communication poses severe problems for the simulation of orogenic propagating precipitation systems (e.g., [120]).
Several approaches have been taken to tackle the challenges associated with grid refinement.Within stochastic parameterizations, for example, convective fluctuations about the ensemble average response of convection to large-scale destabilization can be introduced.Due to the lack of time-scale separation, a memory of subgrid states from previous time steps is required.Examples of stochastic schemes include the work by Plant and Craig [121] or Sakradzija et al. [122].Horizontal communication between adjacent model grid boxes can for example be achieved by introducing cellular automata into a traditional one-dimensional entraining plume model (e.g., [123]).The introduced cellular automata enable convective cells to travel against the mean flow.Many stochastic schemes have been designed and evaluated in an idealized context over flat terrain, so information about their performance in mountainous terrain is still lacking.For a review on stochastic parameterizations see, e.g., the work of Berner et al. [124] and the discussion by Serafin et al. [125].
A more pragmatic approach to modify the convective mass flux within the gray zone is to ramp down the activity of the convection scheme to achieve a convergence towards a CRM as the grid is refined.As the updraft area in a grid box approaches the grid size, the parameterized subgrid convection gradually decreases (e.g., [30,126,127]).In practice, however, the parameterization scheme for convection is either switched off entirely, or only the shallow convection parameterization is retained when approaching the kilometer-scale.As shown in Panosetti et al. [12], for example, a scheme for shallow convection can have detrimental effects on simulated mountain circulations and the associated convection in kilometer-scale simulations, as it overestimates the vertical moisture transport and thereby removes the moisture necessary to build up shallow clouds.
A further challenge in the parameterization of convection is the separate treatment of convection and turbulence by individual schemes.In reality, moist convection and turbulence are tightly coupled.Thus, unified boundary-layer schemes that treat the local effect of turbulent mixing and the non-local effect of convective transport are being employed (e.g., [96]).The non-local convective-transport term of theses schemes fundamentally suffers from the requirement of scale adaptivity and needs to be adapted across the gray zone, as discussed above in Section 4.2.

Explicit Representation of Convection
At kilometer-scale resolutions, deep convective clouds can be explicitly resolved and thus the parameterization of convection can be switched off.Instead of removing the instability, as convection parameterization schemes do, this approach allows convective instability that can lead to the vertical development of clouds, strong downdrafts and gust fronts, and severe weather such as thunderstorms and rain-showers.Heat, moisture and momentum are thus distributed vertically by the convective-scale circulations, updrafts and downdrafts, which are explicitly resolved.In addition to switching off the convection parameterization, there are additional advantages of explicit representation of convection, such as better representation of interactions with complex topography, surface fields and boundary layer processes, which are of key importance for developing convection.
This approach is currently used for numerical weather forecasting at fine resolutions, but with recent model developments and advances in computational power, it has become feasible to use it for climate studies on regional, see e.g., [29,128,129] and continental scales [19,130].To represent long-term climate, the aforementioned studies performed simulations over decade-long time scales.Some studies have even applied km-resolution to global scales [18,131,132], but due to the computational costs they were restricted to simulations of a few weeks.To reduce the computational costs of global convection-resolving simulations, two-or three-dimensional cloud-resolving models can be embedded within each grid column of the global model [133,134].This approach is referred to as superparameterization.
It is still, however, not clear what grid resolution is sufficient to represent resolved scales well.Weisman et al. [135] demonstrated that the use of horizontal grid spacing larger than 4 km overestimates the vertical mass transport produced by convection due to the inability of the coarse-resolution simulations to properly represent non-hydrostatic effects.Thus, it is generally accepted that the grid-spacing of 4 km or smaller is sufficient to explicitly resolve deep convection e.g., [129], even though the statistical properties of it do not yet converge (see Section 5.4) [15,136].
Explicit convection at horizontal resolutions of 4 km can represent large convective systems well and leads to improved model simulations versus coarser resolution models (e.g., [137]), but is not enough to represent small showers.
The largest improvements are found in the representation of the diurnal cycle of precipitation, and the frequency and intensity of heavy summer precipitation [11,19,128,129,138].As an example, Figure 7, taken from Ban et al. [128], shows the comparison of convection-resolving and convectionparameterizing simulations with observations.It can be seen that long-standing problems-such as too early onset of precipitation, underestimation of heavy precipitation, and overestimation of the frequency of light precipitation-are substantially reduced when convection is explicitly resolved.These improvements are more pronounced at hourly time scales, while at daily time scales, the two approaches can give similar results of relatively good performance.However, as can be seen, in the case of the convection-paramterizing model, this is for the wrong reasons.To further demonstrate and explain the difference between different modeling approaches, we look into a few snapshots of summer afternoon precipitation over the greater Alpine region (Figure 8) at 50, 12 and 2 km horizontal resolution in comparison with observations.The sub-daily precipitation observations, available over Switzerland [139], show that summer afternoon precipitation is characterized by convective cells with high intensity, which are realistically reproduced in the convection-resolving model at 2 km horizontal resolution.On the other side, convection-parameterizing models at 50 and 12 km horizontal resolution show widespread precipitation with too low intensity.Although most of the studies on high resolution climate modeling that explicitly resolve convection are focusing on precipitation, further improvements are also found for the simulation of snow-pack [140], clouds [141], and representation of local wind systems [142].It has also been shown that km-resolution models can reproduce severe convective events such as tornadoes, damaging thunderstorm wind gusts, hail (e.g., [143]), meso-scale convective systems [144] and tropical cyclones [145].
In addition to a better representation of convection and their associated spatial scales, explicit treatment of convection can also alter the climate change signal of heavy precipitation events.Over complex mountainous areas such as the Alps, simulations with explicit treatment of convection show a smaller increase in heavy precipitation than in the simulations with parameterized convection [128], while the changes in the seasonal mean precipitation appear robust across the two modeling approaches [29,128,146].

Structural and Bulk Convergence
As indicated above, kilometer-scale resolution is insufficient to fully resolve convective clouds.Typically, a computational mesh with a horizontal grid spacing of ∆x can adequately resolve structures with wavelengths exceeding about 6∆x [147].Thus, a horizontal grid resolution of 2 km will merely resolve wavelengths larger than 12 km.It is evident that such models cannot fully represent the key dynamical aspects of convective clouds (such as updrafts, downdrafts, cold-air pools, entrainment and detrainment).Proper convergence requires substantially higher resolution, as updrafts are often as narrow as a kilometer [45].More specifically, updrafts simulated at km-scale are too laminar.However, the promising results of km-scale models-in both numerical weather prediction and climate applications-suggests that convection is represented to some extent.
For a quantitative assessment of these characteristics one can employ the terminology of "bulk" and "structural" convergence.Here, bulk convergence addresses horizontally-averaged variables related to the water and energy budgets, such as the moist diabatic heating rates, while structural convergence addresses the statistics and scales of individual clouds and updrafts.A recent study has investigated these two types of convergence properties in idealized simulations of diurnal convection over land, using resolutions from 8 km to 500 m [148].It was found that bulk convergence generally begins at ∆x = 4 km, while structural convergence is not yet fully achieved even at the highest resolutions considered (Figure 9).These results are consistent with [41], who indicated the presence of bulk convergence at resolutions around 2 km in real-case simulations over the Alps.The results are also consistent with a number of studies who found that higher resolution is needed for structural convergence [15,149,150].The question of bulk convergence is of paramount importance to climate research.There is the well-founded hope that global simulations at kilometer-resolution could help to better constrain key uncertainties in climate change projections, such as the global climate sensitivity.For such applications, structural convergence is not needed, but an adequate representation of the underlying feedbacks requires some kind of bulk convergence.Current assessments of compute power and model development indicate that global kilometer-scale climate models could become available within the next 10 years [20], while achieving structural convergence at resolutions of 50-100 m would require much longer [151].

Topographic Datasets
In the early days of numerical modeling, topographic datasets were not yet readily available, and certainly not in the digitized versions that are commonly available in a multitude of resolutions today.Rather, the first orographic maps used in atmospheric models were read out of topographic charts.For instance, Berkofsky and Bertoni [152] constructed global topographic datasets from aeronautical charts at a scale of 1:1,000,000.The topographic elevations were averaged "by eye" over one-degree squares.At the time, even the aeronautical charts were incomplete, and values "in unexplored regions were obtained by reasonable interpolations and extrapolations of known values in adjacent regions" [152].The data were then averaged over five-degree squares and tabulated in two tables published in the article (Figure 10a).Even in the 1970s, publication of topographic datasets in the form of tables was still common in the meteorological literature.For instance, La-Valle et al. [153] published a polar stereographic map with a grid spacing of about 90 km for the European-Mediterranean region, on a grid that was intended for an early regional NWP model.
Not surprisingly, the heights of topographic obstacles were lowered dramatically by the procedures followed.For instance, in the two studies mentioned, the peak height of the Alps amounted to 1460 m and 1850 m, respectively, rather than the true peak height of 4808 m.
One of the early (then high-resolution) global datasets was ETOPO5.This is a gridded dataset of worldwide bathymetry and topography derived from several different sources at a resolution of 5 arc-minutes (about 10 km).The dataset was first published on a CD in 1988 by the US National Geophysical Data Center (NGDC).At the time of ETOPO5, better datasets did actually exist in many regions, but it was difficult to get access to such data due to traditional copyright and security restrictions imposed by most providers of topographic data.
A major jump forward came with the GLOBE dataset, the first version of which became publicly available in 1998 [154].It provided global coverage at a spatial resolution of 30 arc-seconds (about 1 km).Later versions of this dataset benefited from the exploitation of remote sensing, such as the joint Space Shuttle Radar Topography Mapper mission in 1999, which provided complete global coverage with interferometric synthetic aperture radar (SAR) data.In 2015, NASA publicly released a refined 1 arc-second dataset (about 30 m resolution) from its Shuttle Radar Topography Mission (SRTM), which is also based on SAR [155].the hand-digitized topography of North America at a resolution of 5 degrees latitude and longitude from [152].The data are tabulated in units of 100 feet.(b) Surface pressure in an early simulation of Alpine Lee Cyclogenesis from [156], with the red line indicating the wall mountain.

Traditional Coordinate Systems and Challenges
The main challenge arising with mountains in atmospheric models is the representation of topography in the discretized numerical grid, in simulations spanning a large range of resolutions.Historically, the bottom boundary of the computational domain has been transformed to align with the terrain, referred to as "terrain-following" coordinates.This transformation has several advantages, as detailed below, but can also introduce numerical errors with steep terrain.In particular, depending upon the horizontal resolution employed, terrain-following coordinates implicitly smooth the topography, thereby significantly reduce its actual height, and thus underestimate the blocking effect of high topography.This problem was recognized already in the 1970s, and the use of "wall mountains" was proposed instead of terrain-following coordinates (for an early example, see Figure 10b).A second challenge is grid quality when spanning multiple scales.The grid aspect ratio (i.e., ∆x/∆z) has been shown to affect simulation accuracy, and yet it is common to refine the horizontal resolution while maintaining the same vertical grid.Use of these traditional methods and challenges that arise in the gray zone are described below.

Pressure-Based Terrain-Following Coordinates
Traditional pressure-based coordinates are intimately tied to the use of the hydrostatic approximation.These coordinates are thus not amenable to km-scale resolutions, but serve to highlight the key challenges and solutions when representing topography in atmospheric models.In fact, two fundamental ideas were introduced with pressure coordinates, and elements of these ideas are relevant throughout all coordinate formulations available today.The first idea is the terrain-following coordinate (discussed below) and the second idea is the step coordinate system to be discussed in Section 6.3.1.
Pressure-based terrain-following coordinates were introduced by [157].In (x, y, p) space, the computational domain has an irregular lower boundary condition, as pressure surfaces intersect the ground, in particular in the presence of topography (Figure 11a).In addition, the lower boundary in p-coordinates is time-dependent, as the surface pressure p s = p s (x, y, t) is changing in time.Terrain-following coordinates may be defined through σ = const, where in the simplest case σ is This formulation yields a monotonous function from σ = 0 at the model top (p = 0) to σ = 1 at the surface (see Figure 11b).With this coordinate transformation, the vertical velocity changes from The resulting implications on the governing set of equations are surprisingly mild, although the mathematical details of the coordinate transformation are complex (e.g., [158]).For instance, the horizontal momentum equation in the x-direction changes from its pressure-based version Here, D/Dt is the total derivative (expressed in the respective coordinate system), f v the Coriolis term, and φ ≈ gz the geopotential.Note that the horizontal derivatives in the pressure-gradient term are taken on p and σ-surfaces, respectively (as indicated by the subscripts).
This approach offers several advantages.First, it maps the atmospheric domain under consideration to a rectangular computational mesh whose data structure is well suited for implementation in numerical models.Second, the terrain-following coordinate transformation yields a simplification of the lower boundary condition (because the transformed vertical velocity in computational space vanishes on the bottom boundary where σ = 1).Third, this alignment of the topography with σ = 1 also simplifies the coupling of the dynamical core of atmospheric prediction models with boundary and surface-layer parameterization schemes.Despite these advantages, the σ-coordinate transformation also introduces some challenges.In particular, while the pressure-gradient term is represented as one single term in the p-based system in Equation (11), two terms are required in the σ-based system in Equation (12).In cases without topography, the second term is very small.However, in regions of complex topography, both the terms are much larger than the pressure gradient itself, and this leads to potentially serious truncation errors-a problem referred to as the "pressure gradient problem".In the horizontal direction pressure varies very slowly with the synoptic-scale pressure gradient (i.e., ∂p/∂x ≈ 10 hPa/1000 km = 10 −2 hPa/km) while in the vertical it varies very quickly with the hydrostatic equation (∂p/∂z = −gp/(RT) ≈ 100 hpa/km).As a result, there are large cancellations between the two terms, and (relatively) small truncation errors in each of the two terms in Equation ( 12) may lead to serious inaccuracies in numerical estimates of the pressure gradient term.To reduce this kind of truncation error, most models introduce a hydrostatically balanced reference state, and the prognostic equations are formulated for the perturbation quantities.However, this approach does not fully eliminate the problem, as deviations from the reference state may be significant, for instance in the vicinity of strong inversions or horizontal temperature gradients.Several approaches have been proposed to further reduce these difficulties (e.g., [159,160]) but a fully satisfactory solution does not exist.Inaccuracies can also arise in the calculation of horizontal gradients for other terms, such as (implicit or explicit) horizontal diffusion (e.g., [161]).

p = const
There are several important generalizations of the σ-coordinate systems.First, most atmospheric models use an uneven spacing of the model levels, usually to enhance the representation of processes in the lower boundary.This approach also enables a cost-effective representation of the upper atmosphere in models with a high model top.Second, many models use hybrid coordinate systems.The most common system employs a transition from σ to p-coordinates [162,163].This improves the representation of stratospheric and mesospheric processes, where the flow is usually quasi-horizontal and where residual topographic patterns in the coordinate surface are undesired.A further generalization has been proposed by Zhu et al. [164].They developed a formulation with smooth transitions from σ to θ to p-coordinates, where θ denotes the potential temperature.Such a system has attractive dynamical properties in the stratosphere.

Height-Based Terrain-Following Coordinates
The basic form of height-based terrain-following coordinates was introduced by Gal-Chen and Somerville [165].The idea closely follows that of the σ-coordinates, but with a z-based formulation.The coordinate surfaces are defined by ẑ = const, with ẑ(x, y, z) Here, H denotes the height of the model top, h(x, y) the topographic height, and ẑ is the new vertical coordinate ranging from 0 at the surface to H at the model top.
For the transformation, the governing equations are rewritten in a vector-invariant form, which always retains their shape irrespective of the underlying coordinates.For instance, the continuity equation in z-coordinates reads ∂ρ ∂t + ∇ • (uρ) = 0 ( In vector-invariant notation this equation remains formally valid in any curvilinear coordinate framework.More specifically, with an alternate vertical coordinate ẑ(x, y, z), the transformed continuity equation can be expressed as where û = (u, v, ŵ) denotes the three-dimensional wind in transformed coordinates, with the transformed vertical velocity defined in the usual way as the total derivative ŵ = D ẑ Dt (16) and where the transformed density is ρ = J −1 ρ ( Here, J = ∂ ẑ ∂z is the Jacobian (or more precisely, the determinant of the transformation's Jacobian).It can be seen that Equations ( 14) and ( 15) are formally identical, but the meaning of the different terms is changing.For simple equations (such as the continuity equation), the transformation is rather straightforward, but for more complicated equations (such as the viscous terms of the Navier-Stokes equation) the transformation can become very cumbersome.
Figure 12a shows the distribution of vertical coordinates in a vertical cross section using the above-discussed approach.One serious disadvantage is apparent: while at low levels the coordinates have to follow the rough underlying topography by design, there are topographic imprints on the coordinate levels even in the upper troposphere and lower stratosphere.This problem can be alleviated partially by hybrid coordinates, but not completely.As the coordinate formulation (Equation ( 13)) relies on a linear decay of terrain features with height, the spectrum of vertical variations in coordinate height will contain small-scale variations, even at elevated levels.It is evident that the associated grid deformations are detrimental to the quality of the numerical solution.Consider for instance a meteorological situation where the upper-level flow is quasi-horizontal and not experiencing the underlying topography (as is actually often the case).In such circumstances, it is inconvenient from a numerical perspective to use a conventional terrain-following coordinate, as a quasi-horizontal flow will-in computational space-ascend and descend through the computational mesh with the underlying topographic variations.
The difficulty under consideration can be quantified using the order-of-accuracy measure of numerical analysis.This measure is used to assess one of the most important properties of numerical algorithms-namely, how the truncation error converges to zero as the computational mesh is refined [166].Theoretical analysis shows that the grid deformations may actually dominate the truncation error, in particular when large-scale smooth flow properties are advected across small-scale grid deformations [167].It is thus of advantage to use terrain-following surfaces that do not exhibit a linear decay of terrain features with height (as is the case with standard terrain-following coordinates), but rather a faster (e.g., exponential) decay with height (see Figure 12b).Several versions of this idea have been proposed [163,[167][168][169] and implemented in a number of limited-area and global models [167,[170][171][172] .
Terrain-following height coordinates are currently a preferred method in mesoscale research models, as well as numerical weather prediction and climate models.The main advantages have already been discussed in the context of σ-coordinates: a rectangular mesh in computational space, simplified implementation of lower boundary conditions, and an efficient coupling to the planetary boundary layer and soil schemes (see Figure 12).Several of these advantages are particularly appealing in lower-resolution models, where terrain-following coordinates have small slopes and are quasi horizontal.With the advent of high-resolution models, the balance of advantages and disadvantages may shift.

Alternative Gridding Techniques
As atmospheric models are used at increasingly fine resolutions, very fine-scale details of the terrain can be captured.For mountainous terrain, this means that the terrain slopes represented by the model become steeper.At very high resolutions, there is interest in capturing the effects of vertical or near-vertical features, such as those formed by cliffs, escarpments, and even the built environment.As mentioned in the previous section on terrain-following coordinates, models using those coordinate systems are challenged by numerical errors which increase in the presence of steep slopes.In an effort to develop atmospheric models that can handle arbitrarily complex terrain, alternatives to terrain-following coordinates have been developed (for more detail, see [10]).
Alternative methods commonly applied when transitioning from simulating the mesoscale to microscale include wall mountains, step coordinates, masking methods, and immersed or embedded boundary methods.A section of a terrain-following grid is shown in Figure 13, along with alternative methods.This figure focuses on details of a peak of Granite Mountain in Utah, USA, whose full topography is shown in Figure 14.

Wall-Mountains and Step Coordinates
The second key innovation introduced with pressure coordinates mentioned above relates to the introduction of wall-mountains.The basic idea goes back to Egger [173], who argued that σ-coordinates would do well with representing large-scale lifting and ascent associated with major mountain ranges such as Greenland and Antarctica, but not with orographic blocking associated with steep mid-latitude mountain ranges such as the Alps.Indeed, a high knife-edge mountain would be smoothed to zero in σ-coordinates, when in fact it would stop all air from flowing normal to it.Egger [173] designed a wall-mountain topography, by enforcing the condition v n = 0, where v n is the component of the horizontal wind normal to the mountain ridge.Distinct from σ-coordinates, which approximately conserve the mountain volume (or even increase it if some form of envelope topography is used), wall mountains have zero volume but represent the blocking effect.The approach was successful and led to the first numerical simulation of Alpine lee cyclogenesis [156].A sample of the results is shown in Figure 10b.Of course, this approach is not an option for a weather forecasting model.In particular, in real-case settings the flow response of a mountain range may involve both blocked and lifting regimes depending upon environmental conditions, while the wall-coordinate is only suitable for the former.
[174] designed a "step" mountain coordinate system (also referred to as the "eta" coordinate).In this system, the mountains are constructed from cuboid-shaped blocks with vertical walls, thereby allowing application of the v n = 0 condition at the lateral edges of the building blocks.This was motivated by the difficulties of the σ system with the pressure gradient term, associated with the non-cancellation of errors in the two terms in Equation (12).With step-mountain coordinates, these difficulties do not occur, as the coordinate surfaces are quasi-horizontal.The step-mountain coordinate has served in the Eta operational model at the National Center for Environmental Prediction, but it has difficulties in representing gravity waves and produces errors along the terrain surface where artificial vorticity is introduced at step corners [175,176].
In simulations of flow over urban areas, step-coordinates can be referred to as masking methods.In these simulations, building topography is fit to the grid and boundary conditions are applied within the domain.This method is easy to implement and frequently gives good results, but introduces errors both in the approximation of the terrain and with artificial flow features appearing in step corners.

Immersed and Embedded Boundary Methods
As with step coordinates and masking methods, a non-conforming grid is used with immersed or embedded boundary methods, meaning that the grid is not transformed to conform to the terrain (Figure 13).Instead, the terrain passes through the fluid domain, and the effect of the terrain is represented on the flow through special treatment at the boundaries.Immersed boundaries use a "fictitious" force added to the governing equations to represent the effects of the boundary, while with embedded boundary methods the numerical discretization is modified at cells intersected by the boundary.Both methods yield a more accurate representation of the terrain than use of step coordinates and masking methods, where errors are introduced by the zeroth-order block terrain representation, as shown for example, in simulations of the Salt Lake valley using a horizontal resolution of 850 m and step coordinates [177].
Use of immersed and embedded boundaries first appeared in the engineering literature, where the methods were used to represent objects with complex geometries.Peskin [178] first used the immersed boundary method (IBM) to model two-dimensional flow through the valve of a heart, and the methods have subsequently been applied to a large variety of engineering applications [179,180].Non-conforming gridding techniques have been employed in atmospheric applications as well, especially as higher-resolution simulations allow explicit representation of complex topographies with steep slopes.Early application of immersed or embedded boundaries to environmental applications include simulations of flow over smooth hills in ocean or atmospheric models [181][182][183].Senocak et al. [184] used an immersed boundary to model neutral atmospheric boundary layer flow over flat terrain, focusing on correctly modeling surface stress.Simulations over the urban terrain of Baltimore, MD were carried out in Tseng et al. [185].Lundquist et al. [9,186] implemented the immersed boundary method in the WRF model and used the method to simulate flows over urban and mountainous terrain.Flow over Bolund hill, the site of a field and model intercomparison study [187,188], has been used to validate implementations of a IBM in several different numerical models [189][190][191][192].
Although immersed and embedded boundary methods have been used successfully, several challenges remain for using non-conforming grid methods for atmospheric applications.Engineering applications frequently use high-resolution grids near objects with complex topology, employing advanced gridding techniques to achieve such high resolutions (e.g., unstructured grids, block grid refinement, and adaptive mesh refinement).This high resolution of the boundary layer allows introduction of interpolated boundary conditions with negligible error.Use of non-conforming methods can be more challenging in atmospheric applications than in engineering applications, partially because atmospheric models tend to lack flexible gridding features, but also because the high-Reynolds number flows in atmospheric applications have coarse resolution near boundaries, comparatively.This coarser resolution can introduce significant interpolation errors when alternative boundary condition treatments are used, and therefore these alternative gridding techniques are more accurate for highly resolved flows than for simulations within the gray zone.In comparison to engineering applications, additional unique challenges also arise with the parameterization of atmospheric physics.For example, parameterizations of surface fluxes are based on Monin-Obukhov similarity theory and commonly rely on use of a reference height.With terrain-following coordinates the first grid point is located at a relatively uniform height above the terrain surface, which works well with similarity theory.In contrast, non-conforming grid methods lack uniform grid spacing above the terrain, thus accurate implementation of models based on similarity theory can be difficult.
Recently, these challenges have begun to be addressed in the literature.Early implementations to combine a stress-type boundary condition with the immersed boundary method showed success in simulations of a neutral atmospheric boundary layer; however, due to use of flat topography, the stress condition was simplified to a one-dimensional implementation [184,193].Later methods for combining IBM with a stress boundary condition extended the method to three dimensions, accommodating complex terrain.Several different implementations of IBM with a shear-stress boundary condition have been developed and validated using observations from the Askervein or Bolund hill field experiments [189][190][191][192]194].These simulations, however, were conducted at the microscale using horizontal resolutions ranging from 1 to 5 m, rather than within the gray zone.DeLeon et al. [192] simulated flow over Askervein hill at slightly coarser horizontal resolutions, ranging from approximately 12 to 25 m.Comparisons with observations were noticeably degraded on the coarse mesh.Bao et al. [195] implemented multiple IBM methods into the WRF model and found degraded performance across implementations for flow in an idealized valley using a horizontal resolution of 90 m, indicating that development of an IBM method for the gray zone is still an open area of research.
Additional challenges beyond accurate specification of surface fluxes exist when using non-conforming grids for atmospheric simulations.One is that atmospheric physics models (e.g., radiation, land-surface, microphysics, etc.) require modification to work with alternative grids, because the bottom grid level no longer aligns with the topography.Therefore, each physics model requires modification to correctly ingest values located at the immersed boundary rather than at the lowest grid level and also to neglect portions of the computational domain which are below the terrain surface.It is also important to keep in mind that, even when atmospheric physics models are correctly coupled to the non-conforming grid, they may not work well over complex terrain or within the gray zone.Lundquist et al. [9] coupled an IBM to a subset of atmospheric physics parameterizations (e.g., radiation and land-surface) within the WRF model, providing important surface fluxes of heat and moisture for simulations of valley flows.Arthur et al. [196] extended this work by coupling a topographic shading routine to the IBM and simulated katabatic flow over Granite Mountain during the MATERHORN field campaign using horizontal resolutions of 25 and 50 m (shown in Figure 14).While surface fluxes of heat and moisture are parameterized in these works, a no-slip boundary condition for velocity was used due to the previously mentioned difficulties with accurate implementation of a surface momentum flux at immersed boundaries for simulations within the gray zone.
Further complicating the use of alternative gridding methods is that a boundary condition is required for each prognostic equation and sometimes for other model quantities, such as eddy viscosity or shear stress.For atmospheric models, this can mean setting immersed/embedded boundary conditions for dozens of variables, often with varying types of boundary conditions, leading to a more complex implementation than found in engineering applications.Additionally, although these methods are suitable for representing resolved terrain-features, smaller scale terrain frequently exists and may require special treatment.Tseng et al. [185] demonstrated that 6-8 grid points are required across a building for the immersed boundary method to produce acceptable results.Under-resolved features could be represented as virtual terrain or buildings with areas of increased drag [197].
Finally, challenges remain with obtaining realistic lateral boundary conditions for models using alternative grids, and often these simulations use idealized flow conditions at lateral boundaries instead.Wiersema et al. [35] developed a method for nesting domains using immersed boundary coordinates into those using terrain-following coordinates, and performed a multi-scale simulation of flow over urban terrain (∆x = 2 m) nested within a mesoscale simulation (∆x = 6 km) using a five domain configuration.The urban scale simulation ran concurrently with the mesoscale simulation, receiving time-varying lateral boundary conditions, and validation against the Joint Urban 2003 field campaign in Oklahoma City was performed.This simulation used the WRF model with its framework for concurrent nesting, along with vertical mesh refinement [53] and an immersed boundary method using a surface stress parameterization [194].
Currently, it is clear that terrain-following coordinates produce large errors when the grid becomes skewed over steep complex terrain.This is a significant challenge for producing accurate simulations within the gray zone.Several methods for handling complex terrain have been developed ranging from more accurate discretization methods, to alternative terrain-following coordinates that flatten quickly with altitude, to the alternative coordinate systems described here.Alternative coordinates, specifically immersed and embedded boundaries, have shown great promise in improving the accuracy of simulations in the gray zone over complex terrain, especially for near surface quantities such as wind speed and direction.While these methods show promise, they have primarily been employed in simulations using a dynamical core alone at microscale resolutions.Development of alternative coordinate systems for both fully featured atmospheric models and for use in the gray zone is the subject of ongoing research.

Discussion and Future Needs
Weather predictions at regional and fine scales increasingly play a critical role in everyday life, affecting economic, personal, and environmental decisions, among others.The atmospheric simulations used to make these predictions traverse multiple gray zones in the transition from the mesoscale to the microscale, where traditional modeling approaches may fail.In complex terrain applications, the gray zones include the representation of turbulence and convective processes and also the representation of topography, among others.Because these multiscale forecasts are the drivers not only for weather prediction but also for air quality models, regional climate simulations and wind energy applications, the impacts of improved modeling techniques in the gray zone are potentially huge.
Current mesoscale simulations have already surpassed the typical resolutions they were designed for, crossing the threshold beyond which traditional parameterization schemes may be suitable.At the same time, fine scale models are now becoming feasible over much larger regions, where mesoscale variability must naturally be incorporated, especially over complex terrain.It is clear that atmospheric simulation is a multi-scale endeavor, and one which will see tremendous opportunities in the number and type of applications as modeling techniques conquer these multiple gray zones.It is evident that for the next 10-20 years most operational and research applications of atmospheric models-in the areas of numerical weather prediction, air quality and climate modeling-will refine their spatial resolutions, but most will still not fully resolve all of the key underlying processes.In other words, they will be operating somewhere in the "gray continuum".Experience presented in this article shows that awareness of these challenges is of key importance in advancing the accuracy and predictive power of these models.
The scale-aware formulations discussed in this review paper for grid nesting and turbulence closures from 10 km to 10 m horizontal resolutions across the gray zone or terra incognita will move models toward a consistent representation of turbulent and resolved motions as grid resolutions increase with advances in computing power.For turbulence closures, scale-aware parameterization means incorporating the grid length scale and enabling a smooth transition from PBL to LES-type closures.Moreover, the active role of unresolved turbulence in the gray zone requires new capabilities, such as representing anisotropy and backscatter, to be built into turbulence closures.The increased complexity of the resolved topography at fine scales will additionally mean that models may benefit from 3D turbulence parameterizations.For convective parameterizations, safely traversing the gray zone may mean ramping down the activity of the convection scheme toward cloud resolving simulations or using stochastic methods.Over complex terrain, the sensitivity of orographic precipitation estimates to grid resolution demonstrate the strong interplay between topography, which is increasingly resolved, and the convective representation, either of which may fall within the gray zone.Finally, underlying all of these issues is the representation of topography at different scales.While terrain-following coordinates serve well at coarse resolutions, at very fine scales, alternative gridding techniques, such as immersed boundary methods, enable simulations of very steep topography.The transition across these scales will no doubt bring innovation in new hybrid coordinate systems that perhaps combine the best parts of both terrain-following and IBM approaches in the gray zone.
Future research efforts in all of these areas must be devoted to continue to ensure that mesoscale models are using cutting edge physics parameterizations and numerical methods.Well-organized intercomparison projects will help the community to expedite the progress and identify the best modeling strategies to provide physically consistent simulations across a wide range of scales over complex terrain.

Figure 1 .
Figure 1.Sketch showing representations of eddies of length scale L ∼ 1 km shown on a vertical grid for: (a) a mesoscale model with horizontal grid spacing ∆ >> O(L); (b) a highly-resolved LES with ∆ << O(L); and (c) the gray zone with ∆ ∼ O(L) .

Figure 2 .
Figure 2. Grid nested domains over the Perdigão domain in Portugal at different horizontal resolutions: (a) 6750 m; (b) 2250 m; and (c) 150 m.Box outline in (c) shows location of fourth nest at 30 m resolution.

Figure 3 .
Figure 3. Horizontal distribution of vertical velocity (m/s) from: (a) a 1200 m run; (b) a 25 m run; and (c) the 25 m run filtered onto a 1200 m grid.The dashed line in (a) is the same size as the domain in (b,c) for comparison.Contours are from the horizontal plane at z/z i = 0.1 at 1500 LST.Note the different color bars.Adapted from Zhou et al. [13].

Figure 4 .
Figure 4. Horizontal contours of vertical velocity from: (a) a periodic 400-m resolution run; and (b) the 400-m run nested within a 1200 m periodic domain.Slices taken at z/z i = 0.1 at approximately the onset of convection for the periodic 1200 m run.Adapted from Zhou et al. [13].

Figure 5 .
Figure 5. Vertical velocity (m/s) at 260 m above ground level for a WRF LES grid at 150 m resolution at 1900 LT.Line contours of topography (white-gray) are shown: (a) a large reference domain without CPM is compared to (b) a small domain with CPM and (c) a small domain without CPM.The dashed lines in (a) show the spatial extent of the small domains in (b,c), which are the same as the domain in Figure 2c.The diagonal lines are used for the w-spectra comparison in Figure 6.

Figure 6 .
Figure 6.w-spectra taken over diagonal lines and three domains described in Figure 5, averaged over five time samples, beginning at 1900 LT.
100].In the data-inferred stochastic model developed by Dorrestijn et al. [100], pairs of turbulent flux profiles are pre-computed based on LES statistics, and the transition probabilities of those flux profiles are conditioned on the resolved grid column states.The stochastic parameterization developed by Sakradzija et al. [101] is based on the EDMF framework, but the mass flux closure is applied to an area larger than a single grid cell to ensure the validity of the QE assumption, while the grid-cell-scale cloud fraction perturbation is generated by subsampling the stochastic cloud ensemble.
u and σ d are the updraft and downdraft fraction, respectively, and w is vertical velocity.Primed variables indicate deviations with respect to the domain mean values.Using the upward mass flux M u = ρ • σ u w u , the downward mass flux M d = ρ • σ u w d , and the mass flux in the environment M e = ρ • (1 − σ u − σ d

1 Figure 7 .
Figure 7.Comparison of precipitation from reanalysis-driven simulations with 12 km (CPM12, dashed blue) and 2.2 km (CRM2, solid blue) horizontal resolutions with parameterized and explicit convection, respectively.Cumulative frequency distribution of (a) daily precipitation and (b) daily maximum 1 h precipitation, expressed relative to the number of wet days.Mean diurnal cycle of (c) mean precipitation, (d) wet-hour frequency, and (e) heavy hourly precipitation (p99H, 99th percentile of all hours).From[128] c Copyright 2015 American Meteorological Society.

1 Figure 8 .
Figure 8. Hourly precipitation amount for a deep moist convection event on 12 July 2006 during 15:00-19:00 UTC (from top to the bottom) over the greater Alpine region obtained from a regional climate model at (first column) 50 km and (second column) 12 km with convection parameterization scheme on, (third column) at 2 km with convection parameterization switched off, and (fourth column) gridded hourly precipitation observations [139].The model results present snapshots from decade-long simulations presented by Ban et al. [11].

Figure 9 .
Figure 9.Comparison of (a) structural convergence and (b) bulk convergence for an idealized simulation of diurnal convection over land, for grid spacings from 8 km to 500 m.(a) For the resolutions considered, structural convergence is not yet achieved, as the average cloud area reduces by a factor 4 when the grid spacing is reduced by a factor 2. (b) However, bulk convergence is achieved near ∆x = 4 km, as demonstrated by the grid-point statistics of convective mass flux (from Panosetti et al. [148]).

Figure 10 .
Figure10.Historical papers on topography in atmospheric models.(a) Section of a table showing the hand-digitized topography of North America at a resolution of 5 degrees latitude and longitude from[152].The data are tabulated in units of 100 feet.(b) Surface pressure in an early simulation of Alpine Lee Cyclogenesis from[156], with the red line indicating the wall mountain.

Figure 13 .
Figure 13.Three grids are shown, focusing on the details of a mountain peak.A subset of a terrain-following grid is shown in (a).Two alternatives are also shown: (b) step (eta) coordinates or masking methods; and (c) immersed or embedded grids.

Figure 14 .
Figure 14.Downslope flow on the east side of Granite Mountain during the MATERHORN field campaign during the evening transition (18:55 MST).The WRF model is used with an immersed boundary method coupled to atmospheric physics including radiation, land-surface, and topographic shading routines.Blue isosurfaces indicate downslope velocity and topographic shading is indicated by darkening of the terrain surface.A real terrain image from Google Static Maps API is used.Observational towers are shown as green dots.An animated version of this figure is available in the online supplemental material from Arthur et al. [196].c Copyright 2018 American Meteorological Society.

Author
Contributions: F.K.C. and C.S. designed the paper outline and contributed text for all sections.All authors provided concepts and text for individual paper sections: N.B., Section 5; K.A.L., Section 6; L.S., Section 5; and X.S., Section 4. Funding: F.K.C. and X.S. are grateful for support from NSF grants ATM-1565483 and 1503860 (Physical and Dynamic Meteorology Program).X.S. also thanks the support of HKUST startup funds.K.A.L. received support from the U.S. DOE Office of Energy Efficiency and Renewable Energy (EERE) Wind Energy Technologies Office.A portion of this work was prepared by LLNL under Contract DE-AC52-07NA27344.