Efficacy of the Cell Perturbation Method in Large-Eddy Simulations of Boundary Layer Flow over Complex Terrain

Abstract

A challenge to simulating turbulent flow in multiscale atmospheric applications is the efficient generation of resolved turbulence motions over an area of interest. One approach is to apply small perturbations to flow variables near the inflow planes of turbulence-resolving simulation domains nested within larger mesoscale domains. While this approach has been examined in numerous idealized and simple terrain cases, its efficacy in complex terrain environments has not yet been fully explored. Here, we examine the benefits of the stochastic cell perturbation method (CPM) over real complex terrain using data from the 2017 Perdigão field campaign, conducted in an approximately 2-km wide valley situated between two nearly parallel ridges. Following a typical configuration for multiscale simulation using nested domains within the Weather Research and Forecasting (WRF) model to downscale from the mesoscale to a large-eddy simulation (LES), we apply the CPM on a domain with horizontal grid spacing of 150 m. At this resolution, spurious coherent structures are often observed under unstable atmospheric conditions with moderate mean wind speeds. Results from such an intermediate resolution grid are often nested down for finer, more detailed LES, where these spurious structures adversely affect the development of turbulence on the subsequent finer grid nest. We therefore examine the impacts of the CPM on the representation of turbulence within the nested LES domain under moderate mean flow conditions in three different stability regimes: weakly convective, strongly convective, and weakly stable. In addition, two different resolutions of the underlying terrain are used to explore the role of the complex topography itself in generating turbulent structures. We demonstrate that the CPM improves the representation of turbulence within the LES domain, relative to the use of high-resolution complex terrain alone. During the convective conditions, the CPM improves the rate at which smaller-scales of turbulence form, while also accelerating the attenuation of the spurious numerically generated roll structures near the inflow boundary. During stable conditions, the coarse mesh spacing of the intermediate LES domain used herein was insufficient to maintain resolved turbulence using CPM as the flow develops downstream, highlighting the need for yet higher resolution under even weakly stable conditions, and the importance of accurate representation of flow on intermediate LES grids.

1. Introduction

Over recent decades, computational power has grown such that it is now becoming routine to perform multiscale atmospheric simulations, in which a larger-scale mesoscale flow field is dynamically downscaled using grid nesting, with finer mesh spacing over regions of interest. Such nesting can be pursued successively to large-eddy simulation (LES) scales, within which the energetically important scales of boundary-layer turbulence are explicitly resolved. Near the inflow boundaries of nested large-eddy permitting domains, however, finer-scale turbulence is lacking because it is not explicitly captured in the coarser bounding domain simulation. The development of turbulence in the nested simulation often requires a long distance from the lateral boundaries, referred to as fetch. Goodfriend et al. [1], in a computational fluid dynamics code, and Mirocha et al. [2], using a numerical weather prediction (NWP) model, explored the transition distance required for turbulence to be generated at a coarse-to-fine grid interface. Both studies found that it takes significant fetch to sufficiently develop turbulence compared to stand-alone periodic simulations at the finer resolution, and that this fetch depends on the chosen turbulence closure model. One way to address this issue is to require very large domain sizes so that there is sufficient distance from the lateral boundary to the region of interest in the flow. This of course results in significantly increased computational costs, considering that up to 75% of the domain can be lost to spurious flow structures according to Mirocha et al. [2].

A recent approach to address this issue is the cell perturbation method (CPM) of Muñoz-Esparza et al. [3], Muñoz-Esparza et al. [4], Muñoz-Esparza and Kosović [5]. The CPM adds small amplitude, $\mathcal{O}$(1 K), perturbations to the temperature field along inflow boundaries to trigger the rapid development of turbulence on the finer nested grid. The CPM has been tested primarily with idealized domains and over flat terrain and found to be successful in greatly reducing the fetch required to develop finer structures, which can develop in as few as 25 points farther into the domain from the perturbation region [3]. This reduces computational cost by reducing the required domain size to develop turbulence on the finer, more expensive grids. Furthermore, the CPM is easy to implement and does not require additional cost due to larger domains or from running a priori simulations as in other turbulence generating methods [6].

An outstanding question for multiscale simulation over complex terrain is whether the topography itself suffices to generate finer scale turbulence in the nested domains, or whether additional turbulence generating techniques are effective. In some applications over complex terrain, a main improvement with nesting to high resolutions is that topography-driven flows such as gap flows can be resolved [7,8]. In many cases, however, it has been found that increasing the grid resolution by nesting does not lead to improved results [9]. Instead, the large flow features from the outer domain carry over into the inner domain, and detailed examination of the flow shows lack of the expected finer scale flow structures on the inner domain. This can be due to insufficient detail in the forcing provided by topography and land cover data [10], but also an insufficient fetch for finer-scale turbulence to develop may be a likely cause. A typical pattern observed is “streaky” roll structures in the wind field [3,11,12] or spurious large-scale convection cells which persist on the inner domain [13,14]. Mirocha et al. [2], Ching et al. [14] found that prohibitively large fetches were required for turbulence to develop in an LES nested in a mesoscale simulation performed by the Weather Research and Forecasting (WRF) model. The development of turbulence was aided by the application of convective forcing or by adding idealized, sinusoidal topography. Nonetheless, 100 s or 1000 s of points, depending on the details of the turbulence parameterization, were required for proper turbulent scales to develop. As seen by Mirocha et al. [2], Muñoz-Esparza et al. [3,4], in some cases it can take up to three-quarters of the inner domain length before finer structures develop. If the inner nest is too small, then it is possible for almost no finer structures to develop, and the benefit of the finer resolution may largely be lost. Warner et al. [15] more than 20 years ago discussed the need to place nested domain boundaries far from the region of interest to minimize spurious influence from the grid nesting boundaries. Now that LES nested within NWP is becoming more commonplace, a new focus on the representation of turbulent motions on the inner grids is warranted.

The detailed turbulent motions at the microscale are of particular interest in applications such as wind energy. Prediction of power generation, blade fatigue, and turbine wake interactions depends critically on the representation of turbulence [16]. Wind energy forecasting over complex terrain was a primary focus of the Perdigão field campaign which took place in Portugal in 2017 [17]. Nested large-eddy simulations performed for the Perdigão experiment also develop unphysical streaky structures despite the complex topography present, as shown later below. At the high computational costs required for LES, it is paramount to develop realistic turbulence as fast as possible in the inner domains.

