## Abstract

Coupled mesoscale–microscale simulations are required to provide time-varying weather-dependent inflow and forcing for large-eddy simulations under general flow conditions. Such coupling necessarily spans a wide range of spatial scales (i.e., ~10 m to ~10 km). Herein, we use simulations that involve multiple nested domains with horizontal grid spacings in the terra incognita (i.e., km) that may affect simulated conditions in both the outer and inner domains. We examine the impact on simulated wind speed and turbulence associated with forcing provided by a terrain with grid spacing in the terra incognita. We perform a suite of simulations that use combinations of varying horizontal grid spacings and turbulence parameterization/modeling using the Weather Research and Forecasting (WRF) Model using a combination of planetary boundary layer (PBL) and large-eddy simulation subgrid-scale (LES-SGS) models. The results are analyzed in terms of spectral energy, turbulence kinetic energy, and proper orthogonal decomposition (POD) energy. The results show that the output from the microscale domain depends on the type of turbulence model (e.g., PBL or LES-SGS model) used for a given horizontal grid spacing but is independent of the horizontal grid spacing and turbulence modeling of the parent domain. Simulation using a single domain produced less POD energy in the first few modes compared to a coupled simulation (one-way nesting) for similar horizontal grid spacing, which highlights that coupled simulations are required to accurately pass the mesoscale features into the microscale domain.

## 1. Introduction

Accurate simulation of the turbulent flow field in the atmospheric boundary layer requires high resolution simulations as well as representation of the impacts of the large-scale features of the flow. One approach is to use coupled mesoscale and microscale simulations, which resolves a wide range of scales of atmospheric motion. Because of the large horizontal footprints of traditional mesoscale grid cells (horizontal grid spacing ~10 km), the effects of turbulence are parameterized in the vertical direction only [e.g., Mellor–Yamada–Nakanishi–Niino level 2.5 (MYNN; Nakanishi and Niino 2006) and Yonsei University (YSU; Hong et al. 2006) schemes] by assuming horizontal homogeneity. Along the horizontal direction, diffusion based on 2D Smagorinsky closure is applied to treat the turbulence. In contrast, for turbulence-resolving simulations (horizontal grid spacing 100 m), turbulence models must parameterize the turbulence in all three orthogonal directions [e.g., Smagorinsky method (Smagorinsky 1963), Lilly method (Lilly 1967)] down into the subinertial range scales. Between the limits of these length scales (i.e., microscale ~100 m to mesoscale ~1 km), parameterization/modeling of turbulence is not represented well to simulate realistic turbulence in either mesoscale or microscale simulations. This range of spatial scales for which existing turbulence models are expected to fail was termed as the “terra incognita” by Wyngaard (2004).

In coupled simulations, the terra incognita may appear in the cascade of the scales of motion and could impact the solution on nested microscale domains. In our previous work using a mesoscale model (Rai et al. 2017a), velocity time series for simulations using the MYNN boundary layer parameterization with the horizontal grid spacing in the terra incognita showed an unrealistic oscillatory behavior. Ching et al. (2014) also documented secondary structures in their simulation for subkilometer horizontal grid spacing. Several others have reported the impact of grid spacing in the terra incognita on their results. For example, Shin and Dudhia (2016) studied the terra incognita using different PBL schemes and found that none of the schemes are scale aware and superior to others. Similarly, Honnert et al. (2011) evaluated the subgrid and resolved portion of variables [e.g., turbulent kinetic energy (TKE), heat and moisture flux, etc.] at different subkilometer scales and found that the size of the structures of the velocity, temperature, and mixing ratio are different. Shin and Hong (2013) studied the resolved and parameterized vertical transport dependence on grid spacing in the terra incognita and found that the resolved TKE and vertical flux increases with the larger wind shear relative to the total TKE and fluxes for a given grid spacing. Zhou et al. (2014) studied the behavior of a single domain model (using both a PBL scheme and an LES model) in the terra incognita and pointed out physical and dynamical problems associated with the terra incognita. Efstathiou et al. (2016) used blending of the turbulence model and a bounding of vertical diffusion coefficient approach to study the impact of horizontal grid spacing on the mean profiles and TKE partitioning during the development of the boundary layer in the morning. Both blending of the turbulence model and a bounding approach showed weaknesses and strengths in improving the results. In recent year, Zhang et al. (2018) and Kurowski and Teixeira (2018) developed the scale-adaptive TKE closure to parameterize the turbulence across the grid resolutions—from large scale to microscale including the terra incognita. All of these studies focused on the model outputs, particularly, within the terra incognita, but they did not evaluate the impact of terra incognita on subsequent nested microscale domains.

Relatively few studies (e.g., Rai et al. 2017a; Muñoz-Esparza et al. 2017) have employed the one-way “online” coupling approach in Weather Research and Forecasting (WRF) Model framework (Skamarock et al. 2008) that combines both mesoscale and microscale domains. A large number of coupled simulations have used the one-way ‘offline’ coupling approach. For example, Talbot et al. (2012) used the WRF model to simulate hydrometeorological variables. Recently, Rodrigo et al. (2016) employed profiles of advection and temperature tendencies obtained from the WRF mesoscale domain to provide mesoscale information to an offline microscale simulation. These works accounted for the use of coupled simulations, but they did not account for the details of the impact of the terra incognita on their results.

The present work utilizes a suite of simulations using one-way online nested domains with horizontal resolution within the terra incognita to investigate the impact of the mesoscale boundary forcing on the microscale domain as diagnosed by the spectral energy of the resolved flow field. The results indicate that the amount of resolved turbulence and the structure of the turbulence simulated on the microscale domain primarily depends on the type of turbulence modeling used for the microscale domain itself and that the solutions are insensitive to the size of the horizontal grid spacing of the parent domain (i.e., mesoscale or within the terra incognita limit). Section 2 provides details of model configurations and physics options used for the simulations. In section 3, a brief description of the analysis tools, such as autospectral density function (ASDF) and proper orthogonal decomposition (POD), are presented. Section 4 presents a discussion of the results in terms of spectra, turbulence kinetic energy, and POD mode energy that are used to investigate turbulence properties of the simulation, followed by conclusions and summary in section 5.

## 2. Model configuration

A model configuration consisting of a parent domain (D01) and two sequentially embedded nests (D02 and D03) as shown in Fig. 1 were used for performing simulations using the WRF Model (version 3.7.1). Nested domains were used to downscale the horizontal grid spacing () from a length comparable to the midafternoon boundary layer depth () to the order of tens of meters. Domains D02 and D03 were shifted to the north to allow a large upstream fetch from the inflow boundary (i.e., south and east side) that provides sufficient space to spin up boundary layer turbulence, which we term as the “spinup distance” herein. Although the grid size of two domains D01 and D02 is fixed, the domain size and grid spacing of the innermost domain D03 varies in size (Fig. 1, right)—the larger domain has a grid spacing of 0.24 km and the smaller domain a of 0.04 km. The domains cover an area of the southern Great Plains that includes the Scaled Wind Farm Technology (SWiFT) site near Lubbock, Texas (Kelley and Ennis 2016). The terrain on the inner domain is relatively flat compared to the outermost domain due to the outer domains including the eastern portion of the Rocky Mountains. This SWiFT facility includes the Texas Tech University tower (black solid dot in Fig. 1, right), which collects data using 10 sonic anemometers at heights ranging from 0.9 to 200 m above the surface. Simulations were performed for three cloud-free days with varying midafternoon . These three cloud-free days provide three different boundary layer depths with simple meteorology and impact of other physical forcing. Each simulation was initialized at 1500 UTC (0900 CST) and run for 9 h. However, only the data from the period of 1400 to 1700 CST are used in the analysis as that period corresponds to a well-developed and quasi-stationary convective boundary layer. The model output was saved every 10 s for domain D01 and D02, and every 10 or 3 s from the domain D03.

