Implementation of a synthetic inflow turbulence generator in idealised WRF v3.6.1 large eddy simulations under neutral atmospheric conditions

A synthetic inflow turbulence generator was implemented in the idealised Weather Research and Forecasting large eddy simulation (WRF-LES v3.6.1) model under neutral atmospheric conditions. This method is based on an exponential correlation function and generates a series of two-dimensional slices of data which are correlated both in space and in time. These data satisfy a spectrum with a near “$-\mathrm{5}/\mathrm{3}$” inertial subrange, suggesting its excellent capability for high Reynolds number atmospheric flows. It is more computationally efficient than other synthetic turbulence generation approaches, such as three-dimensional digital filter methods. A WRF-LES simulation with periodic boundary conditions was conducted to provide prior mean profiles of first and second moments of turbulence for the synthetic turbulence generation method, and the results of the periodic case were also used to evaluate the inflow case. The inflow case generated similar turbulence structures to those of the periodic case after a short adjustment distance. The inflow case yielded a mean velocity profile and second-moment profiles that agreed well with those generated using periodic boundary conditions, after a short adjustment distance. For the range of the integral length scales of the inflow turbulence (±40 %), its effect on the mean velocity profiles is negligible, whereas its influence on the second-moment profiles is more visible, in particular for the smallest integral length scales, e.g. those with the friction velocity of less than 4 % error of the reference data at $x/H=\mathrm{7}$. This implementation enables a WRF-LES simulation of a horizontally inhomogeneous case with non-repeated surface land-use patterns and can be extended so as to conduct a multi-scale seamless nesting simulation from a meso-scale domain with a kilometre-scale resolution down to LES domains with metre-scale resolutions.

1 Introduction

Atmospheric boundary layer flow involves a wide range of scales of eddies, from quasi-two-dimensional structures at the mesoscales to three-dimensional turbulence (normally with higher Reynolds number, i.e. Re∼10⁸–10⁹) at the microscale (Muñoz-Esparza et al., 2015). The Weather Research and Forecasting (WRF) model (Skamarock and Klemp, 2008) is capable of simulating atmospheric systems at a variety of scales. At the mesoscale and synoptic scales, the WRF model allows grid nesting for downscaling from 10–100 to 1–10 km using a fully compressible and non-hydrostatic Reynolds-averaged Navier–Stokes (RANS) solver (Skamarock and Klemp, 2008), which captures the behaviour of mean flows only. At the microscale, a large eddy simulation (LES) can be activated in the WRF model (WRF-LES), enabling users to simulate the characteristics of energy-containing eddies in the atmospheric boundary layer. Challenges remain in downscaling from the mesoscale (resolutions down to 1 km, capturing mean information only) to the LES scale (tens of metres or below, capturing additional turbulence information) (Doubrawa et al., 2018; Talbot et al., 2012; Chu et al., 2014; Liu et al., 2011), such as specifying the appropriate inflow conditions for an LES domain and the sub-grid-scale turbulence schemes for the “grey-zone” resolution, to which neither planetary boundary layer (PBL) nor LES parameterisation schemes can apply well. Consequently, microscale and mesoscale flows are typically studied separately. Most LES models of atmospheric boundary layer flow at the microscale use periodic boundary conditions and simplified large-scale geostrophic forcing for idealised simulations. However, the use of periodic boundary conditions implicitly assumes that atmospheric fields and the underlying land use have repeated periodic features. This assumption may be unrealistic for real landscapes where land-use patterns and the atmospheric phenomena coupled to them can be very heterogeneous. Therefore, such periodic WRF-LES simulations are restricted to studies of the atmospheric boundary layer flow with a single domain (e.g. Zhu et al., 2016; Kirkil et al., 2012; Kang and Lenschow, 2014; Ma and Liu, 2017) or the outermost domain of either one-way nested cases (e.g. Nunalee et al., 2014) or two-way nested cases (e.g. Moeng et al., 2007). Here we implement a well-tested synthetic turbulence inflow scheme (Xie and Castro, 2008) in the WRF-LES model (v.3.6.1), in which the meso-scale model could provide the mean flow information as the input of the synthetic turbulence inflow scheme. This scheme provides a step towards enabling WRF's capability of nesting micro-scale turbulent flows within realistic meso-scale meteorological fields.

Dhamankar et al. (2018) reviewed three broad classes of methods to generate the turbulent inflow conditions for LES models, mainly for engineering applications. The first class is the library-based method, which relies on an external turbulence library to provide inflow turbulence. The turbulence library can be based on either (a) the precursor or concurrent simulation (e.g. Munters et al., 2016) on the same geometry to a main LES simulation; or (b) a pre-existing database (e.g. Schluter et al., 2004; Keating et al., 2004) from experiments or computations (on a different geometry to a main LES simulation). Although this method is usually limited to specialised applications, it can provide good-quality inflow turbulence. The second class is the recycling–rescaling-based method (e.g. Lund et al., 1998; Morgan et al., 2011), in which the velocity field is recycled from some suitably selected downstream plane back to the inflow boundary plane. Although this method may be effective in producing well-established turbulence, there are some limitations, e.g. the requirements of an equilibrium region near the inlet and a relatively large domain. The turbulence profile determined by the geometry of the precursor simulation can be added on the top of any given mean profile, which could be modified and varied in time for more realistic applications. The third class is the synthetic turbulence generator, which includes a variety of methods such as the Fourier transform-based method (e.g. Kraichnan, 1970; Lee et al., 1992), proper orthogonal-decomposition-based method (e.g. Berkooz et al., 1993; Kerschen et al., 2005), digital-filter-based method (e.g. Klein et al., 2003; Xie and Castro, 2008; Kim et al., 2013), diffusion-based method (e.g. Kempf et al., 2005), vortex method (e.g. Benhamadouche et al., 2006) and synthetic eddy method (e.g. Jarrin et al., 2006). The synthetic turbulence generator has the potential to be used for a wide range of flows. Due to the imperfection of the synthetic turbulence, which is not directly derived from generic flow equations, these methods normally require some inputs and a certain adjustment distance for turbulence to become well-established. For more information about the above synthetic turbulence generation methods, we recommend Tabor and Baba-Ahmadi (2010), Wu (2017), and Bercin et al. (2018).