This work focuses on the benefits of CPM to reduce the large fetch required to generate finer turbulent motions on nested domains over complex terrain. The Perdigão field site is used as a test case. CPM has previously been tested in real case simulations of two wind farms over relatively flat terrain by Muñoz-Esparza et al. [18], Arthur et al. [19]. Based on spectra of vertical velocity, Muñoz-Esparza et al. [18] showed unperturbed LES, $\mathsf{\Delta }=$ 90 m, took 45 km of fetch to sufficiently develop turbulence, while applying perturbations reduced this to 9 km. Arthur et al. [19] found marginal benefits on mean flow fields from using CPM; however the structure of the turbulence was certainly enhanced with much finer scales as seen in velocity contours. Arthur et al. [19] applied CPM between domains with 250 m grid spacing and 50 m grid spacing, and again between the 50 m grid and a final nest with 10 m grid spacing. However, this fine LES was intentionally designed with 300 grid points of fetch before the model performance versus field observations was evaluated, making the additional effects of CPM difficult to extract. Their analysis focused on near-neutral and slightly unstable conditions at a field site with terrain smoothly varying by roughly 40 m in elevation over 5 km.

Here, we focus on the Perdigão field site which consists of two approximately parallel ridges and complex hilly terrain. There is a single wind turbine on one of the ridges, and one of the overall goals of the project is to provide data to improve simulation capabilities, including the generation of accurate turbulent structures as well as the impacts of terrain and atmospheric stability on wake behavior in complex terrain [17]. In this paper, we prioritize testing a grid nesting framework to consistently and efficiently generate appropriate turbulent motions on the finer grids. The CPM is used to trigger turbulence on a 150-m resolution LES domain, nested within a 2250-m resolution mesoscale domain. Other work is pursuing finer nests to 50 m and 10 m resolutions [20,21], but our focus here is on the 150 m domain where spurious numerically generated rolls are apparent and the model struggles to develop finer structures.

The goal of this work is to study the effect of the CPM on the resolved turbulence when combined with the effects of complex topography. The previous real case studies that have performed LES with CPM have used fine meshes and large domains, which require significant computational costs [18,19]. Other studies have considered turbulence in ‘coarse’ LES over complex terrain but without CPM [12,22]. We will examine the relative impacts of topography and the CPM in multiscale WRF simulations using observations from Perdigão. Both mean flow and turbulence characteristics will be examined using the CPM relative to two different resolutions of the underlying terrain to examine impacts of terrain steepness and smaller-scale variabilty. The model results are compared to field observations from the Perdigão field campaign (Section 4.1). We then analyze the performance of the CPM by comparing different domain sizes (fetches) and turbulence energy spectra (Section 4.3, Section 4.4 and Section 4.5). We discuss the computational savings made possible by the CPM, as well as the sensitivity of these results to the time of day and to the resolution of model topography. The CPM is successful at reducing fetch and is especially cost-saving during the weakly convective period during which unperturbed simulations require the longest fetch to fully develop turbulence.

2. Background

Several methods have been proposed to address challenges faced by grid nesting and the related generation of appropriate turbulent length scales. These range from turbulence recycling methods [23], which artificially increase the fetch by ‘recycling’ information on the turbulence within the domain to its inflow boundary, to the purely stochastic methods, which instigate turbulence through perturbations not based on fundamental physical equations [24]. Other techniques, the prepared methods, use a priori simulations with periodic boundary conditions to then augment the inflow of the nested simulations so that physically generated turbulence is present at the boundary [6]. However, performing a priori simulation also has a computational cost that could negate the gains of the turbulence generation method. Using engineering fluid mechanics codes, the wind energy community has long employed many of these methods for synthetic turbulence generation [25]. Due to the complexity of meteorological models, not all techniques are readily portable to NWP codes. For example, the existence of multiple inflow boundaries makes it unclear how to generalize recycling methods. Furthermore, the a priori simulations used in prepared methods generally require periodic boundary conditions. In addition to issues of induced periodicity, periodic boundary conditions could not be used for nested simulations which must be forced by their parent at the boundary.

Given these issues with recycling and prepared methods, Mirocha et al. [11] provide an initial formulation of a potential turbulence generation method suitable for NWP applications by applying stochastic perturbations to all the horizontal velocity and potential temperature fields. Muñoz-Esparza et al. [3] improved this initial perturbation technique by applying stochastically generated perturbations to only the potential temperature field over a series of small cells, with discontinuities at the cell boundaries, and generalized the approach via a series of sensitivity studies over a range of wind speeds [4], and later stability regimes [5,18]. Accelerating the development of turbulence in this fashion has shown great efficacy in numerous applications.

Muñoz-Esparza et al. [3] compare the CPM to other stochastic perturbation methods, some requiring a priori simulation, by perturbing only the potential temperature field. Because synthetic turbulence based on momentum perturbation can violate the governing equations, perturbing only the potential temperature is the preferred practice. They find CPM develops realistic turbulence over comparable or fewer points of fetch than the other methods, while retaining a level of simplicity for ease of implementation.

Muñoz-Esparza et al. [4] extend the CPM approach by proposing that the magnitude of potential temperature perturbations for neutral flows be based upon the perturbation Eckert number and the perturbation time scale. This formulation for perturbation magnitude was later modified for stable conditions by Muñoz-Esparza and Kosović [5], suggesting a Richardson number approach instead that improved performance during stable conditions. Here, we do not include the stability correction and instead use the formulation of Muñoz-Esparza et al. [4] to investigate whether it suffices for the weak stability conditions that are the focus of this paper.

3. Methods

3.1. Cell Perturbation Method

Following Muñoz-Esparza et al. [4], pseudorandom perturbations are applied to the potential temperature field near the inflow boundaries over “cells” of 8 × 8 grid points in the horizontal plane. Three such cells, for a total width of 24 grid points parallel to the inflow boundaries, are applied at each vertical level from the first vertical level up to a height of 0.9$z_i$, where $z_i$ is the planetary boundary layer height provided by the mesoscale parent simulation. The maximum height at which perturbations are applied is increased from $\frac{2}{3}{z}_i$ used by Muñoz-Esparza et al. [3,4] in case of shallow boundary layers.

The perturbation temperature, ${\tilde{\theta }}_{\mathrm{p}}$, is drawn from a uniform distribution in the range $[-{\tilde{\theta }}_{\mathrm{pm}},{\tilde{\theta }}_{\mathrm{pm}}]$. Figure 1 illustrates these perturbations for an example domain. We follow Muñoz-Esparza et al. [4] who suggest that the range of the uniform distribution be determined from a perturbation Eckert number:

$$Ec=\frac{U_g^2}{c_p{\tilde{\theta }}_{\mathrm{pm}}}$$

(1)

where $U_g$ is the geostrophic wind, computed as the average wind speed at 1.1$z_i$, and $c_p$ = 1004.6 J kg${}^{-1}$ K${}^{-1}$ is the specific heat capacity of air taken as a constant. Muñoz-Esparza et al. [4] determined an optimal Eckert number, $Ec=0.2$, such that

$${\tilde{\theta }}_{\mathrm{pm}}=\frac{U_g^2}{0.2c_p}.$$

(2)

Another non-dimensional parameter is used to determine the frequency with which new perturbations should be applied to the boundary, the perturbation time scale:

$$\Gamma =\frac{t_p{U}_w}{d_c}$$

(3)

where $t_p$ is the perturbation time period, $U_w$ is representative of the weakest winds at the inflow boundary, and $d_c$ represents a length scale over which the perturbation should be advected before new perturbations are applied. Muñoz-Esparza et al. [4] suggest values of $\Gamma \approx 1$ so that perturbations are not either missing for long periods of time, in the case of $\Gamma \gg 1$, nor rapidly superposed, in the case of $\Gamma \ll 1$. In this study, the CPM uses a value of $\Gamma =0.75$ favoring a slightly more rapid refresh of the perturbations. Additionally, $U_w$ is calculated as the average inflow wind speed at the second vertical level, moved up from the first vertical level in Muñoz-Esparza et al. [4] in case of low near surface wind speeds over complex terrain. The length scale is computed as

$$d_c=\frac{24\mathsf{\Delta }x}{cos\widehat{\varphi }}$$

(4)

where $\widehat{\varphi }$ is the average angle of the wind relative to the inward normal of the inflow boundary. As such, $d_c$ is equal to the width of the perturbation area when flow is normal to the inflow boundary, but $d_c$ becomes larger for skewed angles of approach. Note, a boundary is considered an inflow boundary if it has net inflow at 1.1$z_i$, when averaging across the entire boundary. That is, perturbations are not applied at every instance where inflow occurs over the whole domain but are only applied along the full length of inflow boundaries. For example, if the primary flow direction is from the south-west, then the southern and western boundaries will include perturbations along the entire length of the boundaries (Figure 1).

3.2. Numerical Simulation Setup

The current simulations use the Weather Research and Forecasting (WRF) model version 3.9.1.1 [26] developed by the National Center for Atmospheric Research (NCAR). The WRF model allows multiscale modeling through its nesting capabilities, so that we can generate various microscale LES configurations nested in the same mesoscale parent models. The first two domains, d01 and d02, are mesoscale simulations performed at 6.75 km and 2.25 km, respectively. The microscale simulations, with the prefix d03 used for various configurations, are run at a resolution of 150 m and use an LES turbulence closure [27,28]. A relatively large grid nest ratio of 15:1 (2250 m:150 m), is used between d02 and d03 because a smaller nest ratio would result in grid spacing within the turbulent gray zone, in which it is unclear whether a PBL scheme or LES turbulence closure should be used [28]. It should be noted that a more traditional nest ratio of 3:1 also produced streaky structures in previous work from Mirocha et al. [11] and Zhou et al. [13]. Further details for each domain are given in

Table 1. Information on nested simulation domains include simulation name, horizontal grid spacing, number of grid points, time step, turbulence closure method, resolution of topography, and use of CPM. All domains use the same vertical grid spacing, stretching from 21 m near the bottom surface to 170 m at 2 km a.s.l., and up to 400 m aloft, selected to closely follow the levels used by the ECMWF forcing data below 5000 Pa.

Domain	$\Delta \mathit{x},\Delta \mathit{y}$	Nx·Ny·Nz	dt [s]	Closure	Topography Resolution	CPM
d01	6.75 km	141·141·89	30	MYNN	2 arcmin	off
d02	2.25 km	181·181·89	10	MYNN	30 arcsec	off
d03_30s	150 m	241·241·89	0.5	TKE 1.5	30 arcsec	off
d03_30s_cpm	150 m	241·241·89	0.5	TKE 1.5	30 arcsec	on
d03_30s_ref	150 m	481·481·89	0.5	TKE 1.5	30 arcsec	off
d03_3s	150 m	241·241·89	0.5	TKE 1.5	3 arcsec	off
d03_3s_cpm	150 m	241·241·89	0.5	TKE 1.5	3 arcsec	on

, and Figure 2 illustrates the nest locations. The CPM is applied only at the boundary of the microscale simulations, where they receive inflow information from d02, the parent domain. These simulations are labeled with the suffix _cpm.

The dynamic core of the WRF model is based on a third order Runge–Kutta time advancement scheme, with split time stepping to handle acoustic modes, as well as fifth-order horizontal advection and third-order vertical advection. Physical parameterizations employed in the current simulations include the Noah land surface model [29], the Rapid Radiative Transfer Model for longwave radiation [30], the Dudhia shortwave radiation model [31], and the WRF Single-Moment 5-class microphysics scheme [32]. All domains use a Mellor-Yamada level-3 surface layer parameterization. The mesoscale simulations, d01 and d02, also use the PBL scheme designed to be coupled with this surface layer model [33], while d03 uses a TKE 1.5 LES turbulence closure [27].

Of the land cover data provided by NCAR for use with WRF, the finest available resolution is 30 arcsec, about 1 km at middle latitudes, from the United States Geological Survey (USGS). Our simulations also incorporate a higher resolution, 100 m, land cover data set, the Coordination of Information on the Environment (CORINE), which is then transformed into USGS land use types [34,35].

Topography data is made available at resolutions no higher than 30 arcsec with the standard WRF download. This Global 30 Arc-Second Elevation (GTOPO30) product is used for d03_30s, d03_30s_cpm, and d03_30s_ref. As the LES domain is more finely resolved than GTOPO30 topography, finer resolved topography data can be ingested. Two other LES cases, d03_3s and d03_3s_cpm, ingest 3 arcsec (≈90 m) topography from the Shuttle Radar Topography Mission (SRTM) [36], to test the sensitivity of the CPM to the input topography resolution used only over the smallest domain. Figure 3 shows this 3 arcsec resolution topography next to the 30 arcsec resolution topography for the smaller d03 LES domains.