A number of simulations with horizontal grid spacing (), within mesoscale, terra incognita, and microscale domains (i.e., for cases Ti1–Ti7) were performed. Table 1 provides the configuration details for and domain sizes of the simulations. There are three groups of : (i) 2.2–3.84 km (, mesoscale) for domain D01, (ii) 0.28–1.44 km (, terra incognita) for domain D02, and (iii) 0.04, 0.20, and 0.24 km (~microscale) for domains D02 and D03. The initial condition for domains D01–D03 and the lateral boundary forcing for domain D01 were supplied from North American Regional Reanalysis (NARR) (Mesinger et al. 2006) data at 3-h intervals and 32-km horizontal grid spacing. Model output from domain D01, in turn, provides the lateral boundary forcing to the inner domains D02, and so on. Note that, in all nested runs (i.e., Ti4–Ti7), the geographic size of all domains remains similar (although not identical), and only changes. Use of similar domain size and lateral boundary forcing but varying in simulations provides better comparisons of the model behavior associated with changes in . In all nested domain simulations (e.g., Ti1–Ti7), the grid refinement ratio (GRR) ranging from 2 to 11 has been used for downsizing the horizontal grid spacing from the parent (i.e., outer) to the child (i.e., inner) domains. In addition, simulations R1–R3 were performed using a single domain for comparison with the results from the nested simulations in a configuration that is similar to driving the microscale simulation with only large-scale forcing. It is noted that the domain sizes of R1 and R2–R3 are equal to the domain size of D03 and D01 of the cases Ti1–Ti7, respectively. Each simulation used the same stretched vertical grid spacing—from 10 m at the surface to 90 m at an altitude of 3 km and then to 950 m at the top of the model. The extent of vertical stretching for altitudes smaller than allows the grid aspect ratio (i.e., ratio of vertical to horizontal grid spacing) to be less than or approximately equal to 1 near the surface.

The simulations for D01 use a PBL parameterization [e.g., MYNN 2.5-level (Nakanishi and Niino 2006) or YSU (Hong et al. 2006)] schemes to model turbulence for cases Ti1–Ti7. However, the two inner domains, D02 and D03 use either a standard PBL parameterization (e.g., MYNN or YSU) or three-dimensional subgrid parameterization [1.5-order TKE closure (Lilly 1967)] or the combination of both. Finer grid spacing can resolve more turbulence irrespective of the turbulence model. The grid spacing and type of turbulence model in the simulation are application specific. For instance, a cloud-resolving simulation [e.g., convective updraft (Fan et al. 2017)] uses relatively smaller grid spacing (1 km) with a standard PBL parameterization. However, the frequencies associated with resolved turbulence could depend on the type of turbulence model used. Thus, using a PBL scheme in the innermost domain allows us to compare the turbulence statistics resulting from that scheme with those from the three-dimensional subgrid parameterization. Details of other physics schemes and dynamics options used in the simulations are shown in Table 2. In subsequent discussions, the PBL parameterization with MYNN and YSU PBL schemes are termed as “Myn” and “Ysu” schemes, whereas the three-dimensional turbulence model is called the “Lilly” model. The naming convention for each simulation reflects the combinations of run name (e.g., Ti6), horizontal grid spacing on the innermost domain (e.g., 0.04 km), and turbulence parameterization/modeling schemes (e.g., Myn, Lilly). For instance, name “Ti6–0.04km(Myn–Lilly–Lilly)” represents case Ti6 for coupled simulations consisting of three nested domains with inner domain 0.04 km, and Myn parameterization on D01, and the Lilly model on D02 and D03.

## 3. Method of data analysis

A number of different second-order statistical moments (e.g., variance and TKE) were used to assess the behavior of the turbulent flow on the microscale domain when using a parent domain with grid spacing in the terra incognita. In this section, the methods for analysis of the turbulent flow, such as an ASDF and POD energy are described in detail.

### a. Autospectral density function

The ASDF was calculated either using time series from a given location (i.e., one-dimensional spectra in frequency space) or using two-dimensional data from a horizontal plane at a given instance in time (i.e., two-dimensional spectra in wavenumber space). To calculate the one-dimensional spectra, a time series of a given variable was divided into a number of blocks with equal lengths and each block was then windowed with a Welch method (Welch 1967). For each *k*th block time series, the Fourier transform was applied to obtain the spectral energy density as

where is the complex conjugate of , and *T* is the period of each block. Blocks of spectra were averaged arithmetically to obtain the mean spectra. Similarly, to compute the ASDF in two-dimensional wavenumber space from the planar data, a Fourier transform was applied to the planar data obtained from the seventh model layer (i.e., 95 m above the surface). As in one-dimensional spectra, a Welch window was applied to the planar data before applying a Fourier transform. The spectral energy density of the planar data at a given instant of time was obtained as

where, is the complex conjugate of . Here, and represent the wavenumber for the *x* and *y* directions, respectively. The two-dimensional spectra for a given instant of time were smoothed by averaging the individual spectra of the entire analysis time. A detail procedure for obtaining spectra can be found in Bendat and Piersol (2011).

### b. Proper orthogonal decomposition

POD (also known as orthogonal eigenfunction, principal component analysis, etc.) is a powerful mathematical tool for determining the partitioning of energy among various structures occurring within. The formulation and application of POD is described more thoroughly in Lumley (1967), Sirovich (1987), Berkooz et al. (1993), and Smith et al. (2005), therefore only a short description is presented here. To determine the energetic structures, POD identifies an optimal basis function (*ϕ*) from the turbulent flow field by applying maxima to the functional . The symbol *ν* represents the variable associated with the flow field, whereas the operators , , and signify the operations for time-averaging, scalar, or inner product, and norms. The calculus of variation to the above functional leads to an eigenvalue problem for the interval that can be defined as

where *ϕ* and *λ* are called eigenvector (basis function) and eigenvalues, respectively. The Kernel is given by correlation matrix , where *z* and *t* are the independent variables in space and time.