Several other methods have been developed to generate inflow turbulence for atmospheric boundary layer flow in nested WRF-LES models. Mirocha et al. (2014) introduced simple sinusoidal perturbations to the potential temperature and horizontal momentum equations near the inflow boundaries. This method can speed up the development of turbulence and generally has a satisfactory performance in the nested WRF-LES domains, providing promising results. Muñoz-Esparza et al. (2014) extended the perturbation method of Mirocha et al. (2014) and proposed four methods, i.e. the point perturbation method, cell perturbation method, spectral inertial subrange method, and spectral production range perturbations, to generate perturbations of potential temperature for a buffer region near the nested inflow planes. The cell perturbation method was found to have the best performance regarding the adjustment distance for the turbulence to be fully developed. It has the advantages of negligible computational cost, minimal parameter tuning, not requiring a priori turbulent information, and efficiency to accelerate the development of turbulence. Muñoz-Esparza et al. (2015) further generalised the cell perturbation method of Muñoz-Esparza et al. (2014) under a variety of large-scale forcing conditions for the neutral atmospheric boundary layer. The perturbation Eckert number (describing the interaction between the large-scale forcing and the buoyancy contribution due to the perturbation of potential temperature) was identified as the key parameter that governs the transition to turbulent flow for nested domains. They found an optimal Eckert number to establish a developed turbulent state under neutral atmospheric conditions. These methods impose temperature perturbations at specific length scales and timescales related to the highest resolved wave number in the LES. It was demonstrated in Muñoz-Esparza et al. (2015) that a distance of about 15 boundary layer depths is required to allow the flow to be fully turbulent when the temperature perturbation method is adopted in the one-way nesting WRF model. It is to be noted that the temperature perturbation method was introduced for mesoscale to microscale coupling approach where smooth mesoscale flow (no resolved turbulence) forces microscale flow by using the one-way nesting approach in WRF. Muñoz-Esparza et al. (2014) stated “the perturbation method is to provide a mechanism that accelerates the transition towards turbulence, rather than to impose a developed turbulent field at the inflow planes as the synthetic turbulence generation methods pursue”, and “the use of temperature perturbations presents an alternative to the classical velocity perturbations commonly used by most of the techniques.” The optimisation and generalisation of these methods would require intensive testing. Muñoz-Esparza and Kosovic (2018) extended the cell perturbation method of the inflow turbulence generation to non-neutral atmospheric boundary layers. Instead of adopting temperature perturbations in the original cell perturbation method, Mazzaro et al. (2019) further explored the random force perturbation method (FCPM) in the multiscale nested domains.

Due to its accuracy, efficiency, and, in particular, the capability for high Reynolds number flows, the synthetic inflow turbulence generator (Xie and Castro, 2008) has been implemented and tested on codes developed for engineering applications, such as Star-CD (Xie and Castro, 2009) and OpenFOAM (Kim and Xie, 2016), and the micro-scale meteorology code PALM (PALM, 2017; Maronga et al., 2020). This study focuses on an implementation of this synthetic inflow turbulence generator (Xie and Castro, 2008) in the idealised WRF-LES (v3.6.1) model under neutral atmospheric conditions. In this paper, Sect. 2 describes the methodology of the WRF-LES model and the technique of the synthetic inflow turbulence generator, Sect. 3 presents the results of the WRF-LES model with the use of the synthetic inflow turbulence generator, and Sect. 4 states the conclusions and future work.

2 Methodology

2.1 WRF-LES model

The atmospheric boundary layer is simulated by the compressible non-hydrostatic WRF-LES model, which computes large energy-containing eddies at the resolved scale directly and parameterises the effect of small unresolved eddies on the resolved field using sub-grid-scale (SGS) turbulence schemes (Moeng et al., 2007). The Favre-filtered equations are as follows (Nottrott et al., 2014; Muñoz-Esparza et al., 2015):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{\mathit{\rho }}}{\partial t}}+{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{\mathit{\rho }}{\stackrel{\mathrm{̃}}{u}}_j}{\partial x_j}}=\mathrm{0},\text{(2)}& {\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i}{\partial t}}+{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i{\stackrel{\mathrm{̃}}{u}}_j}{\partial x_j}}=\mathit{\upsilon }{\displaystyle \frac{{\partial }^{\mathrm{2}}{\stackrel{\mathrm{̃}}{u}}_i}{\partial x_j\partial x_j}}-{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\rho }}}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{p}}{\partial x_i}}-{\displaystyle \frac{\partial {\mathit{\tau }}_{ij}}{\partial x_j}}+{\stackrel{\mathrm{̃}}{F}}_i,\end{array}$$

where i (or j) =1, 2, 3, represents the component of the spatial coordinate, ${\stackrel{\mathrm{̃}}{u}}_i$ is the filtered velocity, x_i is the spatial coordinate, t is the time, $\stackrel{\mathrm{̃}}{p}$ denotes the filtered pressure, $\stackrel{\mathrm{̃}}{\mathit{\rho }}$ is the filtered density, υ is the fluid kinematic viscosity, τ_ij are the SGS stresses, and $\stackrel{\mathrm{̃}}{F_i}$ represents external force terms (normally involving the Coriolis force caused by the rotation of the Earth and the large-scale geostrophic forcing).