For all domains, the simulations are initialized with data from the European Centre for Medium-range Weather Forecast (ECMWF) model [37]. The ECMWF High-Resolution Forecast (HRES) model with ≈9 km horizontal resolution and 137 vertical model levels is used to incorporate the highest vertical resolution forcing data available. This is not a standard option in WRF and required modifying the WRF preprocessing system [38]. The model is initialized at 0600 UTC on 20 May 2017. Time is given in UTC for all the analysis, which corresponds with local standard time at the field site (Portugal observes Western European Summer Time such that UTC+1 is the local time on the simulation date). The first 5 h of the simulation are relegated to model spinup and excluded from analysis. The ECMWF data provide the lateral boundary conditions for the coarsest grid at relatively fine, 1 h, intervals throughout the duration of the simulation. For the nested domains, lateral boundary conditions are provided through one-way nesting.

A common strategy to ensure fully developed turbulence in the inner nested simulation domain is to increase the domain size or to locate the area of interest far away from the inflow boundary, i.e., not have it centered within the domain [12,18,39,40]. Both of these strategies serve to increase the fetch or distance from the region of interest to the grid boundary. We therefore also test a larger domain size and compare this to results from smaller domains with and without CPM. This reference LES domain, d03_30s_ref, is twice as large in each horizontal dimension compared to the other LES domains. For the reference simulation, the GTOPO30 topography is used because it is readily available with the WRF download and covers the entire globe. The 30 arcsec topography is likely to be used for the purposes of generating turbulence in a coarse LES nest while high resolution topography is usually only used for finer nests covering less extensive areas. The larger domain, d03_30s_ref, is positioned such that the smaller LES domain fits within the south-western quadrant of the reference simulation domain. Two subregions presented later in Section 4, occupy the same physical location in every LES. Given easterly winds on this simulation date, the reference domain has significantly more fetch than the smaller LES at these locations. The d03_30s_ref simulation requires about 5 times more CPU-hours compared to the d03_3s_cpm simulation. Meanwhile the computational cost difference between same sized domains with or without CPM is entirely negligible.

4. Results and Discussion

4.1. Mean Flow Fields

Instrumentation deployed during the Perdigão field campaign included GPS Advanced Upper-Air Sounding System (GAUSS) radiosondes as well as 195 three-component sonic anemometers and 55 temperature-humidity sensors mounted on masts of heights ranging from 2–100 m above ground level (a.g.l.) [17]. These towers sampled at 20 Hz but we analyze data averaged to 5 min intervals for comparison of the WRF model to a subset of the tower data from masts at 2 m, 10 m, and 100 m above ground level. The focus of the comparisons to field data is on mean flow fields, wind and temperature data, to provide a basic validation of model performance and context for the investigation of turbulence in the following section. The simulated time period is 20 May 2017. This was a dry, clear day with high clouds and a low pressure system setting up to the southwest of Portugal. Surface winds veered from northeasterly in the morning to easterly in the afternoon, and veered back to northeasterly over the evening. Wind speeds were low, <10 m s${}^{-1}$, until night time, when those accelerated to moderate, <15 m s${}^{-1}$, values.

This section describes comparisons of model results to field data from radiosonde and tower data from the valley ridges. The towers within the valley are not included because the valley is unresolved in simulations with standard 30 arcsec topography. For much of the extent of the valley, the two ridges are merged into a single ridge (Figure 3). Error statistics at 2 m and 10 m a.g.l. presented in

Table 2. Biases and root-mean-square errors (RMSE) averaged over the time period between 1100 UTC on 20 May 2017 and 0000 UTC on 21 May 2017 comparing LES model fields interpolated to the locations of towers, excluding those within the unresolved valley. Biases are averages of the difference taken as model data minus observed tower data. The LES data at 2.5 min output interval is first downsampled as described in the body text.

	T2 [K]		U10 [m s${}^{-1}$]		$\mathit{\varphi }$10 [°]		TKE10 [m${}^2$ s${}^{-2}$]
Simulation	Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE
d03_30s	−0.45	1.15	−1.28	2.52	0.2	25.2	−0.87	1.27
d03_30s_cpm	−0.46	1.14	−1.29	2.53	1.2	25.3	−0.86	1.26
d03_30s_ref	−0.34	1.11	−1.29	2.54	2.9	24.6	−0.86	1.25
d03_3s	−0.69	1.41	0.39	1.84	2.0	24.6	−0.66	1.09
d03_3s_cpm	−0.70	1.41	0.37	1.86	1.4	23.7	−0.65	1.09

reflect errors calculated across 17 towers on the ridges only. These are known as towers tNW03, tNW11, tSE04, tSE13, rSW01, rSW02, rSW03, rSW04, rSW05, rSW06, rSW07, rSW08, rNE01, rNE02, rNE03, rNE06, and rNE07 using the naming convention from the field campaign. The LES data at 2.5 min output interval is first downsampled to the 5 min averaging interval of the tower data using the trapezoidal rule. Adopting the convention that $u_i^{2n}$ corresponds to the velocity field after 2n2.5 min simulated time, the downsampling technique for velocity is

$$\overline{u_i^{2n}}=\frac{1}{4}\left(u_i^{2n-1}+2u_i^{2n}+u_i^{2n+1}\right)$$

(5)

for integer, n. Potential temperature is averaged with the same approach as with velocity. Turbulent kinetic energy is calculated as

$$TK{E}^{2n}=\overline{TK{E}_{\mathrm{sg}}^{2n}}+\frac{1}{2}\sum _{i=1}^3\frac{1}{3}\sum _{m=-1}^1{\left(\overline{u_i^{2n}}-u_i^{2n+m}\right)}^2$$

(6)

where the first term represents average subgrid TKE, prognosticated by the LES closure [41], and the second term represents the resolved turbulence. Heights of 2 m and 10 m a.g.l. are used because they are commonly reported for model validation. Specifically, temperature taken at 2 m, T2, along with wind variables taken at 10 m, U10 for wind speed, $\varphi$10 for wind direction, and TKE10 for turbulent kinetic energy, are reported. Additional instrumentation was also installed at 100 m a.g.l. for two of these towers, tSE04 and tSE13. These two towers are used to calculate the error statistics in

Table 3. As in Table 2 but for data at 100 m above ground level.

	T100 [K]		U100 [m s${}^{-1}$]		$\mathit{\varphi }$100 [°]		TKE100 [m${}^2$ s${}^{-2}$]
Simulation	Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE
d03_30s	−0.45	0.62	−0.86	1.55	−9.8	17.8	−0.40	0.90
d03_30s_cpm	−0.44	0.61	−0.90	1.55	−9.5	16.5	−0.35	0.85
d03_30s_ref	−0.38	0.61	−0.86	1.67	−7.4	15.6	−0.32	0.97
d03_3s	−0.87	0.93	0.12	1.04	−9.8	16.9	−0.25	0.95
d03_3s_cpm	−0.86	0.93	0.13	0.97	−9.8	16.2	−0.25	0.83