In the present work, we have applied the direct (i.e., classical) method of POD (Lumley 1967; Smith et al. 2005) to time series of vertical profiles of wind components to obtain the spatial POD modes and turbulence kinetic energy per unit mass (*λ*) for each spatial POD mode . Shah and Bou-Zeid (2014) used snapshot POD to study the behavior of large-scale structures in the atmospheric boundary flow. In this work, there are POD spatial modes for grid points along the *z* direction and wind components. Solving the eigenvalue problem [i.e., Eq. (3)] provides eigenvalues (*λ*), which is the estimate of kinetic energy per unit mass (i.e., *λ*) of the flow field.

The value of *λ* decreases as the number of spatial modes increases. The number of eigenvalues needed to represent the kinetic energy per unit mass of the turbulent flow field depends on amount different coherent structures of different length and time scales in the turbulent flow field. In general, the original flow fields in our simulations can be reconstructed using the first few most energetic POD modes and coefficients. Herein, we use only the quantity *λ* to evaluate the energy in each spatial mode of different simulation cases studies simulated here.

## 4. Results and discussion

This section provides a detailed discussion of results obtained by applying the analysis tools described in section 3 to the simulated data. The observed data are also used to evaluate the simulated data using their POD mode energy. Model results are then discussed briefly in terms of variances, ASDF, and POD mode energy to verify the quality of the simulation. Results are discussed using a range of factors: boundary layer depth, horizontal grid spacing, grid refinement ratio, spinup distance, turbulence parameterization/modeling schemes, and number of nested domains.

### a. Model evaluation

To obtain the simulated data, multiple nested domains were used that downscale the horizontal grid spacing from mesoscale domain (kilometer) to microscale domain (meter). Of these two extreme values of , the outermost domain used PBL schemes for parameterization of turbulence in mesoscale mode, whereas the innermost domain used the three-dimensional turbulence modeling (i.e., Lilly model) in microscale LES mode. In coupled simulations, we could use intermediate domains within the terra incognita to avoid large impractical grid refinement ratios. Subsequent sections present the analyses of the model output from these three limits—mesoscale, terra incognita, and microscale.

We analyze observed data to evaluate the simulated data in terms of velocity and POD spatial mode. The time–height contour plot constructed from the simulated and observed wind speed at SWiFT tower for the time period 1000–1200 CST (3 May 2014) is shown in Fig. 2. It should be noted that this figure is from a case in our previous study (Rai et al. 2017b) and is not listed in Table 1. More details on the model configurations and analyses of TKE and spectra for this data can be obtained from Rai et al. (2017b). In Fig. 2a, the first and second row (from the top) are a time–height plot of wind speed for the simulated and observed data (at 1-Hz frequency) for the first 0.2 km above the surface. The simulated data were obtained from the innermost domain (with 0.04 km) that used the Lilly turbulence model. The time–height plot for simulated data shows a slight underprediction of wind speed near the surface ( 50 m) compared to the observed data. Results from both simulated and observed time–height contours reveal similar details of turbulence structures. A quantitative comparison is conducted using the POD method. Results from the POD energy analysis provide useful information on the velocity correlations across multiple points. POD energy for the first 20 spatial modes (Fig. 2b) were computed using time–height data of the three velocity components. The first two modes from the simulated data underpredicted the POD energy compared to that from observed data. POD energy of spatial modes from both simulated and observed data show a close agreement. However, for higher POD spatial modes ( 3), the simulated data show a slight underprediction of POD energy compared to the observed data. These POD energy differences indicate that the simulations underpredict the energy in the higher modes. One reason for the larger POD energy found in observed data compared to simulated data for the higher POD modes could be due to the disparate effective resolution, which is 7 for the simulated data (Skamarock 2004) and 2 for the observed data (due to Nyquist frequency). Cumulative POD energy shown in the inset of Fig. 2b depicts the behavior of energy distribution in another way over the spatial modes—the first three modes of both observed and simulated data contains of the POD energy, and the POD energy of the observed data is spread across all spatial modes.

### b. Mesoscale limit

Simulations with horizontal grid spacing within the mesoscale limit parameterize all scales of motions smaller than the grid spacing. To see the results from the mesoscale limit, time–height contours and time series of wind speed at the location of the SWiFT tower are plotted (Fig. 3) for three different simulations (i.e., R2, R3a, and Ti6) for three different dates (see in Table 1). Note that each simulation uses the Myn PBL scheme for parameterizing the turbulence in the vertical direction and has a different boundary layer depth ()—varying from 1.6 to 3.2 km, which we will show as a length scale to help generalize our results. The boundary layer depth presented here is diagnosed from the standard output of the mesoscale model. In Fig. 3a, columns from far left to right show panels of time–height contours for decreasing horizontal grid spacing . Similarly, rows from top to bottom display the wind speed contours for cases with decreasing values of . Note that the wind speed seen in the time–height plot is directly associated with the scales resolved by WRF (i.e., 7 ). The contours plot along the diagonal panels are from the simulations whose is comparable to . In these cases, the time–height diagram shows no fluctuations in the wind speed within the boundary layer. Similarly, time series data taken from 0.5 km above the surface at the center of the domain (Fig. 3b) for these cases clearly show smoothly varying wind speed. However, there are secondary structures present within the boundary layer in the off-diagonal panels for the cases where . This model behavior clearly indicates that the mesoscale simulation is sensitive to /, and the mesoscale limit begins at . Results for off-diagonal panels are discussed in the following subsection.

The spectral representation of a given velocity component estimates the frequencies of the resolved turbulence structures. Figures 4a and 4b present the spectra of the *u*- and *w*-velocity components for simulations with varying and Myn/Ysu PBL parameterization within the mesoscale limit. All spectra were computed from time series data from 95 m above the surface for the time period 1400–1700 CST (i.e., three blocks of data, each block of 1 h) using the method described in section 3. Note that the simulated data were obtained from domain D01 for all cases and the grid spacing is approximately equal to (2.4 km). The results show that spectra of the *u*-component velocity for all cases are identical and there are no turbulence structures, indicating that the spectral energy is only associated with mesoscale features. Spectra for the *w*-component velocity for varying cases show some variations in their spectral energy, but the energy level is much smaller compared to that energy of the *u*-component velocity (five orders of magnitude less). This indicates minimal or no turbulence structures in the spectra, which is also supported by the results shown in Fig. 3 because those time–height and time series plots from cases with comparable to displayed no oscillations.

POD energy derived from data at multiple locations of all velocity components may provide a better representation of energy contained in the turbulent flow than the analysis of spectral energy derived from data at single point. The POD energy for different spatial modes presented in this section was computed using time series of a vertical profile of three velocity components up to 5 km from the surface for the time period 1400–1700 CST using the method discussed in section 3. We have analyzed the simulated wind speed/direction, heat flux, and boundary layer depth for this 3-h time period. The statistics shows that the first- and second-order moments remain similar along horizontal direction during this time period. This shows that flow to be quasi-stationary and the POD can be used to analyze the simulated data. Figure 4c shows POD energy for the first 21 spatial modes for cases that use Myn- and Ysu-PBL schemes and . That plot reveals that POD energy in all modes for all cases are nearly identical, except for a small deviation in the first mode. The similarity of POD energy for all cases results from the fact that the size of is approximately equal to or greater than . This indicates that the energy containing scales of motion scale with , and 7 is within the mesoscale limit. Looking at POD energy in the first few modes, the first mode alone can represent approximately 99% of the energy of the flow, whereas the first three modes can represent almost 100% of the energy of the flow. In addition to the vertical profile, an inflow plane was also used to compute POD energy from velocity components. The results show similar POD energies (not shown here) as that from a vertical profile when normalized by the first POD mode energy. The segregation of such high POD energy into the first three modes indicates that the flow structure is dominated by the large-scale flow (mesoscale).

