## Abstract

Realistic multiscale simulations involve coupling of mesoscale and large-eddy simulation (LES) models, thus requiring efficient generation of turbulence in nested LES domains. Herein, we extend our previous work on the cell perturbation (CP) method to nonneutral atmospheric boundary layers (ABLs). A modified Richardson number scaling is proposed to determine the amplitude of the potential temperature perturbations in stable ABLs, with −1.0 overall providing optimum turbulence transition to a fully developed state (fetch reduced by a factor of 4–5, compared to the unperturbed cases). In the absence of perturbations, turbulence onset is triggered by a Kelvin–Helmholtz instability, typically occurring in the vicinity of the low-level jet maximum. It is found that a turbulent length scale can be used to more accurately estimate the optimum , where *q* is the turbulence kinetic energy, and *N* is the Brunt–Väisälä frequency. In convective ABLs, a perturbation amplitude based on mixed layer temperature variance scaling is proposed: . For that criterion to be optimum, the ratio , where is the wind speed at the top of the capping inversion, and is the convective velocity scale, needs to be incorporated: . This allows us to account for the competing roles of the surface thermal instability and the mean flow advection. For 10, the development fetch is reduced by a factor of 6, while when 3, the use of the CP method does not provide a significant advantage in the ability to generate turbulence, provided a smooth mesoscale inflow.

## 1. Introduction

Realistic and computationally affordable multiscale simulations of atmospheric flows involve coupling of mesoscale and large-eddy simulation (LES) models. Such coupling requires efficient generation of three-dimensional turbulence at the inflow boundaries of the nested LES domains, which are forced with “smooth” mesoscale model fields. Mesoscale inflow does not contain any resolved eddies at the scales of atmospheric boundary layer (ABL) turbulence that the LES model can explicitly resolve with typical grid spacings of 10–100 m. Different approaches to generate turbulence in atmospheric LESs forced by a mesoscale inflow have been proposed, including the use of precursor horizontally periodic simulations (e.g., Golaz et al. 2009; Sauer et al. 2016; Munters et al. 2016), recycling techniques (e.g., Mayor et al. 2002; Nakayama et al. 2012), synthetic inflow generators (e.g., Xie and Castro 2008), and the recently proposed perturbation methods (Mirocha et al. 2014; Muñoz-Esparza et al. 2014). From the latter approach, the cell perturbation method has shown promise in efficiently generating realistic turbulence during a diurnal cycle in full-physics atmospheric simulations (Muñoz-Esparza et al. 2017), as well as other semi-idealized scenarios, including ocean–island interactions and cloud formation (Jähn et al. 2016) and sea breeze over an urban coast (Jiang et al. 2017).

The cell perturbation (CP) method applies stochastic perturbations to the potential temperature field that, through local buoyancy effects, break the two-dimensionality of the incoming flow and accelerate the development of three-dimensional turbulence. Perturbations are arranged as three horizontal cells of 8 × 8 grid points in the region near the inflow boundaries of the LES domain. This is done in order to effectively affect fully resolved modes in the Kolmogorov’s inertial range without being rapidly damped by the dissipative nature of the discretized advection term at higher wavenumbers. The CP method was generalized for neutrally stratified ABLs by Muñoz-Esparza et al. (2015), with all the required parameters related to available mesoscale information. The generalization for neutral ABL conditions uses the proposed perturbation Eckert number = 0.2 to derive the optimum range of potential temperature stochastic perturbations for each case, where is the geostrophic wind (or large-scale pressure gradient), and is the specific heat capacity at constant pressure.

Inflow generation in neutrally stratified conditions has historically received primary attention (e.g., Lund et al. 1998; di Mare et al. 2006; Xie and Castro 2008; Golaz et al. 2009; Nakayama et al. 2012). In the case of the CP method, the optimum = 0.2 results in an efficient formation of hairpin-like vortices that rapidly transition to fully developed streaks in the neutral ABL. In contrast, the effects of atmospheric stability on the transition mechanism and turbulence development have not been explored. Despite that, the perturbation Eckert number has been satisfactorily applied for determining the amplitude of the potential temperature perturbations in both convective and stable regimes with the CP method (Jähn et al. 2016; Muñoz-Esparza et al. 2017; Jiang et al. 2017). Nevertheless, modifications of existing methods through introduction of scalings appropriate for stable and convective ABLs could potentially result in more efficient generation of fully developed turbulence.

Herein, we first extend and then evaluate CP method to nonneutral ABLs. We focus on understanding the mechanisms that are responsible for turbulence transition and development in coupled mesoscale–microscale simulations under different stability conditions. The remainder of the paper is structured as follows. The numerical model, simulation setup, and cases of study are described in section 2. A modified Richardson scaling is proposed and tested in five stable ABLs in section 3, including analysis of turbulence statistics (variance, turbulent fluxes, and higher-order moments) and energy spectra and cospectra. In addition, an explanation for the transition in stably stratified conditions is proposed, which is then used to further refine the proposed scaling. A similar analysis for four characteristic convective ABLs is presented in section 4, proposing an alternative scaling using a thermal variance relationship. This relationship is further optimized based on the nature of the identified transition mechanisms. Finally, section 5 is devoted to summarizing the main findings of this work.