at 100 m a.g.l., which is closer to typical hub heights of large wind turbines, and also avoids near-surface extrapolation errors in the WRF output. Extending the naming convention above, these variables are abbreviated as T100, U100, $\varphi$100, and TKE100 for temperature, wind speed, wind direction, and turbulent kinetic energy at 100 m a.g.l., respectively. Time series of these variables at tSE04 are shown in Figure 4 and Figure 5.

Varying the resolution of the topography results in larger differences across models for these variables than application of CPM does. Notably, the bias of wind speed flips sign depending on the topography resolution (Table 2) due to a bifurcation of models during the ramp up of wind speed seen during the evening, when winds increase from ≈6 m s${}^{-1}$ at 1700 UTC to ≈12 m s${}^{-1}$ by 2030 UTC at tSE04 and 100 m a.g.l. (Figure 4). After this ramp up, simulations with 30 arcsec topography underestimate the wind speeds by ≈2 m s${}^{-1}$ while those with 3 arcsec topography overestimate wind speed by ≲1 m s${}^{-1}$ for the remainder of the night. The temperature time series in Figure 5 also show more sensitivity to topography resolution. In this case it is the finer resolved topography models that have a greater magnitude bias earlier in the simulated period. All models develop an approximately 2 K cold bias throughout the daytime that persists through the nocturnal period. The TKE time series show that the TKE from the model is generally lower than that measured from the towers. Using CPM and finer topography do not have a strong effect on the TKE far from the boundary, at the location of the towers, but the biases are slightly reduced overall.

Model results are also compared to the GAUSS radiosonde observations from soundings released on the valley floor at approximately 1115 UTC (Figure 6), 1715 UTC (Figure 7), and 2315 UTC (Figure 8). These times of day represent three different atmospheric stability regimes: strongly convective, weakly convective, and weakly stable, respectively. All the LES domains appear to perform similarly, as confirmed by error statistics (

Table 4. Biases and root-mean-square errors (RMSE) comparing LES results to GAUSS radiosonde data from soundings released from the center of the Vale do Cobrão (red ‘x’ in Figure 2) at 11:13:11 UTC, 17:16:25 UTC, and 23:13:28 UTC. Biases are averages of the difference taken as model data minus observed radiosonde data interpolated to the heights of the LES model levels below 5 km above sea level.

Simulation		d03_30s		d03_30s_cpm		d03_30s_ref		d03_3s		d03_30s_cpm
		Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE
U [m s${}^{-1}$]	1113	−1.21	2.09	−1.13	2.04	−1.18	2.05	−1.11	1.94	−1.19	2.06
$\varphi$ [°]	1113	−2.71	29.32	−3.05	33.77	1.09	28.99	−7.12	37.18	−9.64	38.04
$\theta$ [K]	1113	−0.63	0.85	−0.62	0.87	−0.64	0.81	−0.55	0.68	−0.57	0.82
Q${}_v$ [g kg${}^{-1}$]	1113	−0.33	0.68	−0.23	0.64	−0.08	0.48	−0.19	0.59	−0.19	0.57
U [m s${}^{-1}$]	1716	0.29	1.87	0.53	1.62	0.39	1.32	0.20	1.34	0.20	1.39
$\varphi$ [°]	1716	−0.97	27.21	−2.12	25.65	−3.34	21.75	1.19	23.19	4.68	25.55
$\theta$ [K]	1716	−0.92	1.25	−1.03	1.26	−0.91	1.36	−0.87	1.23	−0.95	1.21
Q${}_v$ [g kg${}^{-1}$]	1716	0.25	0.66	0.28	0.57	0.14	0.69	0.22	0.72	0.27	0.64
U [m s${}^{-1}$]	2313	−0.76	2.27	−0.76	2.27	−0.87	2.34	−0.62	2.05	−0.64	2.05
$\varphi$ [°]	2313	−0.21	6.90	-0.39	7.10	0.47	6.83	4.29	26.15	3.98	25.77
$\theta$ [K]	2313	−1.31	1.46	−1.31	1.46	−1.27	1.41	−1.17	1.31	−1.17	1.31
Q${}_v$ [g kg${}^{-1}$]	2313	0.25	0.66	0.25	0.66	0.24	0.66	0.16	0.64	0.16	0.64

). For example, biases for wind speed are consistently under 1.5 m s${}^{-1}$, though root-mean-square errors (RMSE) can be larger in part due to fine-scale oscillations not resolved by the model. Similarly, biases of wind direction are quite low, less than 10${}^{°}$ in magnitude, though root-mean-square errors are higher, ≈30${}^{°}$ across model configurations. In these profiles, only the simulations with the finer 3 arcsec topography input include the strong veer observed at the lowest elevations. For simulations with coarser 30 arcsec topography, the solution is not even defined until nearer the altitude marked ‘Ridge height’ on the sounding profiles in Figure 6, Figure 7 and Figure 8. This height is given as 473 m above sea level (a.s.l.), corresponding to the base of tSE04, the tower whose time series from 100 m above this height are shown in Figure 4 and Figure 5. As in the tower data, all models show a cold bias at each sounding time and tend to underpredict the boundary-layer depth during the daytime. Despite this, the general diurnal cycle is qualitatively captured well. At 1113 UTC, below a capping inversion around 1 km a.s.l., there is a well mixed layer with nearly due easterly flow. By 1716 UTC, the capping inversion is around 2 km a.s.l and the well mixed layer has grown in the profiles of temperature and water vapor while the wind fields are similar. The wind speeds do increase, to a maximum of 17 m s${}^{-1}$, near the bottom of a residual layer around 1 km a.s.l at 2313 UTC. Temperature profiles are generally stable below the residual layer at this later time. Though all the models generally perform similarly, the profiles of the nocturnal boundary layer at 2315 UTC show particularly notable agreement across models.

Results in Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 and in Table 2 and Table 4 show that the WRF simulations compare similarly well to the radiosonde and surface tower data for all cases. Notably, the simulations with and without CPM are similar in the time series, profiles, and error statistics considered so far while more significant difference are due to the topography resolution. Negligible differences across models may seem to suggest that the CPM has no discernible effect on the simulations. It is important to recognize that so far the variables used to validate the model are all mean fields, taken far from the domain lateral boundaries. In the following section, examining the model’s representation of turbulence, it is shown that the CPM does indeed affect simulations in important ways by reducing the fetch required for turbulence to fully develop.