### c. Terra incognita

As discussed in the previous section (see Fig. 3), the time–height contours of wind speed begin to show secondary structures for grid spacing smaller than . Time series of wind speed taken from 0.5 km above the surface also produces oscillation for smaller than (i.e., in off-diagonal panels). This behavior in the time series of wind speed was also observed in our previous study (Rai et al. 2017a). The PBL parameterization that is developed for the domains in the mesoscale limit (i.e., ) may not work as because the parameterization is based on estimating an ensemble averages of eddies within a grid box. For smaller grid spacing, the model may begin to partially resolve large turbulent structures. These results indicate that domains with horizontal grid spacing smaller than falls into the terra incognita.

Contour plots of instantaneous velocity provide qualitative information regarding turbulence structures in the flow as well as the quality of the simulated flow field. Figure 5 depicts the contour plots of the *w* component of velocity extracted from the seventh model layer (i.e., 95 m above the surface) at two instants in time—1400 and 1600 CST. It should be noted that all cases shown in Fig. 5 lie in the terra incognita as their is smaller than . In this study, even if the PBL schemes and LES closures are ill-defined at the horizontal resolution used, we have analyzed the simulated data in terms of flow visualization, second-order statistical moment, and spectral analysis in varying grid resolution in the terra incognita to help us better understand the representation of boundary layer turbulence using the various closure methods. In coupled mesoscale–microscale simulations, grid spacing of one and/or more domains fall into the terra incognita while keeping the grid refinement ratio within moderate value (e.g., 10). This study investigates the effect of grid spacing of the parent domain within the terra incognita on the microscale domain using both Lilly model and PBL schemes. Panels in Fig. 5a are a snapshot of the *w* component of velocity from the simulation cases that use Myn- or Ysu-PBL schemes for varying [i.e., from 1.44 km (top) to 0.32 km (bottom)]. The snapshots at 1600 CST are a telescopic view of a smaller area (see square in panels at 1400 CST) shown to present the flow structures more clearly. The snapshots from both instances reveal that finer structures are captured as is decreased (panels from top to bottom). Cases that use different PBL schemes but the same [i.e., Ti2–0.96km (Myn–Myn)/(Ysu–Ysu): rows 2–3] predicted different magnitudes of the *w*-component velocity. The impact of changing can also be seen on the spinup distance from the inflow boundary, which is reduced as becomes smaller. Figure 5b shows the contour snapshots obtained for the same domain configuration and velocity component as Fig. 5a but using the Lilly model for modeling subgrid-scale turbulence. Again, contour panels in the first and second columns represent snapshots of velocity for 1400 and 1600 CST. As for the previous velocity contours, the results show that, as expected, the turbulence structures became progressively more refined as the value of becomes smaller. In all cases, the flow structures modeled with the Lilly model in the innermost domain show similar structures—elongated, with a streak-like orientation along the mean wind direction (i.e., southeasterly). These streak-like structures do not occur in the entire boundary layer depth, rather they occupy the first 15% of . In addition, the wind speed and surface heat flux during the simulation time are relatively stronger (i.e., 9 m s^{−1} and 375 W m^{−2}), resulting in the stability parameter for these cases to be 40, where *L* is the Obukhov length. Usually horizontal rolls form in the moderately convective condition () (Weckwerth et al. 1999; Salesky et al. 2017). It implies that the streak-like structures herein are resulted most likely from fast and slow velocity streaks in the surface layer and/or the artifacts translated from the parent mesoscale domain. Compared to the cases that use mesoscale PBL schemes, the flow structures for cases that use the Lilly model are found to be completely different, primarily, in terms of their shapes and orientations. More about the shape and orientation seen in these flows will be discussed in the “microscale limit” section using two-dimensional spectra.

As discussed previously, domain D02 in the terra incognita begins to resolve scales of motions that are smaller than the horizontal grid spacing . To evaluate the model output in the terra incognita, spectra of *u* and *w* components of velocity from the simulations that use Myn/Ysu-PBL schemes with varying are shown in Fig. 6a (*u* spectra) and Fig. 6b (*w* spectra). The grid spacing of domain D02 varies from 0.28 to 1.44 km (all are within the terra incognita) and the domain was forced by domain D01 that uses the Myn- or Ysu-PBL scheme. The spectra show a continuous increase of resolved turbulence for decreasing , from a meager amount in the case Ti1 to a maximum amount in case Ti5. Between the two PBL schemes, the Myn scheme simulations results in more TKE at smaller scales than those for the Ysu scheme (e.g., Ti2 run) given the same value of . Low turbulence energy associated with simulations using Ysu-PBL scheme indicate that the Myn-PBL scheme is more impacted by the terra incognita than the Ysu-PBL scheme. As expected, the spectra for the *w*-component velocity display lower energy than those for the *u*-component velocity. To determine the effect of three-dimensional turbulence modeling on the model output, the Lilly model was substituted for the Myn/Ysu-PBL scheme in the domain D02 for the above cases. The spectra of the *u*- and *w*-component velocities for these cases [Fig. 6d (*u* spectra) and Fig. 6e (*w* spectra)] show similar behavior for the spectral energy for varying cases as was found for cases that use the standard mesoscale PBL scheme, such as increasing turbulence energy density with decreasing . However, the shape of the curve of spectral energy density and the amount of resolved turbulence are different. For example, cases that use the Lilly model predict more spectral energy up to the effective resolution of WRF. Similarly, the *w*-component velocity contains more resolved spectral energy when the Lilly model is used in D02 compared to when a PBL scheme is used. The PBL scheme designed for the mesoscale simulation can overestimate the SGS motion, producing less energy in the resolved scales of motion. This is the reason why the magnitude of *w* component of velocity reduces when using the PBL scheme compared to the LES closure. Shin and Dudhia (2016) and Zhang et al. (2018) also observed similar results for vertical velocity when they used the PBL scheme. The opposite behavior of PBL scheme applies in Lilly model—with more energy in the resolved scales of motion.