## 2. Methodology

### a. Large-eddy simulation model and simulation setup

We perform simulations of idealized ABLs using the Weather Research and Forecasting (WRF) Model v3.8.1 (Skamarock and Klemp 2008; Skamarock et al. 2008). WRF is a compressible, nonhydrostatic code that solves the Navier–Stokes and energy equations for high Reynolds numbers (no viscous term) with filtering of acoustic modes. Equations are discretized over a staggered Arakawa C grid, with a terrain-following hydrostatic-pressure vertical coordinate. Time integration utilizes a third-order Runge–Kutta scheme, while the advection term is discretized using the hybrid even–odd-order numerical scheme from Kosović et al. (2016). This hybrid scheme improves effective spectral resolution (4), compared to the standard fifth-order upwind schemes typically used by WRF (7, Skamarock 2004). This is particularly important for turbulence-resolving simulations, especially for stable ABLs, where relevant phenomena occur at small scales.

To provide homogeneous inflow conditions to the LES domain, we utilize two domains with flat terrain that communicate via one-way nesting. The parent domain uses a one-dimensional PBL parameterization (Nakanishi and Niino 2006), and lateral boundaries are set to periodic in both directions. In the presence of the horizontally homogeneous forcing prescribed as a vertical sounding profile, this setup results in a uniform flow solution in the *xy* plane across the domain. The nested LES domain receives Dirichlet boundary conditions for all the prognostic equations provided by the mesoscale parent domain. A similar configuration has been used in previous studies targeting the development of turbulence generation methods with the WRF Model (e.g., Mirocha et al. 2014; Muñoz-Esparza et al. 2014).

Implicit filtering of the governing equations is used, which assumes that the attenuation of high wavenumbers present in the differencing schemes acts as a low-pass filter, and thus, the grid size acts implicitly as a filter. For the parameterization of subgrid-scale (SGS) effects on the momentum equations in the LES domain, we use the nonlinear backscatter and anisotropy (NBA) model from Kosović (1997), which is implemented in the open release of the WRF Model (Mirocha et al. 2010). From the two variants of the NBA scheme, we use the one that employs an additional prognostic equation for subgrid-scale turbulence kinetic energy (TKE) (Lilly 1966, 1967), following Deardorff (1980) and Moeng (1984). SGS thermal effects are parameterized using a countergradient approach, and the subgrid length scale for both momentum and heat is corrected to take into account possible length scale reductions due to local stable stratification (Deardorff 1980). Monin–Obukhov similarity theory (MOST) is applied at the first vertical grid point to provide surface boundary conditions for horizontal velocity components (vertical velocity is set to zero) and temperature (Jiménez et al. 2012). No cloud, microphysics, radiation, or land surface processes are considered.

### b. Description of the simulated cases

The number of grid points used on the LES domain for both sets of simulations is 1800 × 900 × 150, in the streamwise, spanwise, and vertical directions, respectively. For the stable ABL cases, the LES domain extent is , , = 13 500, 6750, 750 m, uniformly distributed in the three spatial directions ( = = 7.5 m, = 5 m). At the surface, a cooling rate boundary condition is imposed. This is a common choice when modeling stable ABLs with large-eddy simulation models (e.g., Kosović and Curry 2000; Basu and Porté-Agel 2006; Huang and Bou-Zeid 2013) and allows for a better representation of stronger than weakly stable stratifications (e.g., Basu et al. 2008; Gibbs et al. 2015). Five different stable boundary layers (SBLs) are simulated by varying the geostrophic wind = 5, 10, 15 m s^{−1} (uniform and time invariant) and the surface cooling rate *λ* = 0.25, 0.50, 0.75 K h^{−1}. This setup and range of cooling rates is similar to other idealized LES studies of the stable ABL (e.g., Beare et al. 2006; Zhou and Chow 2011; Allaerts and Meyers 2018). The mesoscale parent domain is run for 10 h of physical time to allow the model to stabilize and reach quasi-steady conditions, prior to initializing the nested LES domain. After that, both domains run simultaneously for an additional 2 h. Table 1 collects the main parameters describing the SBL conditions of the simulations, including resulting boundary layer height , friction velocity , surface heat flux , and Obukhov length . Boundary layer height is estimated as the height at which the wind speed profile is 1.05 before reaching geostrophic conditions. In all the cases, an initial profile with *θ* = 300 K is prescribed from the surface up to *z* = 100 m, followed by a temperature inversion of 0.01 K m^{−1} above.