4.2. Turbulence Structure

With 150 m grid spacing, some of the eddies in the inertial subrange are explicitly resolved in the LES. Differences in the resolved turbulence produced by the different model configurations can be quite striking, despite the differences in the mean flow fields examined above being relatively small. The following sections examine the turbulent structure of the LES simulations over the Perdigão domain at various times during the day. The analysis is focused on three time periods coincident with the three radiosonde releases discussed in the previous section. Each time is representative of a different stability class: strongly convective (1115 UTC), weakly convective (1715 UTC), and weakly stable (2315 UTC). As seen in the radiosondes (Figure 6, Figure 7 and Figure 8), the mean horizontal flow in the boundary layer is easterly above the ridge with ∼200 m relief. The model’s seventh vertical level, at ∼200 m a.g.l., lies within this region of easterly flow though its elevation a.s.l. varies with the terrain. This level is selected for the investigation of turbulence through contours and spectra from the vertical velocity field. In LES, the resolved component of vertical velocity is nearly equal to its turbulent fluctuation, because the spatial average of vertical velocity tends to zero as the average is taken over larger areas. The contours of vertical velocity can therefore be used to visually assess part of the turbulent kinetic energy of the flow. Turbulent structures appear as regions of higher magnitude vertical velocity, updrafts and downdrafts.

While the contours of vertical velocity are helpful to get a qualitative sense of the effects of the CPM, a more quantitative analysis is aided by computing the kinetic energy spectra. For this purpose, one-dimensional spectra are computed from N points of $u_i$ ($u_1=u$ for zonal, $u_2=v$ for meridional, or $u_3=w$ for vertical velocity) on a transect of length $L=N\mathsf{\Delta }s$ in the s-direction ($s=x$ for West–East, or $s=y$ for South–North). In example, for w, following Durran et al. [42], the discrete spectral density, $E_{ws}$, is defined such that

$$\frac{1}{2}{\langle w^2\rangle }_s=\frac{1}{L}\sum _{j=1}^N\frac{w_j^2}{2}\mathsf{\Delta }s=\sum _{m=1}^{N/2+1}{E}_{ws,m}\mathsf{\Delta }{k}_s$$

(7)

where $\mathsf{\Delta }{k}_s=\frac{2\pi }{L}$, and ${\langle ·\rangle }_s$ indicates averaging in the s-direction throughout the manuscript. To achieve this normalization, the discrete spectral densities are calculated as

$$E_{ws,m}=\frac{\mathsf{\Delta }s{\left|{\widehat{w}}_m\right|}^2}{2\pi N\left(1+{\delta }_{m,\frac{N}{2}+1}\right)}$$

(8)

where $\widehat{w}$ is the discrete Fourier transform of w following the normalization conventions of numpy’s fft or rfft routine, and $\delta$ is the Kronecker delta operator [42]. In this normalization convention, the forward transform is not normalized such that

$${\widehat{w}}_m=\sum _{j=1}^N{w}_{j}exp\left[\frac{-2\pi i}{N}(m-1)(j-1)\right].$$

(9)

4.3. Weakly Convective Case

At 1715 UTC, turbulent structures are visually apparent in the contours of vertical velocity in Figure 9 and Figure 10. Each of the simulations contains a region of smoother flow near the eastern inflow boundary, where vertical velocity is nearly zero and larger-scale structures persist over longer distances than they do in the rest of the domain. From the vertical velocity contours, it appears that finer scale turbulence develops consistently only after ≈20 km of fetch from the eastern boundary in d03_30s_ref. That is, turbulence appears fully developed only over $x\in$ [0,50] km, which excludes the eastern 22 km in the domain. Similarly, the larger structures diminish only after about 16 km of fetch in d03_30s for this easterly flow case, i.e., limiting the region with well developed turbulence to $x\in$ [0,20] km. There are some regions which develop particularly streaky velocity structures, apparent in the south-east of d03_30s_ref, $x\in$ [60,70] km and $y\in$ [0,30] km. These streaky structures have been observed in numerous previous studies [12,43,44]. While convective rolls likely do develop in weakly convective and sheared boundary layers, they are also likely over-represented in numerical simulations. The presence of these rolls indicates under-developed model turbulence that influences local turbulence characteristics and shear related mixing [45]. In flat terrain cases, these structures can persist up to 54 km before they begin to dissipate [3,4]. In the example here, most of the largest streaky structures break up as the flow passes over the rolling hills in the domain. Furthermore, these structures do not occur in the same physical locations in the reference domain as they do in the smaller simulations. Rather, they are distinctly tied to the distance or fetch from the inflow boundary.

One way to quantify the fetch required for turbulence to fully develop is to compute the spectral densities, defined above, at varying distances from the inflow boundary. Following Muñoz-Esparza et al. [18], Figure 11 computes such spectra from the vertical velocity field shown in Figure 9a over N = 128 points in the South–North direction with $y\in$ [8.4,27.6] km to evaluate the particularly ‘streaky’ region noted above. A resolved inertial subrange, over which spectra are nearly parallel to a −5/3 power law line, is expected to develop with increasing fetch as more turbulent kinetic energy is produced. We also expect energy spectral densities to be increasing over farther fetch as turbulence develops, as seen over the first 40 points of fetch, after which the signal is less clear over the next 40 points of fetch. Measured this way, the spectra suggest turbulent energy at 60 points of fetch is more intense than at 80 points of fetch. This region of the domain, however, is precisely the region over which turbulent streaks seem to be breaking up in the velocity contours. During these easterly wind conditions, spectra computed in the y-direction therefore do not seem to reliably capture what is observed in the contours to be the transition from larger scale streaks to finer scale turbulent structures. As the rolls are oriented with a spin around the x-axis and thus produce strong oscillations in the y-direction, the correspondingly large y-direction spectra are misleading in the current context.

To remedy this, another method of measuring the turbulence is required. Simply computing two-dimensional spectra is one possibility that has been tested (not shown). However, computing one-dimensional spectra in the x-direction produces a clearer signal of the change in turbulence over fetch, because the roll structures are enlarged in the x-direction given easterly winds. For this reason, we prefer a pseudo-two-dimensional technique in which spectra are calculated for N = 64 points in the x-direction before being averaged over 180 y-direction indices, and over 13 time indices (30 min at 2.5 min intervals). The horizontal velocities have their linear tendencies, from the difference between the 1st and 64th points, removed before Fourier decomposition. Two subregions are selected to compute spectral density at different distances from the inflow boundary. One subregion, outlined with solid black lines in Figure 9 and all following contour plots, has 40 points of fetch before the first points in d03_30s, d03_30s_cpm, d03_30s, and d03_30s_cpm. This region will be referred to as the short fetch subregion. The other subregion, referred to as the long fetch subregion, is outlined with a dashed black line and has 120 points of fetch for the same LES domains.