Another way of evaluating the resolved turbulence in the simulated flow field is the use of POD analysis. The POD energy calculated for domain (i.e., D02) with grid spacing in the terra incognita for different cases (i.e., Ti1–Ti5) that use Myn- or Ysu-PBL schemes is presented in Fig. 6c. These results show that, for all cases, the POD energy in the first few modes (10) of domain D02 is excited compared to that of domain D01 due to the decrease of . The case Ti1 shows the largest increment of POD energy in the first mode, possibly because of the smallest grid refinement ratio (i.e., 2), which allows the translation of more flow structures into domain D02. The POD energy starting from the second mode onward does not increase systematically with the decreasing , possibly due to the varying grid refinement ratio, size of , etc. As was found for the results from the spectra, there is some variation in the POD energy between the Myn- and Ysu-PBL schemes. Replacing Myn/Ysu-PBL schemes with the Lilly model in domain D02 produces similar results for the first mode (Fig. 6f) regardless of the turbulence model used (i.e., PBL scheme and Lilly model). This indicates that both turbulence models inherit the mesoscale features originated from domain D01. However, POD energy is highly excited for the mode 2 in cases that use the Lilly model, and the energy is distributed among the large number of spatial modes (15).

### d. Microscale limit

#### 1) Flow visualization

The horizontal grid spacing of the intermediate domain (between mesoscale and microscale domains) may fall into the terra incognita and may produce secondary structures in the boundary layer. To study the effect of the terra incognita on microscale result, we analyzed the model output from the microscale domain (of two horizontal grid spacing 0.24- and 0.04-km cases) from the coupled simulation, where the microscale domain was forced by the parent domain with grid spacing within the terra incognita. Usually, on a microscale domain a three-dimensional model (e.g., Lilly model) is used to represent the effects of unresolved scales of motion on the resolved fields. Herein, in addition to the Lilly model, PBL schemes (i.e., Myn and Ysu) were used for modeling the turbulence in the microscale domain to examine the model output due to two different turbulence modeling approaches. Although PBL parameterization with grid spacing 100 m is not common in the modeling community; the cloud-resolving WRF simulations that uses grid spacing 1 km with standard PBL parameterizations is common for studying many weather events, such as the convective updrafts (Fan et al. 2017) and tornado outbreaks (Fierro et al. 2012).

Vertical and horizontal contour plots of wind speed help evaluate flow structures resulting from two different types of turbulence modeling—PBL schemes and the Lilly model. Figure 7a shows the *x*–*y* contour plots (i.e., horizontal slice) of an instantaneous *w* component of velocity at two instants of time (1500 and 1700 CST) obtained from the microscale domain with 0.24 km using either the Myn-PBL scheme or the Lilly model and forced by parent domain D01 or D02 with varying . This allows us to study the impact of forcing of parent domains within the terra incognita on the microscale domain. Contours for *x*–*y* planes were extracted from the seventh model layer (95 m above the surface). The plots for 1700 CST are an enlarged view of the area shown in the square box on the panels for 1500 CST. The enlarged view is derived from the northern region of domain D02 to minimize the impact of the spinup distance on the results. The *x*–*y* slices in rows 1–6 evaluate the results caused by variation of spinup distance from inflow boundaries, turbulence modeling schemes, and grid refinement ratio. Note that the original domain size of the simulations shown in rows 5–6 at 1500 CST are larger than is shown in the plot. The flow structures presented in the first two rows display different characteristics as they use different approaches for modeling turbulence—Myn-PBL scheme and the Lilly model. The structures in the simulations using the Lilly model show streamwise elongated updraft and downdraft orientating along the mean wind direction. On the other hand, structures in the case using the Myn-PBL scheme show cellular-like structures, and the magnitude of velocity is smaller compared to the result from Lilly model. The differences in velocity magnitude and flow structures between the results using the Lilly model and PBL schemes is associated with the parameterized SGS motion in two turbulence modeling approaches. Compared to the MYN PBL (i.e., local) scheme, the YSU PBL (i.e., nonlocal) scheme overestimates the coherent structures (more mixing) in SGS motion. This overestimate results in decrease of the resolved motion that reduces the magnitude of resolved up/downdrafts in YSU scheme than in Myn scheme. Similar results were observed in the study done by Shin and Hong (2015) and Shin and Dudhia (2016). To study the discrepancies seen in results produced by the cases that used the Myn PBL scheme and Lilly model, we examined the horizontal and vertical wind speed and found that the case that used the Myn PBL scheme produced smaller horizontal wind speed (1.5 m s^{−1} less in the seventh model layer) than the case that used the Lilly model. However, the wind direction and surface flux were found to be similar in both simulations. The larger horizontal wind speed and fluctuations of vertical velocity with relatively moderate heat flux in Lilly model favored the formation of streak-like structures near the lower boundary layer. On the other hand, smaller wind speed and vertical velocity in the case that used the PBL schemes produces the cell-like structures, particularly for cases with smaller grid spacing. At 1500 CST, the contours in rows 1–2 display a larger spinup distance from the inflow boundaries (i.e., from the south and east sides) due to the larger grid refinement ratio (i.e., 16; 4 each between D01 and D02 and D02 and D03) compared to that for contours in rows 3–4. The structures look very similar (allowing for the spinup distance) among the cases that use the same turbulence model, which indicates that the results are insensitive to the number of parent domains and the grid refinement ratio at a sufficient distance from the inflow boundary. Moreover, contour plots in rows 5–6 show results from the domain that is forced by a single parent domain within 7 of the mesoscale limit, whereas contour plots in rows 1–4 are from the domain forced by two nested domains with a 7 of the immediate parent domain within the terra incognita limit. However, the flow structure discerned from the panels in odd or even rows are similar. This suggests that the flow in the microscale domain is insensitive to whether its parent domain is within mesoscale or terra incognita limits.

A microscale simulation with smaller horizontal grid spacing is expected to resolve more turbulence and reveal more details of the flow structure. Contour plots (in *x*–*y* and *y*–*z* planes) of the instantaneous *w* component of velocity for cases with 0.04 km are shown in Fig. 7b. The panels in the left column show *x*–*y* (horizontal slice) contour planes at two instants of time: 1500 and 1700 CST extracted from the seventh model layer (95 m above surface) for different cases. Again, the velocity plots at 1700 CST are an enlarged view of the area shown in the square box for 1500 CST. Contours shown in the first four top rows are from cases that use the Lilly model with different grid refinement ratios. Flow structures seen in the top four rows are similar. For example, the flow field for all cases at time 1700 CST is covered with streamwise elongated updraft and downdraft structures along the mean wind direction. This indicates that there is no direct effect of grid refinement ratio on the flow structures (e.g., cases Ti4 and Ti6). The horizontal slice planes shown in the bottom two rows display different patterns compared to each other or to the top four rows. Although the case R1(Lilly) (in the fifth row) uses the Lilly model, the structures are different from the nested domain case (i.e., from the top four rows). Note that case R1 uses the NARR reanalysis data to provide forcing in the microscale domain. The stability parameter for this case R1 is 25 and the streak-like structure along the mean wind direction occupies more than half of lower height of , indicating the possible formation of horizontal roll. As observed in the case with 0.24 km (in Fig. 7a), the case Ti6(Myn–Myn–Myn) (in the sixth row) shows irregular cellular-like flow structures. Contours in the vertical *x*–*z* plane (Fig. 7b, right column) slice through the middle of the domain of plots in Fig. 7b (left) at 1700 CST. The vertical slices for the first top four rows clearly show similar vertical structures, whereas the contours in the bottom two rows show different flow structures within the boundary layer. The cases shown in the last two rows resolved smaller energy compared to the cases in the top rows and the flow structures are different as well as they seem to be less realistic, suggesting that they are poor microscale simulations. This result highlights the importance of utilizing the nested runs and a three-dimensional turbulence model in microscale simulations.