For the closure of Eq. (2), τ_ij is parameterised using a SGS model. In this study, the 1.5-order turbulence kinetic energy (TKE) SGS model is used,

$$\begin{array}{}\text{(3)}& {\mathit{\tau }}_{ij}=-\mathrm{2}{\mathit{\upsilon }}_{\mathrm{sgs}}{\stackrel{\mathrm{̃}}{S}}_{ij},\end{array}$$

where ${\stackrel{\mathrm{̃}}{S}}_{ij}$ is the filtered strain-rate tensor and calculated as

$$\begin{array}{}\text{(4)}& {\stackrel{\mathrm{̃}}{S}}_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i}{\partial x_j}}+{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_j}{\partial x_i}}).\end{array}$$

υ_sgs denotes the SGS eddy viscosity and is defined as

$$\begin{array}{}\text{(5)}& {\mathit{\upsilon }}_{\mathrm{sgs}}=C_{k}l{k}_{\mathrm{sgs}}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where C_k is a model constant, and l is the SGS length scale, which equals the grid volume of size (Δ) under neutral conditions (Deardorff, 1970),

$$\begin{array}{}\text{(6)}& \mathrm{\Delta }=(\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z{)}^{\mathrm{1}/\mathrm{3}}.\end{array}$$

k_sgs is the SGS TKE with the transport equation:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\displaystyle \frac{\partial k_{\mathrm{sgs}}}{\partial t}}+{\displaystyle \frac{\partial }{\partial x_i}}\left(k_{\mathrm{sgs}}{\stackrel{\mathrm{̃}}{u}}_i\right)& =-{\displaystyle \frac{{\mathit{\upsilon }}_{\mathrm{sgs}}}{\mathit{Pr}}}{\displaystyle \frac{g}{{\mathit{\theta }}_{\mathrm{0}}}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{\mathit{\theta }}}{\partial z}}+\mathrm{2}{\mathit{\upsilon }}_{\mathrm{sgs}}{\stackrel{\mathrm{̃}}{S}}_{ij}{\stackrel{\mathrm{̃}}{S}}_{ij}\\ & +\left(\mathit{\upsilon }+{\mathit{\upsilon }}_{\mathrm{sgs}}\right){\displaystyle \frac{{\partial }^{\mathrm{2}}{k}_{\mathrm{sgs}}}{\partial x_i\partial x_i}}-C_{\mathit{\epsilon }}{\displaystyle \frac{k_{\mathrm{sgs}}^{\mathrm{1.5}}}{l}},\end{array}\end{array}$$

where $\stackrel{\mathrm{̃}}{\mathit{\theta }}$ is the filtered potential temperature, Pr is the turbulent Prandtl number, and C_ε is a dissipation coefficient (for more details about the parameterisation see Moeng et al., 2007). Without loss of generality, the “$\stackrel{\mathrm{̃}}{}$” notation for all filtered variables is omitted hereafter.

2.2 Synthetic inflow turbulence generator

The synthetic inflow turbulence generator in Xie and Castro (2008) adopted the digital filter-based method and is used in this study. For simplicity, a one-dimensional problem (the streamwise velocity, u, along the x direction) is used as an illustration to describe this method. The two-point velocity correlations R_uu(kΔx) are assumed to be represented by an exponential function:

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\stackrel{\mathrm{‾}}{u_m{u}_{m+k}}}{\stackrel{\mathrm{‾}}{u_m{u}_m}}}=R_{uu}\left(k\mathrm{\Delta }x\right)=\exp \left(-{\displaystyle \frac{\mathit{\pi }\left|k\right|}{\mathrm{2}n}}\right),\end{array}$$

where m, the index that the averaging operator is applied, denotes the mth element of a vector (one-dimensional data series of, for example, the digital-filtered velocity, u, in Eq. 9 below), k is the number of elements for the two-point distance of kΔx, n is related to the integral length scale L=nΔx with the grid size of Δx, and u_m is the digital-filtered velocity,

$$\begin{array}{}\text{(9)}& u_m=\sum _{k=-N}^N{b}_j{r}_{m+k},\end{array}$$

where r_m is a sequence of random data with mean $\stackrel{\mathrm{‾}}{r_m}=\mathrm{0}$ and variance $\stackrel{\mathrm{‾}}{r_m{r}_m}=\mathrm{1}$, N is related to the length scale for the filter (here N≥2n), and b_j is the filter coefficient and can be estimated from