The resulting profiles of potential temperature and wind speed after 10 h of physical time integration of the mesoscale domain are shown in Fig. 1. As the cooling rate is increased from 0.25 to 0.75 K h^{−1}, the boundary layer height progressively decreases, leading to wider low-level wind maxima with slightly enhanced velocities at the center of the jet and larger potential temperature gradients. As a result, a variety of local wind shear and temperature gradient scenarios is produced, which serves as a reasonable database to evaluate the CP method. To judge the adequacy of the grid size utilized, we calculate the Ozmidov length scale , which provides an estimate for the smallest scale affected by buoyancy effects in stably stratified turbulent flows (also included in Table 1). In the previous expression, *ε* is the turbulence dissipation rate, calculated by applying Kolmogorov’s second hypothesis of similarity to the resolved energy spectra (e.g., Lundquist and Bariteau 2015; Muñoz-Esparza et al. 2018), and is the Brunt–Väisälä frequency. The U10C25 case is characterized by an Ozmidov scale of 8.5 m, which decreases as geostrophic forcing is weakened or cooling rate is increased. The smallest Ozmidov scale is 4.7 m, confirming that a vertical grid resolution of = 5 m is appropriate, while some reliance on the SGS model is expected.

While some recent LES studies of the SBL have utilized 5 m (e.g., Sullivan et al. 2016; = 0.39 m), Basu and Lacser (2017) have recently pointed out that when vertical grid spacing is lower than 50, the first grid point where MOST is applied lies well within the roughness sublayer (RSL). Placement of the first vertical at such close distance to the surface would require the use of additional RSL corrections in combination with the MOST surface boundary conditions, which are not a standard in surface-layer models. Violation of such condition can potentially lead to spurious reduction of mixing at the surface, as noted by Basu and Lacser (2017). In addition, several LES studies (e.g., Zhou and Chow 2011; Allaerts and Meyers 2018) have demonstrated that a vertical grid spacing of 5 m is appropriate to simulate weakly stratified ABLs (with similar forcing conditions to the ones we consider) when using a fourth-order energy-conserving finite difference scheme in the vertical direction, which has a comparable spectral resolution to the hybrid even–odd-order difference scheme used herein. Therefore, these aspects further support the choice of = 5 m in our simulations.

For convective conditions, four different ABLs were considered. To set up these cases, two geostrophic wind conditions ( = 5, 15 m s^{−1}) were combined with four surface kinematic heat fluxes ( 0.05, 0.1, 0.3, 1.0 K m s^{−1}) to span a large range of conditions, as seen in Table 2. The parameter determines the nature of the structures in the convective ABL. Typically, represents the transition from convective rolls to cells as increases (e.g., LeMone 1973; Grossman 1982; Weckwerth et al. 1999; Salesky et al. 2017). Therefore, two cases corresponding to the roll regime and two corresponding to the cell regime were simulated, with each case within a particular regime achieved by different strength of the forcing parameters (geostrophic wind and surface heat flux). In addition, Deardorff’s convective velocity scale and the surface temperature scale were calculated, since these will be used later in the analysis. The grid resolution for the convective ABL simulations was set to and m, since the integral length scales (order of ) are much larger than in the stable ABL, and a larger domain is needed to properly represent these coherent structures (, , = 54, 27, 2 km; with 200 grid points in the vertical direction). For the convective boundary layer (CBL) simulations, the parent domain runs for 2 h of physical time, since equilibration rates are much faster than in the stable ABL. In all the simulations (SBL and CBL) a Rayleigh damping layer with a coefficient of 0.1 s^{−1} is placed for to prevent turbulence fluctuations to reach the model top. Roughness length was set to 0.1 m, with a Coriolis forcing parameter 1 × 10^{−4} s^{−1}, representative of midlatitudes (43.3°N).

## 3. Turbulence generation in stably stratified ABLs

### a. Perturbation amplitude based on a modified Richardson number

One aspect of the perturbation Eckert number that complicates its application to mesoscale–LES downscaling of stable boundary layers in real conditions is the use of the geostrophic wind in its definition. While idealized cases are typically forced by a geostrophic wind constant with height (e.g., Beare et al. 2006), SBLs often develop under the presence of baroclinicity. These horizontal temperature gradients originating from fronts or differential heating lead to the development of geostrophic wind shear (i.e., geostrophic forcing varying with height). Moreover, low-level jets occurring during stable conditions can present a strong ageostrophic component when induced by inertial oscillation mechanism (e.g., Vanderwende et al. 2015). This requires estimation of large-scale pressure gradient from mesoscale results, which further complicates the derivation of an accurate yet simple parameterization of the potential temperature perturbation amplitude that works under realistic forcing conditions.

An appropriate determination of the amplitude of the perturbations in the stable ABL requires information about mechanical production and buoyancy destruction of turbulence. In a recent application of the cell perturbation method to an engineering-type channel flow configuration in the presence of no thermal forcing (neutral stratification) by Buckingham et al. (2017), the authors found that thermal perturbations satisfying a Richardson number of 0.25 provided development of turbulence within the shortest fetch. Although only one forcing scenario was considered, results provided similar performance as the optimum = 0.2 for atmospheric flows (fetch 10). Herein, we extend that approach and propose a modified Richardson number scaling, which accounts for both the mesoscale stabilizing buoyancy gradient and the instability generated by the stochastic perturbations in the following manner:

where , , and are the local vertical differences of zonal, meridional, and potential temperature over a grid size, calculated from the incoming mesoscale flow. In Eq. (1), the stochastic perturbation term has a negative sign, since its role is to induce thermal instability, competing with the stabilizing source of the stratified mesoscale potential temperature gradient.