#### 2) Spinup distance

The velocity contours of Fig. 7b (first column) showed that the turbulence for the southeasterly flow developed after a given distance from the upwind boundary of the domain and that the distance varied with the grid refinement ratio. To evaluate the resolved turbulence, spectra along five locations (i.e., L1–L5) with increasing distance from the inflow boundary (Fig. 8a) were calculated for the case with 0.04 km and using the Myn–Myn–Lilly turbulence parameterization (Fig. 8b). All five simulation cases use identical domain configuration and lateral boundary forcing. The spectra for location L1 shows the least resolved turbulence energy among the spectra for the five locations, as expected, whereas the spectra for the locations L3–L5 show a steady growth in spectral energy density compared to L1. This variation of spectral energy density along the domain is caused by the mismatch between the resolved scale of motion of parent domain and the spatial scale (i.e., ) of the microscale domain, as well as by the transition of energy from the large scale to the small scale. The development of the spectral energy density indicates that a certain distance is required to spin up the turbulence. Adding perturbations at the inflow boundary (Muñoz-Esparza et al. 2014) would help to develop the turbulence over a shorter distance. To provide a quantitative result, the spectral curves for frequencies larger than 0.004 Hz (Fig. 8b, shaded region) of all five cases for each location were integrated to calculate the variance. It should be noted that the five cases represent different combinations of PBL schemes and Lilly model in the nested domains. Figure 8c depicts the variance at each location for five cases with different combinations of turbulence parameterization/modeling approach. Results show that the amount of resolved turbulence in the cases that use the Lilly model in the innermost domain is larger than for cases that use PBL schemes for the innermost domains. For cases that use the Lilly model, the turbulence is nearly resolved after location L3, indicating that almost one-third of the distance (for a 10 km × 18 km domain) is needed to fully develop the turbulent flow field. On the other hand, cases that use the Myn- or Ysu-PBL schemes for the innermost domain results in much lower energy in all five locations as demonstrated by their spectra and POD energy. These variations in resolved turbulence for different model configurations confirm that the amount of resolved turbulence in the simulation depends on both the type of turbulence parameterization/modeling and the distance from the inflow boundary. Given that the turbulence is totally developed after 1/3 of the way into the domain from the south, only the data from upper half of the domain is used in the subsequent analyses.

#### 3) TKE, 1D spectra, and POD

The variance of a variable provides an estimate of the amount of turbulence resolved in a simulated flow field. To calculate the profile of TKE for different cases, the sum of variances using the resolved flow field was computed at each vertical grid point using time series of the three components of velocity. Three TKE profiles were calculated using each 1-h time series data (i.e., 1400–1500 CST, etc.) and then averaged. The top two panels in Fig. 9a show the vertical profile for TKE for the cases with 0.24 km that use Myn/Ysu-PBL schemes and the Lilly model. For the Myn-PBL scheme and 0.24 km case (first row), the TKE value from the surface up to for all cases is found to be similar, but it increases toward the top part of the boundary layer for cases Ti2 and Ti6. In contrast, the Ysu-PBL scheme predicts much smaller TKE in the entire vertical profile, indicating that the Ysu-PBL scheme resolves less turbulence than the Myn-PBL scheme. For the real case simulation, Zhang et al. (2018) also found a more intense *w* component of velocity when using Myn scheme than the counterpart Ysu PBL scheme. Similarly, Shin and Dudhia (2016) reported that a strong up/downdraft in the convective condition when using the Myn-PBL scheme. The variation of magnitude in the TKE profile in the same PBL scheme cases may result from different factors including grid refinement ratios. The results for the case with the Lilly model and with 0.24 km shows that the TKE increases throughout the vertical profile, but not systematically, as the changes. The TKE value for all cases tends to converge near an altitude of 0.25. Below this point, the TKE value increases significantly and the shape of the profile is quite different from those of the PBL scheme cases. The cases with horizontal grid spacing of 0.04 km using PBL schemes (Fig. 9d, top panel) show similar but slightly larger TKE profiles compared to the 0.24-km case. The use of the Lilly model (instead of the Myn/Ysu-PBL scheme, Fig. 9d, bottom panel) leads to an increase in the amount of turbulence resolved as in these cases with 0.24 km. The TKE profile above 0.5 shows more convergence for the cases with smaller horizontal grid spacing. The case Ti4(Myn–Lilly–Lilly) predicts the largest TKE throughout the vertical profile, whereas the single domain case R1(Lilly) produced the least TKE for the entire vertical profile, again highlighting a potential shortcoming in the single domain case R1 and importance of coupled simulation.

The spectral representation estimates the turbulence resolved of the simulated flow. Spectra of *u* and *w* component of velocity for the cases with 0.24 km were computed from time series data sampled at 0.1 Hz at 95 m above surface (Fig. 9b). The parent domain D02 for cases Ti1–Ti3 has in the terra incognita, whereas the parent domain D01 for case Ti6 has = , which is in mesoscale limit. The result for the *u*-component velocity in the top panel (that uses the Myn-PBL scheme) shows that, except in the case Ti2(Ysu–Ysu–Ysu), all cases predict similar spectral energy for frequencies greater than 10^{−3} Hz. The discrepancy in energy for the higher frequencies between two cases (i.e., Ti2 and Ti6) that use PBL-Ysu schemes may be due to the number of grids used to downscale the spatial scale to 0.24 km. Not shown here, spectra for *v* and *w* components of velocity also predicted similar spectral behavior as predicted by the *u* component of velocity. Keeping the same domain configurations, but replacing the PBL turbulence schemes Myn- and Ysu-PBL with the Lilly model (Fig. 9b, bottom panel) shows that the spectral energy of most cases are similar, and that there is no underresolved spectral energy associated with the use of the Ysu-PBL scheme (see Fig. 9, top panel). Compared to cases with Myn- and Ysu-PBL schemes, these cases display a larger magnitude of the spectral energy and resolved turbulence up to the effective resolution of WRF (i.e., 7). The case Ti1(Myn–Lilly–Lilly) produced slightly smaller values of spectral energy. It should be noted that, for both PBL schemes and the Lilly model, the child domain in case Ti6 receives forcing from a parent domain with km and in cases Ti1–Ti3 receives forcing from a parent domain with grid spacing in the terra incognita [with (i.e., from 0.48 to 1.44 km)]. This indicates that, for both the Lilly model and the Myn PBL scheme, the horizontal grid spacing of the parent domain, either in the terra incognita or the mesoscale limit, does not affect the results of the child domain. Figure 9e shows spectra for the cases with 0.04 km that use Myn/Ysu-PBL schemes and Lilly model. Because of the fine and representation of smaller eddies in cases with 0.04 km, more turbulence scales are captured than in cases with 0.24 km. These cases exhibit similar spectral energy behavior as the cases with 0.24 m, such as having larger magnitudes of spectral energy for the Lilly model. The case Ti6(Ysu–Ysu–Ysu) shows smaller spectral energy at the higher frequencies compared to other cases that use PBL schemes. In addition to the multiple nested domain cases, the single domain simulation R1—that uses the NARR reanalysis data for its lateral boundary forcing—also uses the Lilly model. Results from the case R1 show a little less spectral energy at the longer wavelengths, but has comparable spectra for shorter wavelength.