The pseudo-two-dimensional energy spectra for each velocity component are computed at the short fetch subregion and long fetch subregion (Figure 12). At both fetches, the velocity fields are taken at 200 m a.g.l. and averaged over 30 min centered around 1715 UTC, the time plotted in Figure 9 and Figure 10 when weakly convective conditions are observed. The efficacy of CPM is demonstrated at the short fetch subregion, where we have the benefit of comparing spectra from near the boundary of the small domains to d03_30s_ref which allows for a much longer fetch before the first transect location. Recall, these smaller LES domains fit within a quadrant of the reference simulation, d03_30s_ref. The regions are drawn for the same physical locations as in the smaller LES, but the reference domain will have an additional fetch of 240 grid points.

As we may expect from the contour plots, short fetch spectra from d03_30s_cpm and d03_3s_cpm agree well with those from d03_30s_ref for a range of wavenumbers, e.g., $k\in [4·{10}^{-3},1.1·{10}^{-2}]$ m${}^{-1}$ for vertical velocity, over which simulations without CPM are missing energy. Though the addition of higher resolution topography, as in d03_3s, increases the energy compared to its low resolution topography counterpart, d03_30s, their spectral densities are missing energy at these wavenumbers compared to the reference simulation spectra. These wavenumbers are not entirely in the resolved inertial subrange, over which the model spectra are nearly parallel to the −5/3 slope line, and include wavenumbers just greater than those of the resolved interial subrange. Particularly for vertical velocity spectra, at the lowest wavenumbers the smaller LES actually have more energy than the reference simulation, perhaps due to effects of the boundary or the presence of the streaks. Similarly, the vertical kinetic energy of the very highest wavenumbers are increased in the smaller LES compared to the reference simulations. This may suggest that the turbulent energy cascade is still developing at the high wavenumbers for the smaller simulations. For the high resolution topography simulations, some of this energy at the highest wavenumbers is likely due to the effects of the terrain. Especially for the horizontal velocity spectra after long fetch, the high resolution topography runs have increased energy at the highest wavenumbers, even as the rest of the spectra collapse (Figure 12). The eventual collapse of spectra from LES with and without CPM suggests the application of CPM does not affect the long fetch characteristics of the turbulence while it does accelerate the development of turbulence. For the weakly convective conditions observed around 1715 UTC, simulations with CPM developed turbulence over 40 points of fetch that is similar to that from a simulation without CPM but with 280 points of fetch.

Muñoz-Esparza and Kosović [5] show the importance of the CPM can be predicted by the ratio of characteristic wind speed at the top of the capping inversion, $U_{ci}$, to a convective velocity scale,

$$w_{*}={\left(\frac{g}{\theta }{z}_i{〈w^{\prime }{\theta }^{\prime }〉}_{\mathrm{sfc}}\right)}^{1/3},$$

(10)

where ${〈w^{\prime }{\theta }^{\prime }〉}_{\mathrm{sfc}}$ is the surface heat flux. They find that for $U_{ci}/w_{*}>5$, the CPM will be especially beneficial. In the current weakly convective case, with $U_{ci}=$ 10 m s${}^{-1}$, $z_i=2$ km and ${〈w^{\prime }{\theta }^{\prime }〉}_{\mathrm{sfc}}=0.05$ K m s${}^{-1}$, the ratio, $U_{ci}/w_{*}=7$, and the expected effects of CPM are apparent, as predicted, at this time. As this time was selected on the basis of the radiosonde launch schedule, and not on maximizing the ratio, $U_{ci}/w_{*}$, the CPM may be even more beneficial later in the evening transition as the surface heat flux decreases.

In summary, the x-direction one-dimensional spectra, computed over short fetch and long fetch regions, demonstrate a clear advantage to using CPM to trigger appropriate length scales. Instead of requiring 280 grid points to develop turbulence under standard grid nesting configurations, nested WRF simulations using CPM can develop over only about 40 points to an equivalent level of turbulence. There is no detrimental effect seen from adding the perturbations, as the level of turbulence at long fetches is unaffected and mean flow differences are minimal.

4.4. Strongly Convective Case

The benefits of the CPM are clearly demonstrated during the evening transition time period selected above. Other times of day show different and sometimes less impact of the CPM due to the effects of atmospheric stability. At 1115 UTC, stronger convection is evidenced by smaller thermal cells in the vertical velocity contours (comparing Figure 13 and Figure 14 to Figure 9 and Figure 10). Under strong convective conditions the flow is seen to develop considerably smaller scales, which challenges the numerical resolution and accuracy by introducing additional high-pass dynamics in the computational model [46]. During this time, the fetch required for streaky structures to develop into fully convective turbulence appears to drop to approximately 12 km in d03_30s_ref and only 6 km in the smaller domains, as seen in the contours of vertical velocity. This is confirmed by the vertical velocity power spectra (Figure 15) of the smaller LES that generally agree with the spectrum from d03_30s_ref, even after limited fetch. Examining power spectra associated with the other two components of velocity, we can see some benefit of the CPM method, but this effect is small compared to that of the weakly convective case. Corroborating the hypothesis of Rai et al. [12], the combined influence of complex terrain and convection are sufficiently strong that finer scale turbulence develops quickly after the inflow enters the domain regardless of CPM being applied or not. The lessened effects of CPM in this strongly convective case compared to the previous weakly convective case are related to the ratio of characteristic wind speed to convective velocity scale. In the strongly convective case, with $z_i=1$ km and ${〈w^{\prime }{\theta }^{\prime }〉}_{\mathrm{sfc}}=0.3$ K m s${}^{-1}$, the ratio, $U_{ci}/w_{*}=2.8$, is less than the threshold, $U_{ci}/w_{*}\lesssim 3$, proposed by Muñoz-Esparza and Kosović [5]. This confirms for a real terrain case what was observed in ideal cases, that the CPM is less impactful when the convective velocity scale is sufficiently large relative to a characteristic wind speed. Furthermore, during these strongly convective conditions, the use of higher resolution topography does not affect the spectra as much.

4.5. Weakly Stable Case