In order for the thermal perturbations to result in an efficient destabilizing effect, the modified Richardson number should be < 0, implying local overturning of the stably stratified inflow conditions and therefore triggering of an accelerated transition to turbulence. In contrast to the generalized perturbation Eckert number for neutral ABLs, a constant scaling results in a distribution of perturbation amplitudes that varies with height, = , since the Richardson number is typically height dependent. For all the simulations with the CP method, we adjust 1, as proposed by Muñoz-Esparza et al. (2015), to allow the signature of the perturbations to be advected out the cell before the next set of stochastic perturbations is introduced.

To have an initial qualitative assessment of the performance of the modified Richardson number scaling, the five mesoscale-to-LES nested SBL cases were simulated using = −1.0 and compared to the reference unperturbed cases. Figure 2 shows horizontal planes of vertical velocity *w* at a representative height = 0.5, where is the height of the low-level jet maximum in each case. In all cases, the simulations using the cell perturbation yield an accelerated transition to turbulence, occurring within a short fetch, with a slight delay in the U10C25 case. On the contrary, in all situations, the unperturbed cases necessitate a fetch of 3–6 km to initiate turbulence development. It is worth remarking that simulations of real-world stable ABLs are typically performed over domains with horizontal extensions of only a few kilometers (e.g., Zhou and Chow 2014; Muñoz-Esparza et al. 2017), given the constraints imposed by the required grid spacing, and for which little to no turbulence would develop in the unperturbed cases.

A series of simulations with the CP method were carried out to span a range of values to quantify the impact of the perturbation amplitude on SBL turbulence generation. Figure 3 includes the streamwise evolution of vertical velocity variance for the unperturbed case, several values (= −0.2, −1.0, −2.0, and −5.0), and the optimum Eckert number of 0.2 at a height = 0.5. Regardless of the amplitude of the potential temperature perturbations introduced, the CP method accelerates the development of turbulence in the nested LES domain. As is decreased, the magnitude of increases, and the turbulence onset is accelerated. For sufficiently large values, an initial peak in starts to develop, increasing both the slope and its maximum value. This effect is clearly observed in all cases for = −5.0. While the generation of three-dimensional turbulence is considerably accelerated for that value, there is a subsequent larger departure from equilibrium that necessitates larger fetches for stabilization with the forcing to be achieved. Ideally, we would like turbulence onset and energy growth to occur as fast as possible, while reaching quasi equilibrium within the shortest fetch (i.e., minimum initial overshoot). In general, = −1.0 provides an efficient development of turbulence over the range of considered forcings. However, some departures from such are found to be more optimal in specific cases.

To assess the performance of the different -based potential temperature perturbation amplitudes, we determine the optimum as the value for which the solution stabilizes within ±5% of the quasi-equilibrium state. Note that the nested mesoscale-to-LES cases we simulate correspond to heterogeneous boundary layers, since the ABL is continuously evolving in response to the forcing conditions. Such heterogeneity precludes the use of periodic LES as reference. In addition, as shown by Muñoz-Esparza et al. (2014), imposing similar forcing conditions in periodic and coupled simulations is challenging. Instead, we found that linearly extrapolating in a backward manner from the northern and eastern boundaries of the LES domain provides sufficiently accurate criterion and clear distinction between the different cases. The optimum values are = −2.0 for U10C25, = −1.0 for U10C50 and U15C75, and = −0.2 for U10C75 and U5C25. The overall optimum value, considering all the cases, is = −1.0. It is worth noting that despite the turbulence onset observed in Fig. 2 for the unperturbed cases, the equilibration processes are delayed in comparison to when the CP method is applied, thus requiring fetches of 8–10 km.

For completeness, the results from = 0.2 are included in the comparison presented in Fig. 3, confirming the validity of the perturbation Eckert number scaling for stable ABLs. However, application to real-world cases includes determination of geostrophic wind and boundary layer height, making the derivation of the optimum perturbation amplitude difficult. This can lead, for instance, to larger perturbations than likely required, as noted by Muñoz-Esparza et al. (2017), in performing multiscale full-physics simulations of a nighttime period during the CWEX13 campaign. In general, the perturbations corresponding to = −5.0 are too strong, leading to deviations from the equilibrium state (i.e., more vigorous turbulence fluctuations). In contrast to the rest of the cases, U15C75, which features the largest geostrophic wind, results in the slowest equilibration rates, with the solutions from the different simulations revealing some variability at the end of the domain. We attribute this effect to the larger mean wind speed, compared to the turbulent counterpart, thus requiring larger distances for the turbulence to interact with the mean flow forcing and converge to the equilibrium solution.