To estimate the amount of fluctuation of variables in the simulated flow using three velocity components at a number of locations, the POD energy of the first 21 spatial modes were calculated for the case with km that uses the Myn/Ysu-PBL scheme (Fig. 9c, top panel). The POD energy in the first mode is found to be similar to in the finding from Figs. 4c and 6c, suggesting that the mesoscale features have been translated to the nested domain. From the second mode onward, as expected, the POD energy increases slightly as decreases. However, the POD energy for the cases that use Ysu-PBL scheme is smaller compared to the cases that use the Myn-PBL scheme. There are variations in the amount of energy resolved in the various cases that use the Myn-PBL scheme. The discrepancy in POD energy for modes can be attributed to the type of PBL parameterization used and the grid refinement ratio. The cases that use the Lilly model (Fig. 9c, bottom panel) show an increase of POD energy in all spatial modes for all cases. Again, the POD energy of the first mode for all cases is similar to that found using the Myn-PBL schemes. Furthermore, there is no noticeable difference in POD energy due to turbulence parameterization/modeling is used to force the parent domain or to the size of the grid refinement ratio used between the domains. Recall that the parent domain for cases Ti1–Ti3 is within the terra incognita, whereas the parent domain for the case Ti6 falls into mesoscale limit. Also note that the parent domain D02 in cases T2(Myn–Lilly–Lilly) and T2(Myn–Myn–Lilly) used the Lilly turbulence model and the Myn scheme, respectively, but their nested domain that used Lilly turbulence model produced similar results. The POD energy for the cases with 0.04 km that uses the Myn/Ysu-PBL scheme and the Lilly model is shown in Figs. 9f. The trend of the POD energy is similar to the cases with that for 0.24 km. The magnitude of POD energy for the first mode is similar to that of cases 0.24 km, but the POD energy in the higher spatial modes is highly excited when using both turbulence models (i.e., PBL schemes and Lilly model). However, the cases that use the Lilly model show a larger amount of resolved turbulence resolved in the POD spatial modes .

The importance of using a coupled nested domain simulation rather than a single domain simulation for 0.04 km can be evaluated using the results (i.e., flow structures, spectra, and POD) of both coupled and single domain simulations discussed above. Note that the domain size and horizontal grid spacing of the innermost domain of the coupled simulation cases (i.e., Ti4–Ti7) are identical to that of case R1 without coupling. Both the coupled and single domain simulations use the NARR reanalysis data as the forcing. The only difference is that the innermost domain in the coupled simulation receives the lateral forcing from the flow field of a parent domain that is dynamically evolving, whereas the single domain case uses interpolated lateral forcing from much coarser data (i.e., NARR with grid spacing of 32 km) than its domain size. The POD energy for the single domain case R1 using the Lilly model (Fig. 9f, bottom panel) predicts the least energy for all POD spatial modes. Similarly, the *x*–*y* contour plot of wind speed (Fig. 7a, fifth row) reveals that the single domain case R1 captures the wind direction but the simulated flow structures are much different than for the coupled simulation cases. The spectra at 95 m above the surface shows that, except at the lowest frequency, case R1 is able to resolve the turbulence similar to the coupled simulation cases. Thus, the results from examining POD energy, *x*–*y* contours of wind speed, and spectra at lowest frequencies indicate that the single domain case R1 lacks both resolved turbulence and mesoscale features. The smaller amounts of energy in the first few modes of POD analysis are associated with the grid resolution of the reanalysis data, the domain size of microscale domain, and the grid refinement rations (from 32 km to 40 m). The size of the microscale domain used in case R1 is small compared to the horizontal resolution of the reanalysis, with only one vertical profile of reanalysis data. In addition, the larger refinement ratio used for this case requires a spinup distance larger than the small microscale domain. Use of large microscale domain having hundreds of grid points could improve results as it contains more reanalysis data and provides large spinup distance, but it would be very computationally expensive. This indicates that there is not much of a mesoscale footprint present in the R1 case. The flow structure and POD energies in the R1 case in the higher POD modes might be improved by adding perturbations near the boundary or using periodic boundary conditions. However, the POD energy of the first mode cannot be augmented unless the microscale domain is forced by the flow field evolved dynamically in the mesoscale domain. This indicates the importance of coupling between mesoscale and microscale domains for simulating highly resolved atmospheric flow.

#### 4) 2D spectra

Flow structures in the microscale domain vary with the type of turbulence modeling schemes used in the parent domain (refer to Fig. 7). Between the flow structures simulated using the PBL schemes and the Lilly model, it is difficult to say which scheme or model produces more realistic flow structures. As discussed earlier for Lilly model case, the stability parameter 40 and the occupancy of streak-like structures in the first 15% of do not support that the structures are horizontal rolls. In this context, two-dimensional spectra of velocity can reveal the expected flow structures in the microscale simulation. Figure 10 shows the two-dimensional spectra derived from the *u* component of velocity for the horizontal plane at 95 m above the surface. These spectra display the spectral energy in terms of wavenumbers and in the *x* and *y* direction, respectively. Data from the subregion (21.6 km × 21.6 km) of the northwestern portion of domain D03 were used for computing spectra to minimize the effect of spinup distance on the simulated flow field (see Fig. 8). The data can be considered statistically horizontally homogeneous since the plane of data is nearly horizontal and the boundary layer turbulence is fully developed for the period considered for this analysis. Two-dimensional spectra were calculated using plane of data saved every 3 s for a period of half an hour and then time averaged for smoothing. Figures 10a–h display spectra for the cases with 0.24 km but with different combinations of turbulence modeling (i.e., Myn/Ysu-PBL schemes and Lilly model). The contour lines indicate changes of an order of magnitude in spectral energy. The spectra in the top row show the cases that use three nested domains to produce a flow field with 0.24 km in domain D03. Figure 10a displays the spectral energy for the case Myn–Myn–Myn, where all domains use Myn-PBL scheme. The results show that the spectral energy, except in low wavenumber space, decreases in a similar manner in both and . In contrast, the cases that use the Lilly model in domain D03 (Figs. 10b and 10d) indicate that the contours of spectral energy are more squeezed along the mean wind direction than across the mean wind direction. The simulation using Ysu–Ysu–Ysu (Fig. 10c) shows the spectral energy contours that are a little less elongated across the mean wind direction, but the contours are denser at the smaller wavenumbers (larger sales). Between the two turbulence model cases (i.e., Myn-PBL scheme and Lilly model in domain D03), the case with the Lilly model is able to capture more spectral energy in the higher wavenumbers compared to the case with the Myn-PBL scheme. This decay of turbulence and the shape of the spectral contours seen in the Lilly model case are in accordance with the turbulence theory—the ratio of spectra of longitudinal to transverse directions is greater than that in the inertial subrange (Tennekes and Lumley 1972). The compaction of spectral contours in two-dimensional spectra were also mentioned by Gibbs and Fedorovich (2014) in their comparison of spectra from the WRF and in-house LES generated data. Based on these two-dimensional spectra and results presented earlier, the Myn- or Ysu-PBL scheme should not be used for microscale simulation.