At night (2315 UTC), the reference simulation, d03_30s_ref, does not generate apparent small scale turbulent motions even over its extensive fetch; neither do the smaller LES cases. Rather, the velocity contours (Figure 16 and Figure 17) are dominated by what appear to be internal gravity wave structures with alternating streaks of negative and positive vertical velocity. These streaks are aligned in the spanwise direction (and are not aligned in the streamwise direction as are the spurious convective rolls or streaks seen earlier). Simulations with the same topography resolution result in similar spectra (Figure 18) at the short fetch subregion regardless of the application of CPM. This is similar to the strongly convective case in that the application of CPM does not significantly change the model turbulence. Conversely, varied topography results in differences between spectra that persist to the long fetch subregion more notably than in the convective cases. The higher resolution topography simulations have much higher spectral density at all but the lowest wavenumbers. This is reflected in small scale structures apparent in the vertical velocity contours of the high resolution topography cases (Figure 17) that are missing in the low resolution topography LES including the reference simulation (Figure 16). However, this increased energy is likely more reflective of increased wave energy in the mean flow field, rather than turbulent energy, because no inertial subrange develops.

Using higher resolution topography and extending the domain fetch as in d03_30s_ref do not help the model produce an inertial subrange in this stable case. This indicates that the $\mathsf{\Delta }x$ = 150 m resolution may be insufficient for the stable conditions described here. The CPM does not provide significant benefit for simulations of the nocturnal boundary layer at this coarse resolution. Increasing the resolution of the input topography is recommended for developing the fine scale features of the nocturnal boundary layer, but, at the current resolution, these features are most likely associated with the mean gravity waves due to the topography itself rather than with resolved turbulence.

5. Summary and Conclusions

The cell perturbation method (CPM) is used here to accelerate the development of realistic turbulence in nested large-eddy simulations performed using the WRF model over complex terrain. This is, to the authors’ knowledge, the first study in which the effects of CPM are examined for real case simulations with highly complex terrain conditions. The results indicate that CPM is particularly helpful during certain time periods and has an overall positive impact on the solution. The use of higher resolution topography is also explored as a way to develop more realistic turbulence in the nested simulations, in addition to improvements made from applying the CPM.

Nested simulations with and without CPM are performed over the site of the Perdigão field campaign and are compared to a reference simulation with a larger domain. Even with complex terrain, all simulations during convective conditions show spurious streaky structures near the inflow boundaries which persist over a large fraction of the domain. In the simulations without CPM, these spurious structures break down into finer scale turbulence after ≈10–20 km of fetch. This distance is less than the typical distances reported (as much as 54 km of fetch) for flat terrain cases [2,3,4,5,11], but is nonetheless a significant portion of the domain which is lost to spurious flow features. Using CPM shortens the fetch required for turbulence to fully develop. As an example, the fetch required to develop turbulence at 1715 UTC (Figure 9) was reduced from 16 km to 6 km using CPM. This recovery of over 27% of the entire domain size is achieved through CPM with negligible computational cost. If a larger reference simulation were used instead to generate the same level of turbulence, as illustrated with d03_30s_ref, it could require ∼5 times more CPU-hours than a smaller domain using CPM.

Unlike CPM, the use of more finely resolved model topography does not generally accelerate the generation of the fine-scale turbulence desired on the inner domain. Refining the topography can increase energy in the spectral densities of velocity, but this effect is most notable under the weakly stable stratification and is not apparent during sufficiently strong convective conditions. For the weakly convective conditions, for example, increasing the topography resolution is not sufficient for fully developing turbulence on the inner domain after 40 points, or 6 km, of fetch. Application of CPM at this time, regardless of topography resolution, leads to agreement with the reference simulation at the same location. The longest fetches required to fully develop the energy spectra occur during this weakly convective time period and therefore set a minimum domain size for simulations that might include a full diurnal cycle. CPM has more impact than high-resolution topography on turbulent structures during this weakly convective time, making CPM much more cost effective. Furthermore, while there are of course important reasons to implement more finely resolved topography whenever possible, e.g., to resolve a narrow valley such as the Vale do Cabrão, greater benefits are found from implementing the CPM at the same time.

The influence of the CPM is found to depend on the time of day, which reflects a dependence on atmospheric stability. For the convective cases, this dependence is characterized by a ratio of a characteristic wind speed and a convective velocity scale, $U_{ci}$/$w_{*}$[5]. The benefits of CPM are made most clear under weakly convective conditions, such as at 1715 UTC, when $U_{ci}$/$w_{*}$=7 as discussed above. At this time, the simulations using CPM have greater turbulent energy than simulations without CPM which are missing energy compared to the reference simulation. During more strongly convective conditions, at 1115 UTC, the effect of CPM is less pronounced when $U_{ci}$/$w_{*}\lesssim$ 3, which closely agrees with previous reports by Muñoz-Esparza and Kosović [5]. At this time, all simulations develop turbulence and spectra resembling those from the reference simulation soon after the inflow boundary. The combined effects of convection and heterogeneous topography may render CPM unnecessary during these strongly convective conditions (as noted by [12]). With weak stability, such as at 2315 UTC, the use of finer topography resolution seems to be more significant than CPM in generating finer scale structures as seen in the model energy spectra. The stable simulation period would possibly benefit from transition to the Richardson’s number-based CPM approach [5], and more likely from even finer model resolution. Extending the current multiscale approach, additional nests with finer resolved grid meshes should be used for even weakly stable conditions. As finer resolved models would be nested within simulations resembling the coarse LES presented here, there is a practical interest in reducing spurious information on these intermediate coarse domains. We have shown that use of higher resolution topography is a means to increase energy found in velocity spectra in the current complex terrain case, though finer nests are likely to benefit from the application of CPM, as in Arthur et al. [19]. Since the current work does not pursue these finer nests, the role of CPM under stable conditions with complex terrain needs further investigation.

In conclusion, using CPM during the entire simulation period can lead to improvement during certain conditions (e.g., weakly convective) and does not have any detrimental effects during other times of day. Looking forward, applications such as wind energy that rely on weather forecasting will continue to increase in resolution in response to increasing computational power. Eventually, simulations that can resolve narrow valleys and even turbine wakes, such as those of interest to the Perdigão project, will be nested in simulations such as those presented in the current work. As shown here, the CPM is a cost-effective method to generate realistic turbulence without requiring very large domains to provide excessively large fetches from inflow boundaries. Under the appropriate conditions, the benefits of CPM are substantial even over complex, heterogeneous terrain. With essentially no drawbacks in computational cost or the quality of the developed flow, CPM is a promising method for maintaining low fetches and improving the representation of turbulence in future nested simulations.