Streamwise evolution of resolved turbulent kinetic energy , momentum flux , heat flux , and skewness and kurtosis (third- and fourth-order moments) of the vertical velocity and , where is the *j*th moment about the mean, are presented in Fig. 4 for the U10C25 case at 0.5 (comparable results are found for the other cases and heights; not shown). Spatial evolution of *q* and are similar to those of shown in Fig. 3, exhibiting a direct correlation between and the transitional behavior of the flow in the nested LES domain. The evolution of the heat flux reveals some interesting features regarding the transition mechanism with the cell perturbation method in stably stratified ABLs. At first, there is a region in which grows to a positive maximum value, followed by a decrease to zero. Afterward, the heat flux has a large drop, reaching a minimum that is coincident with the peak in *q* and , and of negative sign that represents the cooling flux in the stable ABL. For small perturbation amplitudes, −0.2, the strength of the perturbation-induced convective flux is not sufficient to significantly disrupt the quiescent mesoscale inflow, and although transition to turbulence is accelerated with respect to the unperturbed case, the required fetch is still significant. If the perturbations are too strong, as for = −5.0, the local convective flux is significant, but the ABL displays a “rebound” effect. This originates from the necessity to overcompensate for the effect of the strong perturbations, resulting in cooling heat fluxes larger in magnitude than the equilibrium solution that in turn delays the transition to a fully developed turbulent state. In general, the larger the positive convective heat flux, the larger the departure of the initial peak in TKE and fluxes from the quasi-equilibrium state. For this case (U10C25), = −2.0 is the optimum, exhibiting a fast yet smooth transition of in the absence of an overcooling heat flux peak.

Spatial development of skewness and kurtosis in Fig. 4 reveals the signature of the above-described mechanism. Higher-order moments appear to be more sensitive to the perturbations, as seen in the evolution of for = −2.0, which develops a peak, although it is followed by a rapid reequilibration. Such equilibration transitions from 0 are indicative of the convective instability to slight 0 to balance the previous positively skewed distribution, finally achieving 0. From the evolution of vertical velocity kurtosis, it is observed that the increase in the amplitude of the perturbation magnitude (decrease in ) results in larger values, indicative of a probability distribution with predominance of tails (large-magnitude *w* values), resulting from the strong perturbations applied to the resolved potential temperature field. Progressively, the distributions transition to 3.25, which bears a Gaussian distribution, as expected in nearly isotropic vertical fluctuations in the stably stratified ABL. Of note, for large perturbation amplitudes (i.e., = −5.0), the strong interaction of turbulence with the forcing results in slight changes to the equilibrium state, typically manifesting as enhanced turbulent fluctuations and momentum flux, while the cooling heat flux becomes smaller in magnitude as a result of the enhanced mixing.

### b. Energy spectra, turbulence length scale, and transition mechanism in the absence of inflow perturbations

To analyze the quality of the generated turbulence with the modified Richardson number-based perturbations, streamwise evolution of compensated spectra for the vertical velocity , potential temperature , and the compensated cospectra between the streamwise and vertical velocity components are shown in Fig. 5. Results correspond to the U10C25 case at = 0.5 for the unperturbed case (left panels) and CP method with = −2.0 (right panels). Spectra and cospectra are calculated along the spanwise direction using a Hamming window, with each line representing a 1-km streamwise separation starting at *x* = 0.5 km, and are averaged for 1 h physical time. The CP method simulation exhibits a signature from the imposed perturbations at *x* = 0.5 km, which quickly evolves into an energy distribution that is self-similar to the equilibrium state at *x* = 1.5 km and becomes fully developed by *x * 2.0–2.5 km, displaying the expected Kolmogorov’s −2/3 slope at the inertial range of three-dimensional turbulence. At the very early stages, and display a multipeak distribution, with the lowest minimum at 8. While this result may seem counterintuitive, given the fact that the cells have an 8 × 8 gridpoint footprint, numerical dissipation from the finite-difference discretization rapidly smooths out the edges of these sharp perturbations, resulting in a reduction in the effective scale of the perturbation to 5.5. The development within the inertial range points to a sort of triadic interaction of scales, resulting from the quadratic nonlinearity of the governing Navier–Stokes equations. However, full understanding of the turbulence interactions at these very close distances from the perturbation region is beyond the scope of the current study. The key point of the CP method remains the rapid development of the entire spectrum of turbulent motions to the quasi-equilibrium state. In contrast, the unperturbed simulation requires a fetch of 10 km to generate the proper energy distribution (a factor of 5 times longer than the CP method).

The energy cospectra reveal similar features, with a fully developed spectrum already at 2 km. Formation of mesoscale-like energy content for 3 × 10^{−4 }, mainly affecting both horizontal velocity components and temperature, is clearly observed in and . This large-scale variability arises from the differential cooling that results from the heterogeneous nature of the developing boundary layers simulated herein. The apparent large-scale variability is a consequence of the misalignment of the mean wind and the domain layout due to Coriolis forcing and is consistent with previous findings by Muñoz-Esparza et al. (2014). Spatial evolution of energy spectra of the vertical velocity at = 0.5 is presented in Fig. 6 for the rest of the stable ABL cases. The CP method using the optimum value to derive the amplitude of the potential temperature perturbations shows fully developed energy spectra in all cases by 2 km. In the unperturbed simulations, the development of turbulence occurs in all cases through energy production at wavenumbers in the neighborhood of the inverse of the integral length scale, which form a pronounced peak. Such peak often exceeds the equilibrium state, and it is responsible for the overshoot in vertical velocity variance observed in Fig. 3. The rate of energy growth at larger as well as inertial range scales is delayed, and the spectrum does not equilibrate until all the scales are present and the energy can cascade from the integral-scale peak through the resolved inertial range. Moreover, turbulence onset and its rate of development vary from case to case in a nonintuitive manner.