The two-dimensional spectra can also be used to analyze the impact of using a forcing domain within the terra incognita on the microscale output. The number of nested domains and the grid refinement ratios are different between the cases shown in Figs. 10a–d and 10e–h, although their inner domains have the same grid spacing 0.24 km). Figures 10e–h show spectra for the cases that have two nested domains with a grid refinement ratio of 10. The two-dimensional spectra show that, except for the Ysu–Ysu–Ysu case, all other cases (i.e., Figs. 10e,f,h) produced similar spectra to those shown in Figs. 10a–d. The case Ysu–Ysu (Fig. 10g) with a grid refinement ratio of 10 is able to capture spectral energy in larger wavenumbers. In the above cases, the two nested domain cases use , whereas the three nested domain cases use in the parent domain to force the domain with 0.24 km. However, the results from two nested and three nested domains are very similar. This suggests that the size of the grid spacing in the forcing domain does not play an important role in resolving the spectral energy. However, the inflow distance to develop turbulence in the nested domain is affected by the size of grid spacing of the parent domain. The two-dimensional spectra for the cases with 0.04 km show similar behavior to the cases with 0.24 km in wavenumber space (Figs. 10i–l), such as a compaction of the spectral contours along the mean wind direction. The subregion of (4.8 km × 4.8 km) from the northern part of the domain was used for extracting the data. This subregion from the northwestern of the domain was selected as it provides fully developed flow for the analysis. Because of the smaller grid spacing, more spectral energy is captured for smaller wavenumber space in these cases. The different nature of contour plots in two-dimensional spectra with respect to the turbulence modeling approach in both cases (i.e., with 0.04 and 0.24 km) can also be explained from the instantaneous *w*-component velocity in the horizontal plane as shown in Fig. 7. The cellular-like irregular structures seen in the case with Myn-PBL scheme lacks directionality; therefore, their spectra are symmetric and decay in a similar manner along and wavenumber space. In contrast, elongated and streaky structures oriented along the mean wind direction in the case with the Lilly model are responsible for producing elongated contours of spectra across mean wind direction. Moreover, the shape and orientation of structures among the cases discussed previously imply that the results from the case that use the Lilly model with 0.24 km (among those shown in Fig. 7) may better represent the flow features than cases that use the Myn-PBL scheme given the same .

## 5. Summary and conclusions

A number of coupled mesoscale–microscale simulations with varying horizontal grid spacing falling into three limits—mesoscale, terra incognita, and microscale—were performed using the WRF modeling framework for a case with a convective atmospheric boundary layer. Combinations of two PBL parameterizations and a three-dimensional turbulence model were used for modeling turbulence on a nested domain. High-frequency model output of 3-h window split into three records was used to calculate the spectral energy density, POD energy, and TKE.

Our results grouped into different boundary layer heights showed spurious secondary structures for the horizontal grid spacing smaller than , indicating that a grid resolution comparable to represents the boundary between mesoscale simulations and the terra incognita—and that horizontal grid spacing smaller than fall into the terra incognita. This behavior was more evident when using the Myn-PBL than the Ysu-PBL scheme. To examine the effect of grid spacing using the terra incognita, tests were completed in which a microscale domain was forced by a parent domain with several horizontal grid spacing and with different approaches for turbulence parameterization/modeling. The spectra, POD energy, and TKE showed that the results in the microscale domain do not vary with the type of turbulence parameterization/model and the size of the grid spacing (i.e., terra incognita or mesoscale limit) used in the parent domain. Rather the results depend on the type of turbulence model used in the microscale domain itself. In practice, the cloud-resolving WRF model uses a grid spacing 1 km with the PBL parameterization for modeling the turbulence. However, for a grid spacing of m, our results suggest that care is needed when choosing the turbulence model. The effect of grid refinement ratio is observed only in the spinup distance of turbulence from the inflow boundary. In addition to spinup distance, Mazzaro et al. (2017) found that the grid refinement ratio has impact on many turbulent statistics, such as wind speed, TKE, and stress in the nested simulation.

For a given horizontal grid spacing and turbulence model in the microscale domain, the results using the Myn scheme showed cellular-like structures in the wind speed 95 m above the surface, whereas results using the Lilly turbulence model showed elongated and streak-like structures along the mean wind direction. For the Lilly model case, the stability parameter () was found to be 40 and the streak-like structures occupy the first 15% height of . This indicates that the structures are less likely to be horizontal rolls, rather they are fast and slow streaks of wind. The ratio of wavenumber along and across the streamwise flow in this southeasterly flow case favored the flow structures predicted by the Lilly turbulence model case. This implies that for a microscale domain, simulations using the Myn- or Ysu-PBL schemes would not produce realistic flow structures. Allowing adequate space for turbulence spinup, the results indicate that the grid refinement ratio between the parent and nested domains did not impact the results, provided that a similar turbulence model (either PBL scheme or Lilly model) is used in the nested microscale domain.

The POD mode energy calculated from the vertical profile of three components of velocity for the microscale domain showed very similar amounts of energy in the first two spatial modes, irrespective to the types of turbulence modeling/schemes and the different domains. However, the POD energy at higher spatial modes was larger for the Lilly turbulence model compared to that from the Myn-PBL scheme case. The similarity of POD energy from the first two modes is associated with mesoscale features that were translated from the mesoscale domain to the microscale domain. However, the POD energy of the first few modes for the single domain simulation driven with NARR data showed lower energy compared to that from the nested domain simulation for a given grid spacing. These results clearly demonstrate the importance of coupling simulations in order to fully capture the mesoscale features in the simulated flow field. The behavior of model performance in the terra incognita for the complex terrain would be different than for the flat terrain, which needs to be investigated in the future.

## Acknowledgments

This research was supported by the Office of Energy Efficiency and Renewable Energy of the U.S. Department of Energy as part of the Wind Energy Technologies Office. The Pacific Northwest National Laboratory is operated for the DOE by Battelle Memorial Institute under Contract DE-AC05-76RLO1830. JDM’s contribution was prepared by LLNL under Contract DE-AC52-07NA27344.

## REFERENCES

*Random Data: Analysis and Measurement Procedures.*Wiley Series in Probability and Statistics, Vol. 729, John Wiley & Sons, 640 pp.

*Proc. IBM Scientific Computing Symp. on Environmental Sciences*, Yorktown Heights, NY, IBM.

*Atmospheric Turbulence and Radio Wave Propagation*, A. M. Yaglom and V. I. Tatarski, Eds., Nauka, 166–178.

*A First Course in Turbulence.*The MIT Press, 300 pp.

## Footnotes

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