$$\begin{array}{}\text{(10)}& b_k={\stackrel{\mathrm{̃}}{b}}_k/{\left(\sum _{j=-N}^N{\stackrel{\mathrm{̃}}{b}}_j^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where ${\stackrel{\mathrm{̃}}{b}}_k\cong \exp \left(-\frac{\mathit{\pi }\left|k\right|}{n}\right)$. For a two-dimensional filter coefficient, it can be obtained that

$$\begin{array}{}\text{(11)}& b_{jk}=b_j{b}_k,\end{array}$$

which will then be used to filter the two-dimensional random data at each time step,

$$\begin{array}{}\text{(12)}& {\mathit{\phi }}_{\mathit{\beta }}\left(t,x_j,x_k\right)=\sum _{j=-N_j}^{N_j}\sum _{k=-N_k}^{N_k}{b}_{jk}{r}_{m+j,m+k},\end{array}$$

where β indicates the velocity component. At the next time step, the filtered velocity field is calculated as

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\mathrm{\Psi }}_{\mathit{\beta }}\left(t+\mathrm{\Delta }t,x_j,x_k\right)& ={\mathrm{\Psi }}_{\mathit{\beta }}\left(t,x_j,x_k\right)\exp \left(-{\displaystyle \frac{\mathit{\pi }\mathrm{\Delta }t}{\mathrm{2}T}}\right)\\ & +{\mathit{\phi }}_{\mathit{\beta }}\left(t,x_j,x_k\right){\left[\mathrm{1}-\exp \left(-{\displaystyle \frac{\mathit{\pi }\mathrm{\Delta }t}{T}}\right)\right]}^{\mathrm{0.5}},\end{array}\end{array}$$

where T is the Lagrangian timescale representing the persistence of the turbulence, and ${\mathit{\phi }}_m\left(t,x_j,x_k\right)$ is calculated based on Eq. (12). Xie and Castro (2008) demonstrated that Eq. (13) satisfies the correlation functions in an exponential form in space and in time. The two-dimensional filter in Xie and Castro (2008) is more computationally efficient than a three-dimensional filter.

Finally, the velocity field is obtained by using the simplified transformation proposed by Lund et al. (1998),

$$\begin{array}{}\text{(14)}& {\stackrel{\mathrm{̃}}{u}}_i={\stackrel{\mathrm{‾}}{u}}_i+{\mathit{\alpha }}_{i\mathit{\beta }}{\mathrm{\Psi }}_{\mathit{\beta }},\end{array}$$

where

$$\begin{array}{}\text{(15)}& \left[{\mathit{\alpha }}_{i\mathit{\beta }}\right]=\left[\begin{array}{ccc}{\left({\stackrel{\mathrm{̃}}{R}}_{\mathrm{11}}\right)}^{\mathrm{1}/\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ {\stackrel{\mathrm{̃}}{R}}_{\mathrm{21}}/{\mathit{\alpha }}_{\mathrm{11}}& {\left({\stackrel{\mathrm{̃}}{R}}_{\mathrm{22}}-{\left({\mathit{\alpha }}_{\mathrm{21}}\right)}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}& \mathrm{0}\\ {\stackrel{\mathrm{̃}}{R}}_{\mathrm{31}}/{\mathit{\alpha }}_{\mathrm{11}}& ({\stackrel{\mathrm{̃}}{R}}_{\mathrm{32}}-{\mathit{\alpha }}_{\mathrm{21}}{\mathit{\alpha }}_{\mathrm{31}})/{\mathit{\alpha }}_{\mathrm{22}}& {\left({\stackrel{\mathrm{̃}}{R}}_{\mathrm{33}}-{\left({\mathit{\alpha }}_{\mathrm{31}}\right)}^{\mathrm{2}}-{\left({\mathit{\alpha }}_{\mathrm{32}}\right)}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}\end{array}\right],\end{array}$$

and ${\stackrel{\mathrm{̃}}{R}}_{i\mathit{\beta }}$ is the resolved Reynolds stress tensor, which can be estimated based on measurements or other simulations with periodic boundary conditions. The calculations of α_iβ follow an iterative order: α₁₁, α₂₁, α₂₂, α₃₁, α₃₂, and α₃₃.

Figure 1 Integral length scales prescribed at the inlet used in the inflow BASE case (LS1.0).

2.3 Model coupling and configuration

In this study, we firstly configured a WRF-LES model with periodic boundary conditions in both the streamwise and spanwise directions to obtain prior mean profiles of first and second moments of turbulence, such as the vertical profiles of mean velocity and Reynolds stress components, which are required as inputs by the synthetic inflow turbulence generator. Additional essential quantities as inputs of the inflow generator are three integral length scales in the x, y, and z directions, denoted by L_x, L_y, and L_z, respectively (or L_i, $i=x,y,z$). For the inflow BASE case (denoted by LS1.0), the vertical profiles of L_i are specified as functions of z∕H, where H is the boundary layer height (500 m in this study), shown as Fig. 1, similar to those in Xie and Castro (2008). The streamwise length scale (L_x) is specified based on the mean streamwise velocity profile (〈u〉) and a constant Lagrangian timescale T (prescribed in Eq. 13), i.e. L_x=T〈u〉 using Taylor's hypothesis (turbulence is assumed to be frozen while it is moving downstream with a mean speed of 〈u〉). The spanwise length scale (L_y) is specified as a constant value. The vertical length scale (L_z) is specified as a smaller constant value near the bottom and a larger constant value for the upper domain to be closer to the measured length scales, as explained in Xie and Castro (2008) and Veloudis et al. (2007). We conducted a sensitivity study of integral length scales by varying all three baselines L_x, L_y, and L_z with the ratio of 0.6, 0.8, 1.0, 1.2, or 1.4; these individual cases are denoted by “LS0.6”, “LS0.8”, “LS1.0”, “LS1.2”, and “LS1.4”, respectively, in which “LS1.0” is the base case. The size of the computational domain is $\mathrm{9.98}\mathrm{km}\times \mathrm{2.54}\mathrm{km}\times \mathrm{0.5}\mathrm{km}$ (in the x, y, and z directions), with the resolutions of $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{20}\mathrm{m}$ and stretched Δz (from about 3 up to 27 m). The grid number is then $\mathrm{499}\times \mathrm{127}\times \mathrm{49}$. In order to achieve the constant wind direction vertically, the Coriolis force is not activated in this study. The external driving force is specified as a constant pressure gradient force in Eq. (2), similar to that used in Ma and Liu (2017), resulting in a prevailing wind speed of about 10 m s⁻¹ at the domain top. At the top boundary, a rigid lid (“top_lid” in the “namelist.input” file of the WRF-LES model) is specified, and a Rayleigh damping layer of 50 m is used to prevent undesirable reflections (Nottrott et al., 2014; Ma and Liu, 2017) and to maintain a neutral atmospheric boundary layer.

For the cases with the synthetic turbulence at the inlet and periodic conditions in the spanwise direction, the constant pressure gradient force is not necessary anymore. Instead, a pressure-drop between the inlet and outlet is implicitly derived from the prescribed mean momentum profiles as part of the synthetic inflow and the outflow boundary conditions in the solver. The periodic case is used for the validation of the results from the inflow case. The WRF-LES is solved at a time step of 0.2 s. A spin-up period of 6 h is adopted for all inflow cases to allow turbulence inside the domain to reach quasi-equilibrium. The further 1 h outputs with 5 s interval (approximately the advection timescale of the smallest resolved eddies, which is equivalently twice the grid resolution of 20 m) were used for the analysis. We take advantage of the homogeneous turbulence in the spanwise direction (Ghannam et al., 2015) and calculate all resolved-scale turbulent quantities by averaging in the spanwise direction (the y direction) and in time t over the last 1 h period. This averaging is referred to as “the y−t averaging” hereafter and is denoted by 〈φ〉, for example, for the y−t-averaged φ. For a 4D variable, $\mathit{\phi }(t,x,y,z)$, the y−t-averaged φ is a function of xz, i.e. 〈φ〉(x,z); for a variable defined on the x−y plane, e.g. friction velocity $u_{\ast }(t,x,y)$, the y−t averaging u_∗ is a function of x, i.e. 〈u_∗〉(x).

In the synthetic inflow turbulence generator, a uniform mesh is used with resolutions of Δy=20 m (same as that on the physical inlet of the WRF-LES domain) and Δz=4.2 m (slightly larger than the smallest vertical grid spacing of the WRF-LES domain). The three filtered velocity components at the inlet from the inflow generator are then interpolated onto the vertically non-uniform mesh in the WRF-LES domain. It should be noted that the grid resolution can differ between the inflow patch and the inlet of the WRF-LES domain. The standalone synthetic turbulence generator code in Xie and Castro (2008) was originally run on a single processor, whereas the WRF-LES simulation here is run in parallel mode. It is therefore necessary to ensure that each processor in the parallel mode has the same information of the two-dimensional slice of flow field before each processor can extract the corresponding patch from the same two-dimensional inlet data. In this implementation, the synthetic turbulence generator code is firstly run on the master processor at each WRF-LES time step. The generated inlet data are then passed to other processors. The flow field at the inlet of each corresponding processor was then be updated at every time step. The additional computational time for the inflow case is associated with the synthetic inflow turbulence generator and data passing, i.e. non-parallelisation of the current inflow generator. Increasing the integral length scale would increase the computation time since bigger arrays are constructed and computed for the filtered velocity in the synthetic inflow turbulence generator, as in Eq. (9) for the larger integral length scale.

Figure 2 Horizontal slice of instantaneous streamwise velocity component, u (m s⁻¹), at $z/H=\mathrm{0.1}$ after a simulation time of 6 h: (a) the fully periodic case, (b) the synthetic inflow BASE case (LS1.0), and (c) the inflow case without perturbations at the inlet.

3 Results

3.1 BASE case output

3.1.1 Horizontal slices of instantaneous streamwise velocity component

Figure 2 illustrates the horizontal slices of the instantaneous streamwise velocity component at $z/H=\mathrm{0.1}$ in the periodic case, the synthetic inflow case (LS1.0 in Fig. 1a), and the inflow case without inlet perturbations (with mean information only) after a simulation time of 6 h. The synthetic turbulence structures imposed at the inlet are advected into the domain and are adjusted by the model dynamics at further downwind distances. After an adjustment distance (about $x/H=\mathrm{5}$–10), the inflow case (LS1.0) clearly generates turbulence streaks, which are similar to these in the periodic case. Other quantities that may further demonstrate this adjustment distance will be discussed in the following subsections. This suggests that the synthetic turbulence generated at the inlet can develop into realistic turbulence with well-configured structures from an adjustment distance downwind of about $x/H=\mathrm{5}$–10. For the inflow case without inlet velocity perturbations, there is almost no turbulence generated in the domain even after several hours of simulation. This is consistent with other similar tests using engineering CFD codes with no synthetic turbulence added at the inlet, e.g. Xie and Castro (2008), which confirms that a very long distance (e.g. 100 times the boundary layer thickness) is needed to allow turbulence to develop. This indicates the importance of imposing synthetic turbulence, or at least some form of random perturbations (e.g. Muñoz-Esparza et al., 2015), at the inlet. The inflow case without the inlet velocity perturbations is not presented in later sections.

Figure 3 Spatial variation of $\langle u_{\ast }\rangle /u_{\ast }$ for the periodic case and the inflow case (LS1.0), where 〈u_∗〉 is the y−t-averaged local friction velocity and u_∗ is the $x-y-t$-averaged friction velocity for the periodic case.

3.1.2 Development of local friction velocity

Figure 3 shows the development of the y−t-averaged local friction velocity, 〈u_∗〉(x), for the periodic case and the inflow BASE case (LS1.0), normalised by u_∗, the $x-y-t$-averaged friction velocity for the periodic case. The variation of the local friction velocity is within ±0.5 % of u_∗ along the streamwise direction for the periodic case and is within 1.5 % of u_∗ for the inflow case after a downwind distance of $x/H=$7. There is a larger variation close to the inlet region ($x/H<\mathrm{7}$) for the inflow case. This is because the imposed turbulence on the inflow plane is “synthetic”, of which only the first-order and second-order moments, integral length scales, and the spectra aim to match the prescribed data (Bercin et al., 2018). It must develop over a certain distance in the WRF-LES domain before it can be fully developed “realistic” turbulence.

Figure 4 Horizontal profiles (spatially and temporally averaged) of (a) $\langle u\rangle /u_{\ast }$, (b) $\langle u^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$, (c) $\langle v^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$, (d) $\langle w^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$, (e) $\langle u^{\prime }{w}^{\prime }\rangle /u_{\ast }^{\mathrm{2}}$, and (f) $\langle \text{TKE}\rangle /u_{\ast }^{\mathrm{2}}$ at $z/H=\mathrm{0.1}$ and $z/H=\mathrm{0.5}$ in the periodic case and the inflow case (LS1.0).

3.1.3 Horizontal profiles of mean flow and turbulence quantities

Figure 4 illustrates the y−t-averaged horizontal profiles of the normalised mean streamwise velocity component, normal and shear turbulent stresses, and TKE at $z/H=\mathrm{0.1}$ and $z/H=\mathrm{0.5}$ for the periodic case and the inflow case (LS1.0), respectively. These horizontal profiles show the development of synthetic turbulence along the streamwise direction. There are only slight differences in the normalised mean streamwise velocity component ($\langle u\rangle /u_{\ast })$ between the periodic case and the inflow case. This suggests that the inflow case reproduces successfully the desired mean wind. The curves of normalised streamwise velocity variance ($\langle u^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$) for both cases match well with each other downstream from $x/H=\mathrm{7}$–8, although there is a sudden jump close to the inlet and a subsequent decrease until the location of convergence. The horizontal profiles of normalised cross-stream velocity variance ($\langle v^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$) for the inflow case are in a good agreement after a developing distance of $x/H=\mathrm{10}$–12, compared with those for the periodic case. The development convergence of normalised vertical velocity variance ($\langle w^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$) is achieved after a distance of about $x/H=\mathrm{5}$–10 from the inlet. The development distance of turbulent shear stress ($\langle u^{\prime }{w}^{\prime }\rangle /u_{\ast }^{\mathrm{2}}$) is about $x/H=\mathrm{5}$–15. Since the streamwise velocity variance comprises a large proportion of TKE, the development distance for TKE is similar to that for the streamwise velocity variance, i.e. about $x/H=\mathrm{7}$–8. The distance needed for different quantities to reach a converged state differs from each other, and it is about $x/H=\mathrm{5}$–15.

Figure 5 Spatially and temporally averaged vertical profiles of (a) $\langle u\rangle /u_{\ast }$, (b) $\langle u^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$, (c) $\langle v^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$, (d) $\langle w^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$, (e) $\langle u^{\prime }{w}^{\prime }\rangle /u_{\ast }^{\mathrm{2}}$, and (f) $\langle \text{TKE}\rangle /u_{\ast }^{\mathrm{2}}$ at a series of downwind locations in the inflow case (LS1.0), and the periodic case (also averaged in the streamwise direction).

3.1.4 Vertical profiles of mean flow and turbulence quantities

Figure 5 shows the y−t-averaged vertical profiles of the normalised mean streamwise velocity component, normal and shear turbulent stresses, and TKE at a series of downwind locations, $x/H=\mathrm{0}$, 4, 6, and 10, for the inflow case (LS1.0). Inflow cases are not averaged in the streamwise direction so that the development of turbulence at each downwind location (x∕H) can be investigated. Red lines in Fig. 5 are the spatially (including both in the streamwise and spanwise directions) and temporally averaged vertical profiles for the periodic case. It is noted again that these data for the periodic case are also used as the inputs for a priori turbulence information required by the synthetic inflow turbulence generator. It is also noted that the profiles of the mean velocity and second-order moments at the inlet $(x/H=\mathrm{0})$ are overall in a good agreement with these of the periodic case, which suggests precise settings of the turbulence generator. The profiles of the normalised mean streamwise velocity component ($\langle u\rangle /u_{\ast }$) in the inflow case match closely those of the periodic case. Although the sampled data are limited, this confirms again that the inflow case achieves the desired the mean wind. The normalised streamwise velocity variance ($\langle u^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$) converges towards the periodic profile after $x/H=\mathrm{6}$ as shown in Fig. 5b. Although the vertical profiles of $\langle v^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$, $\langle w^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$ and $\langle \text{TKE}\rangle /u_{\ast }^{\mathrm{2}}$ for the inflow case show small variations between different locations, they are all in a good agreement with the corresponding data of the periodic case. These are consistent with the results shown in Fig. 4. The turbulent shear stress $\langle u^{\prime }{w}^{\prime }\rangle$, which is the cross-correlation between the streamwise and vertical velocity fluctuations, usually converges more slowly than the normal turbulent stresses, e.g. 〈v^′2〉. Overall, the synthetic inflow turbulence generator performs well in terms of the mean flow and the turbulence quantities against the data from the periodic case, as well as the short development distance.

Figure 6 Spectra of streamwise velocity component for a series of downwind locations at the height of $z/H=\mathrm{0.5}$, k is the angular wavenumber, with 〈u〉 and 〈u^′2〉 the spatially averaged mean and streamwise normal turbulent stress, respectively.

3.1.5 Spectral analysis

Figure 6 illustrates the spectra of the streamwise velocity component at a series of downwind locations ($x/H=\mathrm{0}$, 4, 6, and 10) at $z/H=\mathrm{0.5}$ for the periodic case and the inflow case (LS1.0). For each x location, e.g. $x/H=\mathrm{10}$, the spectrum for the inflow case was first calculated from the streamwise velocity component over a time series of 3600 s with an interval of 5 s for five selected sample locations of y_n ($y/H=\mathrm{1.76}$, 2.16, 2.56, 2.96, and 3.36), namely, $\stackrel{\mathrm{̃}}{u}(t,\mathrm{2}H,y_n,\mathrm{0.5}H)$. The spectral data were then averaged over y_n to give the spectra plotted in Fig. 6.

The spectrum for the periodic case is calculated using the same method as that used for the inflow case, with an additional average over the streamwise direction x. The spectrum at the inlet ($x/H=\mathrm{0}$) possesses the broadest range of wavenumbers where eddies exhibit inertial sub-range behaviour, as evidenced by the wavenumber range within which the slope of each spectrum is approximately $-\mathrm{5}/\mathrm{3}$. There is evidence of the tendency in the profiles from the inlet downstream to recover to that of the periodic case. The spectrum drops slightly at high wavenumbers from the imposed spectra at $x/H=\mathrm{0}$ to downwind locations and approaches the spectrum of the periodic case. The slight drop suggests a decay of small eddies due to the SGS viscosities and the numerical dissipation originating from the advection scheme in the WRF-LES model. The spectra in Muñoz-Esparza et al. (2015) drop at lower wave numbers than those in Fig. 6, mainly due to a coarser resolution (than the current one). Our resolution of 20 m in the horizontal direction is much finer than the resolution of 90 m in Muñoz-Esparza et al. (2015). In other words, the size of the smallest eddy (twice the grid resolution) that can be resolved by the LES model is 40 m in our paper vs. 180 m in Muñoz-Esparza et al. (2015). These confirm that synthetic turbulence with an inertial subrange in the spectrum generated by using the Xie and Castro (2008) method is able to be mostly sustained in WRF-LES for a high resolution. It is noted that for a very high resolution, e.g. of the order of magnitude of 1 m, similar to that used in the simulations of PALM (PALM, 2017; Maronga et al., 2020), the inertial subrange in the spectrum is much wider. It is to be noted that Muñoz-Esparza et al. (2015) also tested the Xie and Castro (2008) method in WRF-LES using the same resolution of 90 m as that for the temperature perturbation method. Again this is rather a coarse resolution to test the performance of the Xie and Castro (2008) method when a spectrum is of interest.

Figure 7 Development of local friction velocity (averaged over spanwise direction and time) with various integral length scales. 〈u_∗〉 is the local friction velocity along the streamwise direction, and u_∗ is the $x-y-t$-averaged friction velocity for the periodic case.

3.2 Sensitivity tests of integral length scale in the flow cases

It is not trivial to obtain “accurate” integral length scales of the inlet turbulence generator. Indeed these data are always incomplete. Therefore, it is necessary to conduct sensitivity tests of the integral length scales. Figure 7 shows the influence of integral length scale on the development of local friction velocity. Various integral length scale ratios (ranging from 0.6 to 1.4) to those (L_x, L_y, and L_z respectively) in the LS1.0 case are tested. Note that these three integral length scales (L_x, L_y, and L_z) are in the same ratio as those respectively in the LS1.0 case. For all inflow cases, there is a sudden change near the inlet due to the imposed “imperfect” inflow turbulence. The adjustment distance to well-established turbulence (i.e. within 4 % error) is generally short, i.e. about $x/H=\mathrm{2}$–7 for the studied cases LS0.6–1.4, but seems shorter for the case with the smaller integral length scales. This suggests that the imposed integral length scales for the inflow turbulence slightly affect the convergence to well-developed turbulence. We conclude that a variation of ±40 % in the integral length scale in the cases LS0.6–1.4 yields a variation of less than 4 % in the local friction velocity after about $x/H=\mathrm{7}$, and that the sensitivity of integral length scale on the local friction velocity is not significant in the WRF-LES model if the used integral length scale is within a reasonable range. This is consistent with that in engineering-type CFD solvers in Xie and Castro (2008).

Figure 8 Horizontal profiles (spatially and temporally averaged) of (a) $\langle u\rangle /u_{\ast }$, (b) $\langle u^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$, (c) $\langle v^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$, (d) $\langle w^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$, (e) $\langle u^{\prime }{w}^{\prime }\rangle /u_{\ast }^{\mathrm{2}}$, and (f) $\langle \text{TKE}\rangle /u_{\ast }^{\mathrm{2}}$ at $z/H=\mathrm{0.5}$ with various integral length scales.

Figure 8 shows the effects of integral length scale on the horizontal profiles of the normalised mean streamwise velocity, normal and shear turbulent stresses, and TKE at $z/H=\mathrm{0.5}$. Figure 8a shows that $\langle u\rangle /u_{\ast }$ is slightly greater for the length scale ratio less than 1.0. This is likely due to a slightly smaller u_∗, which is common for smaller integral length scale cases (as shown in Fig. 7). Figure 8b–d and f show that in general the normal stresses, $\langle u^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$, $\langle v^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$, $\langle w^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$, and $\langle \text{TKE}\rangle /u_{\ast }^{\mathrm{2}}$, increase as the length scale ratio increases. This is because small eddies tend to decay faster than large eddies. It is crucial to note that for those with the integral length scales close to those of LS1.0 (the base case) the development distance to converged turbulence is shorter compared to other cases, indicating that the length scales of the base case are reasonable estimations.

Figure 9 Vertical profiles (spatially and temporally averaged) of (a) $\langle u\rangle /u_{\ast }$, (b) $\langle u^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$, (c) $\langle v^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$, (d) $\langle w^{\prime \mathrm{2}}\rangle /u_{\ast }^{\mathrm{2}}$, (e) $\langle u^{\prime }{w}^{\prime }\rangle /u_{\ast }^{\mathrm{2}}$, and (f) $\langle \text{TKE}\rangle /u_{\ast }^{\mathrm{2}}$ at $x/H=\mathrm{10}$ with various integral length scales.

Figure 9 shows effects of integral length scale on vertical profiles of the mean velocity, normal and shear turbulent stresses, and TKE at a typical streamwise location ($x/H=\mathrm{10}$). These profiles are consistent with those in Fig. 8 and draw the same conclusions as from Fig. 8. For all the tested integral length scales, downstream from $x/H=\mathrm{10}$ both mean and turbulent quantities converge to the periodic case. This suggests again that the mean velocity and the turbulent stresses are not very sensitive to the integral length scales if they are not too different from the realistic values. In general, there are slight differences in $\langle u\rangle /u_{\ast }$ between each case. The magnitudes of turbulent quantities for smaller integral length scales are generally smaller than those for larger integral length scales.

Figure 10 Spectra of streamwise velocity component for a series of downwind locations at $x/H=\mathrm{10}$ and $z/H=\mathrm{0.5}$ with various integral length scales; k is the angular wavenumber, with 〈u〉 and 〈u^′2〉 the spatially averaged mean and streamwise normal turbulent stress, respectively.

Figure 10 shows the effect of the integral length scale on the spectra of the streamwise velocity component at $x/H=\mathrm{10}$ and $z/H=\mathrm{0.5}$. For all cases tested in the current study, the spectra with various integral length scales generally match those of the periodic case at a distance of $x/H=\mathrm{10}$ from the inlet albeit with slight changes of the spectrum for small wavenumber turbulence. A very small variation of the spectra is within margins of uncertainty in the calculation of the spectra from the raw data. All spectra show an inertial subrange of $-\mathrm{5}/\mathrm{3}$ slope, which are consistent with those in the references, such as Xie and Castro (2008), indicating the robustness of the synthetic turbulence generator on the generation of an inertial subrange.

4 Discussion and conclusions

A synthetic inflow turbulence generator (Xie and Castro, 2008) was implemented in an idealised WRF-LES (v3.6.1) model under neutral atmospheric conditions. A WRF-LES model with periodic boundary conditions was firstly configured to provide a priori turbulence statistical data for the synthetic inflow turbulence generator. Previous studies (e.g. Xie and Castro, 2008) suggest that it is important to have an approximation of the integral length scales, which are the key inputs of the inflow turbulence generator. The results from the inflow cases were then compared with those from the periodic case. Sensitivity tests were conducted for the response of the local friction velocity, the mean flow, the Reynolds stresses, and the turbulence spectra for the flow cases for varying integral length scales.

The inflow case with the baseline integral length scales generates similar turbulence structures to those for the periodic case after an adjustment distance of $x/H=\mathrm{5}$–15. The WRF-LES model with the inflow generator reproduces realistic features of turbulence in the neutral atmospheric boundary layer. The development of local friction velocity suggests that a downwind distance of about $x/H=\mathrm{7}$ is required to recover the local friction force for the inflow case, which agrees with the findings in Xie and Castro (2008) and Kim et al. (2013). Keating et al. (2004) suggested a development distance of about 20 times the half-channel depth for modelling a plane channel flow. The difference between this value and our results can be attributed to the different synthetic turbulence generation approaches adopted here versus those adopted by Keating et al. (2004). Laraufie et al. (2011) suggested that an increase in the Reynolds number decreases the adjustment distance when a synthetic inflow turbulence generator is used. For our simulated atmospheric boundary layer flow here, the Reynolds number is extremely large. Thus adopting a synthetic inflow turbulence generator for the atmospheric boundary layer should also be advantageous in engineering applications. Regarding the minimum resolution required to generate turbulence synthetically, our presented results confirm that the tested grid resolution sufficiently resolves the important features.

Horizontal and vertical profiles of mean velocity and second-moment statistics further confirm that a short adjustment distance is required for the development of synthetic turbulence. The mean velocity profiles at all tested locations were in very good agreement with the reference data, while the turbulence second-moment statistics profiles were in reasonable agreement with the reference data about $x/H=\mathrm{5}$–15 downwind of the inlet. An accurate estimation of the second-order moments is crucial for the assessment of the synthetic inflow turbulence generator, in particular when the inflow turbulence information is not completely available. We found varying the integral length scale within ±40 % of the value in the base case has a negligible influence on the mean velocity profiles, while the effects of the variation on the turbulent second-order moment statistics are visible, for example the local friction velocity was within 4 % error of the reference data at $x/H=\mathrm{7}$. The synthetic inflow turbulence generator requires additional computational time compared to periodic boundary conditions. This will be certainly improved by running the synthetic inflow generation subroutine in parallel as a future task. This study is focused on the feasibility of implementing the inflow method (Xie and Castro, 2008) in the meso-to-micro-scale meteorological code WRF and the impact of the key variables (i.e. the integral length scales) on the simulated turbulence development inside the domain. This inflow subroutine has previously been implemented in both serial and parallel mode in several codes, including engineering-type codes Star-CD (Xie and Castro, 2009) and OpenFOAM (Kim and Xie, 2016), and the micro-scale meteorology code PALM (PALM, 2017). Although the current implementation in WRF is affordable for a moderate-sized simulation (e.g. resolution of tens of metres), the technical parallelisation of this inflow subroutine in WRF-LES can be the future work for very large simulation domains with high resolutions.

In summary, the synthetic inflow turbulence generator is implemented successfully into the idealised WRF-LES model. The generated synthetic turbulence is correlated both in space and in time in the exponential form. The spectrum of these data shows an inertial subrange of $-\mathrm{5}/\mathrm{3}$ slope, and this again suggests the capability of the method to generate high Reynolds number flows. The tests on WRF also confirm that this method yields a satisfactory accuracy, after having compared the local friction velocity, the mean velocity, the Reynolds stresses, and the turbulence spectra against the reference data. The WRF-LES model with the synthetic turbulence generator provides promising results as evaluated against the periodic case. The limitation of this method is the requirement of a priori turbulence statistics data and integral length scales, which can be estimated by the similarity theory of the atmospheric boundary layer or experimental data. Sensitivity studies have been performed to address this issue, in particular in terms of the effect of the integral length scale. The implementation of the synthetic inflow turbulence generator (Xie and Castro, 2008) can be extended to the WRF-LES simulation of a horizontally inhomogeneous case with non-repeated surface land-use patterns and be further developed for the multi-scale seamless nesting case from a meso-scale domain with a kilometre-scale resolution down to LES domains with metre-scale resolutions. It is also worthwhile to examine the wind spiral case induced by the Coriolis force in the atmospheric boundary layer.