A closer inspection at the turbulent modes is provided by the quasi-equilibrium (*x* = 12.5 km, = 0.5) compensated spectra of vertical velocity corresponding to the optimum shown in Fig. 7. Production of three-dimensional turbulence emerges for wavelengths approaching *k * 2 × 10^{−2} m^{−1}, characterized by a 0 slope. Then, transition to −2/3 inertial range is observed for scales smaller than 50 m. Note that only a limited extent of the inertial range is fully resolved, since implicit filtering arising from the discretization of the advection term starts to remove energy for *k* ≤ 3 × 10^{−2} m^{−1}. This result corroborates the necessity of small grid spacings ( < 10 m) together with accurate numerical schemes and SGS models to represent the relevant turbulent motions in stably stratified ABLs. Also, it further supports that the selected grid spacing is sufficient to resolve production scales and some portion of the inertial range in all cases; therefore, our LES results used as base to propose a stability-aware scaling for the cell perturbation method capture the important small-scale turbulent instability mechanisms in the SBL.

To gain further insight into the turbulence transition process in the absence of perturbations, vertical cross sections of potential temperature at = 0.5 are presented in Fig. 8 for the U10C25, U10C75, and U15C75 cases. In both the U10C25 and U10C75 cases, the onset of turbulence occurs at a similar distance from the western lateral boundary, 3 km, and farther downstream for the U15C75 case, while the vertical location at which such instability appears differs in all cases. In the U10C25 case (Fig. 8a), the instability starts at 350 m; in the U10C25 case, it initiates at 225 m (Fig. 8c); and in the U15C75 case, it initiates at 400 m (Fig. 8c), always near the jet maxima. In all cases, turbulence onset is produced by a Kelvin–Helmholtz (KH) instability, as demonstrated from the structures in the potential temperature distribution (note that KH billows look almost vertically oriented due to the magnified scale of the vertical axis). Then, the instability propagates spatially, expanding in the vertical direction as the flow is advected downstream. The simulations using the CP method with the optimum (Figs. 8b,d,f) display a homogeneous development of turbulence across the entire ABL, which also takes place within a much reduced fetch. In particular, it is worth remarking that for the U15C75 case, a vertical layer of 200 m from the surface remains nearly unmodified from the smooth mesoscale inflow conditions. This will have a negative impact on applications where near-surface phenomena, such as wind energy, wildland fire, or pollutant dispersion, among others, are relevant.

Finally, an explanation for the variability of optimum modified Richardson number from −1.0 in some cases is proposed. The amplitude of the perturbations that provides the smoothest transition to turbulence without resulting in an overshoot near the inflow region and a slower adjustment to the equilibrium state is associated with the right amount of mixing to produce sufficient local overturning of the potential temperature field without having the stratified flow to overcool and restore the stable stratification of the ABL. This can be interpreted as the competition between turbulence kinetic energy and buoyancy suppression, for which a characteristic length scale can be defined. Figure 9a shows vertically integrated , , as a function of the optimum . Stable ABLs with low values are more dominated by buoyancy effects and, therefore, do not require large-amplitude perturbations (i.e., turbulent mixing) to optimally reach their equilibrium state. As becomes larger, the turbulence kinetic energy term dominates; therefore, the equilibrium states can withstand stronger mixing and, hence, reduced values (larger in magnitude). A good linear correlation is found between the optimum and , = 1.602 − 0.414 , with a coefficient of determination 0.93. This relationship can be used to determine the optimum . Examination of the vertical distribution of , shown in Fig. 9b, points to a critical value of 8 m that is favorable for the onset of KH instability in the unperturbed cases. For larger values of that are present at closer distances to the surface, the numerical noise is not strong enough to trigger the instability, while for reduced , buoyant stratification is dominant and prevents the growth of numerical noise.

## 4. Turbulence generation in convective ABLs

### a. Perturbation amplitude based on a thermal variance scaling

As for the stable ABL, an alternative scaling for convective ABLs is explored. The amplitude of the potential temperature perturbations can be derived from mixed layer scaling laws for temperature variance , which correspond to quasi-equilibrium conditions. Herein, we use the temperature variance scaling for convective boundary layers suggested by Sorbjan (1989):

From Eq. (2), and knowing and from the mesoscale inflow conditions, can be derived. Note that this expression accounts for increased near the ABL top as a result of entrainment of warmer air aloft. Nested mesoscale-to-LES simulations of the CBLs summarized in Table 2 were carried out, comparing the unperturbed case to the CP method using the perturbation Eckert number and the temperature variance scalings. Figure 10 shows spatial development of vertical velocity variance in the streamwise direction for the four CBL cases at a height of = 0.5. A variety of situations are found, from cases like U5H3, where perturbations have very reduced impact and turbulence development occurs within a relatively short fetch, to cases such as U15H2, in which the unperturbed simulation requires a fetch of 30 km and the application of the CP method with 0.2 reduces the fetch by a factor of 6 (5 km). However, in the latter case, the temperature variance scaling does not result in a significant reduction of the development distance, compared to the unperturbed simulation.

