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.

A rotating radiative‐convective equilibrium on a sphere is reached using a global atmospheric model with prescribed globally uniform sea surface temperature and no insolation. In such an equilibrium state, multiple tropical cyclone‐like vortices coexist in the extratropics, moving slowly poleward and westward. Many vortices have a lifetime longer than 2 months and travel from the tropics to the polar regions. The typical spacing of simulated tropical cyclone‐like vortices is comparable to the deformation radius, while the production of available potential energy is at a scale slightly smaller than those vortices. It is hypothesized that the growth of tropical cyclone‐like vortices is driven by the self‐aggregation of convection, while baroclinic instability destabilizes any vortices that grow significantly larger than the deformation radius. A weak Hadley circulation dominates in the deep tropics, and an eastward propagating wave number 1 Kelvin‐like mode having a period of 30–40 days develops at the equator. The weak Hadley circulation is found to emerge from an initially quiescent atmosphere due to poleward momentum transport by the vortices.

The sensitivity of stratiform midlatitude orographic precipitation to global mean temperature is investigated through numerical simulations. As a step toward understanding the relative influence of thermodynamic and dynamical processes on orographic precipitation, simple idealizations of Earth’s major north–south mountain chains are considered. The individual terrain elements occupy four islands equally spaced around the Northern Hemisphere of a planet otherwise covered by ocean. Although these mountains have very little influence on the sensitivity of the zonally averaged precipitation to changes in global mean surface temperature, the precipitation on the windward slopes of the ridges is highly sensitive to such changes. When the ridges run between 40° and 60°N, the windward-slope hydrological sensitivity exceeds the Clausius–Clapeyron scaling of about 7% K−1 over the northern half of the barrier, leading to substantial precipitation changes. The annual accumulated orographic precipitation is modified by changes in both the mean precipitation intensity and the changes in the number of hours during which precipitation occurs. The changes in the number of hours with significant precipitation largely results from modifications in synoptic-scale storminess associated with changes in the midlatitude storm tracks. A simple diagnostic model reveals the primary factors modifying the mean orographic precipitation intensity are variations in 1) the moist adiabatic lapse rate of saturation specific humidity, 2) the wind speed perpendicular to the mountain, and 3) the vertical displacement of saturated air parcels above the windward slope. The strong dependence of 2 and 3 on latitude further confirms that changes in midlatitude storminess are a major factor determining the response of orographic precipitation to global warming.

Tropical cyclones (TCs) are expected to produce more intense precipitation under global warming. However, substantial uncertainties exist in the projection of coarse-resolution global climate models. Here, we use deep learning to aid targeted cloud-resolving simulations of extreme TCs. Contrary to the Clausius-Clapeyron (CC) scaling, which indicates a 7% moisture increase per K warming, our simulations reveal more complex responses of TC rainfall. TCs will not intensify via strengthened updrafts but through the expansion of deep convective cores with suppression of shallow cumulus and congestus. Consequently, while localized hourly rainfall may adhere to the CC scaling, precipitation accumulation over city-sized areas could surge by 18%K^-1. This super-CC intensification due to changing TC structure has profound implications for floods and landslides. Estimations using Hong Kong’s slope data confirm this concern and suggest an up to 215% increase in landslide risks with 4-K warming, highlighting amplified threats from compound disasters under climate change.

Kilometer‐scale grid spacing is increasingly being used in regional numerical weather prediction and climate simulation. This resolution range is in the terra incognita, where energetic eddies are partially resolved and turbulence parameterization is a challenge. The Smagorinsky and turbulence kinetic energy 1.5‐order models are commonly used at this resolution range, but, as traditional eddy‐diffusivity models, they can only represent forward‐scattering turbulence (downgradient fluxes), whereas the dynamic reconstruction model (DRM), based on explicit filtering, permits countergradient fluxes. Here we perform large‐eddy simulation of deep convection with 100‐m horizontal grid spacing and use these results to evaluate the performance of turbulence schemes at 1‐km horizontal resolution. The Smagorinsky and turbulence kinetic energy 1.5 schemes produce large‐amplitude errors at 1‐km resolution, due to excessively large eddy diffusivities attributable to the formulation of the squared moist Brunt‐Väisälä frequency ( ). With this formulation in cloudy regions, eddy diffusivity can be excessively increased in “unstable” regions, which produce downward (downgradient) heat flux in a conditionally unstable environment leading to destabilization and further amplification of eddy diffusivities. A more appropriate criterion based on saturation mixing ratio helps eliminate this problem. However, shallow clouds cannot be simulated well in any case at 1‐km resolution with the traditional models, whereas DRM allows for countergradient heat flux for both shallow and deep convection and predicts the distribution of clouds and fluxes satisfactorily. This is because DRM employs an eddy diffusivity model that is dynamically adjusted and a reconstruction approach that allows countergradient fluxes.