### b. Transition mechanisms in the absence of inflow perturbations

To elucidate the physical mechanisms responsible for the observed differences in turbulence generation and its equilibration, the temporal evolution of the potential temperature field is analyzed. Figure 11 displays potential temperature contours in a vertical plane at = 0.5 for the U5H3 case at *t* = 6, 8, 10, 15, and 20 min from the initiation of the nested LES domain. The thermal instability at the surface is rapidly triggered by numerical noise, likely induced by adjustments in the LES domain from the imposed LES solution, resulting in a mushroom-like updraft typical from Rayleigh–Bénard convection as seen at *t* = 6 min. This upward motion coherently develops along the entire spanwise direction and propagates the instability downstream through advection and pressure terms. Development of vertical updrafts occurs downstream of the initial instability, *t* = 8, 10 min. By *t* = 15 min, some of these eddies have reached the top of the ABL, exhibiting local entrainment of inversion layer air parcels. At *t* = 20 min, the flow has completely reached equilibrium, and potential temperature field presents a well-mixed structure. The bottom panel in Fig. 11 shows that while the CP method has a beneficial effect in accelerating the instability process within the entire ABL, the transition to turbulence is dominated by the existence of surface instability and its vertical propagation, further confirming the negligible differences in the spatial development of observed in Fig. 10.

In the U15H2 case, notable differences are observed (see Fig. 12) with respect to the U5H3 case. In spite of the onset of thermal instability originating at the surface at a close distance to lateral boundary of the nested LES domain, the growth of the instability is now delayed in space and remains in the near-surface region, *t* = 18 min. Organized updraft motions start to develop, progressively losing coherence and giving birth to smaller-scale eddies (*t* = 20, 22, 24 min). By the time the flow has reached its quasi-equilibrium state, *t* = 30 min, it can be clearly seen that the thermal instability at the surface does not lead to the onset of turbulence until 7 km. Even at this location, an additional streamwise distance of 7–8 km is required for the instability to propagate across the ABL. In contrast, the CP method helps the thermal instability to propagate by triggering flow fluctuations within the entire ABL and results in a significant acceleration of turbulence development and equilibration, as shown in Fig. 10.

Analysis from the CBL cases clearly indicates a competition between the thermal instability that attempts to propagate vertically and the horizontal wind that tries to advect the instability in the downstream direction. Therefore, we propose the ratio as the parameter that primarily indicates the dominant transition mechanism, where is the wind speed at the top the capping inversion = + , with the thickness of the capping inversion. For similar values, we find that rolls tend to reach equilibrium more rapidly, compared to cells. This can be attributed to the anisotropy of the former, aligned with the wind direction (cf. U5H1 and U15H4 in Fig. 10, which have = 12 and 148, respectively). This parameter can then be used to correct the amplitude of the perturbations , leading to improvements in all cases (see Fig. 10). The use of the CP method in CBLs will be especially beneficial in cases where 5, when the thermal instability by itself is not sufficient to result in an efficient vertical propagation throughout the entire ABL. For illustration purposes, Fig. 13 shows the spatial development of turbulence kinetic energy, momentum and heat fluxes, and higher-order moments for the U15H2 case. The use of the corrected temperature variance scaling results in a rapid transition and stabilization to the equilibrium *q*. Momentum fluxes require slightly larger distances to stabilize, most likely associated with the mean flow transition from the mesoscale to the LES solutions. Nevertheless, heat flux, skewness, and kurtosis of all the velocity components (only *w* is shown) rapidly become fully developed, with the unperturbed simulation requiring significantly larger fetches.

It is worth mentioning that the perturbation Eckert number proposed by Muñoz-Esparza et al. (2015) performs very well in all convective ABL scenarios analyzed herein, with only minor differences with respect to the scaling. We attribute this behavior to the fact that the amplitude of thermal perturbations is related to the geostrophic wind, which is one of the two competing mechanisms in the transition process for convective ABLs. The scaling uses the wind speed at the top of the capping inversion, which is almost equal to the geostrophic wind for the idealized simulations performed herein. While the use of the ABL height may be problematic in stable boundary layers, convective ABLs are typically limited by a well-defined capping inversion, which makes the determination of an accurate , and in turn , much simpler. Finally, in situations where 3, the use of the CP method does not provide a significant advantage; nevertheless, it is still beneficial since it accelerates the transition in the presence of underresolved convective structures that emerge in high-resolution mesoscale simulations (Mazzaro et al. 2017).

Finally, to illustrate the dependence of the amplitude and variability of the potential temperature perturbations within the ABL upon atmospheric stability, the vertical distributions of for all the cases simulated herein are presented in Fig. 14. Overall, the maximum amplitude of the perturbations during stable conditions is <1 K, with a trend to increase in the lowest few meters above the ground, where the shear production is larger. During convective conditions, is also typically 1 K in the bulk of the mixed layer. However, in the CBL, there is a significant increase in thermal variance toward the surface, which in turn augments the strength of the perturbations. Note that the U15H2 and U15H4 cases result in larger values within the mixed layer. This is mainly originated by the correction term in the scaling , which for these cases is 10 and 5, respectively. However, an enhanced thermal instability is required to trigger an effective onset and transition to a turbulent state. Strong perturbations at the surface may create some spurious response from surface physics under some specific conditions, including surface-layer parameterization and land surface model. Therefore, when targeting full-physics simulations, we recommend leaving the first three levels unperturbed, as in Muñoz-Esparza et al. (2017), to minimize any undesired feedback from the potential temperature perturbations. Nevertheless, note that all these surface parameterizations are one-dimensional, and thus its impact is likely to remain within the perturbation region, where model results should not be used anyway.

## 5. Conclusions

The generation of forcing-consistent turbulence in nested mesoscale-to-LES simulations of atmospheric boundary layers is one of the challenges hindering seamless multiscale modeling capabilities. Herein, we have extended our previous work on the cell perturbation method for nonneutral ABLs (stable and convective) and provided insight into the transition mechanisms. This is a relevant aspect, since in most cases the atmospheric boundary layer is influenced by stability effects, with neutral conditions mostly occurring on inland locations during the morning and evening transitions of the diurnal cycle. To that end, a set of nested LES within a mesoscale flow simulation has been performed, carefully designing the experiments to cover a broad range of the conditions that can be simulated using LES of ABLs with large yet affordable computing capabilities.

A modified Richardson number scaling is proposed for the stable ABL, which accounts for buoyancy suppression from the mesoscale inflow (i.e., positive vertical potential temperature gradient) and the instability effect of the thermal perturbations introduced under the CP method. Generally, it is found that −1.0 results in the optimum transition to fully developed turbulence, where equilibrium is reached within the shortest fetch. The optimal value of is determined through analysis of spatial evolutions of vertical velocity variance, turbulence kinetic energy, momentum and heat fluxes, skewness, kurtosis, and energy spectra and cospectra distributions. However, deviations from the optimum value are encountered. For further accuracy, a turbulent length scale , which is the ratio of the square root of turbulence kinetic energy to the Brunt–Väisälä frequency, is proposed to account for this variability, showing a good degree of correlation with the optimum . In the absence of perturbations, the transition mechanism is a Kelvin–Helmholtz instability, typically occurring near the low-level jet maximum, where 8 m. Then, the instability progressively spreads in the vertical direction as the flow is advected within the LES domain. The optimum transition with the CP method is the result of sufficient potential temperature amplitude that can generate moderate local convective overturning. This overturning leads to development fetches of 2 km, while the unperturbed cases necessitate fetches that are 4–5 times larger.

In the context of convective ABLs, perturbation amplitudes based on scaling of potential temperature variance in the mixed layer are proposed: . However, these criteria were initially not found to be optimum. By analyzing the transition mechanisms in convective conditions, we found that in spite of the thermal instability that is always present, larger variations of the required fetch to reach equilibrium are observed for the unperturbed cases. It is identified that there are two competing mechanisms: thermal instability at the surface and horizontal advection. Therefore, we propose the ratio of the wind speed at the top of the capping inversion to the convective velocity scale as the characteristic parameter. By using as scaling for the temperature perturbations, optimum transition to turbulence is achieved. While in cases where 3, the use of the CP method does not provide a significant advantage, the development fetch was reduced by a factor of 6 when 10. Moreover, application of the CP method has proven beneficial in the presence of mesoscale underresolved convective structures.

With the present work, we have provided understanding of the turbulence transition mechanisms in large-eddy simulations nested within nonneutral mesoscale-driven ABLs. Furthermore, we have presented stability-aware scalings of the optimum amplitude of the potential temperature perturbations with the CP method. While recent application of the CP method for a diurnal cycle by Muñoz-Esparza et al. (2017) has shown promise in enabling realistic multiscale modeling of ABL flows, we expect the current extensions of the method to further optimize the inflow turbulence generation aspect of coupled mesoscale–LES. We intend this contribution to enable more realistic multiscale simulations of flows in atmospheric models. We plan to provide the improvements developed herein to the community release version of WRF in the near term.

## Acknowledgments

We would like to acknowledge high-performance computing support from Cheyenne (https://doi.org/10.5065/D6RX99HX) provided by NCAR’s Computational and Information Systems Laboratory. The National Center for Atmospheric Research is sponsored by the National Science Foundation. Computing time was provided by NCAR Strategic Capability Computing Grant NRAL0016. The authors are grateful to Jeremy Sauer and Robert Sharman for providing comments on an earlier version of the manuscript, as well as to two anonymous reviewers for their constructive comments.

## REFERENCES

*22nd Symp. on Boundary Layers and Turbulence*, Salt Lake City, UT, Amer. Meteor. Soc., 3B.5, https://ams.confex.com/ams/32AgF22BLT3BG/webprogram/Paper295892.html.

*Structure of the Atmospheric Boundary Layer*. Prentice Hall, 317 pp.

## Footnotes

^{a}

The National Center for Atmospheric Research is sponsored by the National Science Foundation.

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

^{a}

The National Center for Atmospheric Research is sponsored by the National Science Foundation.