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 10 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 ‘-5/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 a 15 priori 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, 20 whereas its influence on the second-moment profiles is more visible, in particular for the smallest integral length scales, e.g. with the friction velocity less than 4% error of the reference data at x /H=7. This implementation enables a WRF-LES simulation of a horizontally inhomogeneous case with non-repeated surface landuse patterns and can be extended so as to conduct a multi-scale seamless nesting simulation from a meso-scale domain with a km-resolution down to LES domains with metre resolutions.


Introduction
Atmospheric boundary layer flow involves a wide range of scales of eddies, from quasi two-dimensional structures at the mesoscale scales to three-dimensional turbulence (normally with higher Reynolds number, i.e. Re ~10 8 -10 9 ) at the microscale 30 (Munoz-Esparza et al., 2015). The Weather Research and Forecasting (WRF) model (Skamarock and Klemp, 2008) provides the capability 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 km to 1-10 km using a fully compressible and non-hydrostatic Reynoldsaveraged Navier-Stokes (RANS) solver (Skamarock and Klemp, 2008), which captures the behaviours 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. There still remains a challenge for downscaling from a mesoscale simulation (resolutions down to 1 km, capturing mean information only) to an LES scale (tens 5 of meters or below, capturing additional turbulence information) (Doubrawa et al., 2018;Talbot et al., 2012;Chu et al., 2014;Liu et al., 2011), e.g. the appropriate inflow conditions for an LES domain, and the sub-grid scale turbulence schemes suitable for appropriate treatment of the "gray-zone" resolution domain where neither planetary boundary layer (PBL) nor LES parametrisation schemes apply well. Consequently, these two scales of problems are studied separately. Most LES models of atmospheric boundary layer flow at the microscale use periodic boundary conditions and simplified large-scale geostrophic 10 forcing for idealised simulations. However, implicit in the use of periodic boundary conditions is the assumption that atmospheric fields and the underlying landuse have repeated periodic features. This assumption may be unrealistic for real landscapes where landuse 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 outmost domain for the nested 15 cases (e.g. Moeng et al., 2007;Khani and Porte-Agel, 2017;Nunalee et al., 2014). 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. 20 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/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 25 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 a downstream boundary back to the upstream inlet. 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 30 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), digitalfilter-based method (e.g. Xie and Castro, 2008;Klein et al., 2003;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 be well-established. For more information about the above synthetic turbulence generation methods, readers are recommended to read Tabor and Baba-Ahmadi (2010), Wu (2017) and Bercin et al. (2018). 5 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. Munoz-Esparza et al. (2014) 10 extended the perturbation method of Mirocha et al. (2014) and proposed four methods, i.e. 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 15 development of turbulence. Munoz-Esparza et al. (2015) further generalised the cell perturbation method of Munoz-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 20 conditions. Generally speaking, these methods impose "white-noise" perturbations, thus having a flat spectrum, to a variable (e.g. temperature) at the inlet, and the model dynamics will "process" the signals once these signals are advected into the domain, e.g. to dissipate high-wavenumber signals quickly and to adjust low-wavenumber signals gradually. These methods are not the classic inflow turbulence generation methods, which are aimed at providing spatially and temporally correlated wind fields with appropriate power spectra. It is thus not surprising that a large distance of about 20-40 boundary-layer depths 25 (Munoz-Esparza et al., 2015;Mazzaro et al., 2019) is normally required to allow a transition to fully-developed turbulence.
The optimisation and generalisation of these methods would also require intensive testing. Munoz-Esparza et al. (2014) commented that 'the use of temperature perturbations presents an alternative'. Munoz-Esparza and Kosovic (2018) extended the cell perturbation method of the inflow turbulence generation to non-neutral atmospheric boundary layers.

30
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., 2019). 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, Section 2 describes the methodology of WRF-LES model and the technique of the synthetic inflow turbulence generator; Section 3 presents the results of the WRF-LES model with the use of the synthetic inflow turbulence generator; and Section 4 states the conclusions and future work.

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 subgrid-scale (SGS) turbulence schemes (Moeng et al., 2007). The Favre-filtered equations are (Nottrott et al., 2014;Munoz-Esparza et al., 2015): 10 where i (or j) = 1, 2, 3, represents the component of the spatial coordinate, ̃ is the filtered velocity, is the spatial coordinate, t is the time, ̃ denotes the filtered pressure, ̃ is the filtered density, is the fluid kinematic viscosity, are the SGS stresses, ̃ represents external force terms (normally involving the Coriolis force caused by the rotation of the Earth and the large-scale 15 geostrophic forcing).
For the closure of Eq. (2), are parameterised using a SGS model. In this study, the 1.5-order turbulent kinetic energy (TKE) SGS model is used, where ̃ is the filtered strain-rate tensor and calculated as, denotes the SGS eddy-viscosity and is defined as, where is a model constant, ℓ is the SGS length scale and under neutral conditions, ℓ equals the grid volume of size (Δ) (Deardorff, 1970), is the SGS TKE with the transport equation where ̃ is the filtered potential temperature, Pr is the turbulent Prandtl number, is dissipation coefficient (for more details about the parameterisation see Moeng et al. (2007)). Without loss of generality, the "̃" notation for all filtered variables is omitted hereafter.

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 10 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 ( Δ ) are assumed to be represented by an exponential function: where , the index that the averaging operator is applied, denotes the m-th element of a vector (one-dimensional data series 15 of, for example, the digital-filtered velocity, u, in (9) below), is the number of elements for the two-point distance of Δ , is related to the integral length scale = Δ with the grid size of Δ , is the digital-filtered velocity, where is a sequence of random data with mean ̅̅̅ = 0 and variance ̅̅̅̅̅̅ = 1, is related to the length scale for the filter (here ≥ 2 ), and is the filter coefficient and can be estimated from 20 For a two-dimensional filter coefficient, it can be obtained that which will then be used to filter the two-dimensional random data at each time step, where indicates the velocity component. At the next time step, the filtered velocity field is calculated as, where is the Lagrangian time scale representing the persistence of the turbulence, φ ( , , ) 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), 5 where and ̃ 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: 11, 21, 22, 31, 32, and 33. 10

Model coupling and configuration
In this study, we firstly configured a WRF-LES model with periodic boundary conditions in both streamwise and spanwise directions to obtain a priori 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 input by the synthetic inflow turbulence generator. Additional essential quantities as input of the inflow generator are three integral length scales in the x, y and z directions, denoted by , 15 and , respectively (or , i=x,y,z). For the inflow BASE case (denoted by LS1.0), the vertical profiles of are specified as functions of / , where 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 ( ) is specified based on the mean streamwise velocity profile (〈 〉) and a constant Lagrangian time scale T (prescribed in Eq. 13), i.e. = 〈 〉 using Taylor's hypothesis (turbulence is assumed to be frozen while it is moving downstream with a mean speed of 〈 〉). The spanwise length scale ( ) is specified as a constant value. The 20 vertical length scale ( ) 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 baseline , with a same 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", "LS1.4", respectively, in which "LS1.0" is the base case. The size of the computational domain is 9.98 km×2.54 km×0.5 km (in x, y and z directions), with the 25 resolutions of Δ = Δ = 20 and stretched Δ (from about 3 m up to 27 m). The grid number is then 499×127×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 -1 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 5 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 second interval (~ the advection timescale of the smallest resolved eddies, which is equivalently twice the grid resolution of 20 m) were used for the analysis. In this study, by taking advantage of the homogeneous turbulence 10 in the spanwise direction (Ghannam et al., 2015), we calculate all resolved-scale turbulent quantities by averaging in the spanwise (the y-direction) 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, ( , , , ), the y-t averaged is a function of , , i.e. 〈 〉( , ); for a variable defined on the x-y plane, e.g. friction velocity * ( , , ), the y-t averaging * is a function of x, i.e. 〈 * 〉( ). 15 In the synthetic inflow turbulence generator, a uniform mesh is used with resolutions of Δ = 20 (same as that on the physical inlet of the WRF-LES domain) and Δ = 4.2 (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) 20 wass originally run on a single processor, whereas the WRF-LES simulation is run in parallel mode. It is therefore necessary to ensure that each processor in the parallel mode has the same information of the 2-dimensional slice of flow field before each processor can extract the corresponding patch from the same 2-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 25 accordingly. 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.  Fig. 1a), and the inflow case without inlet perturbations (with mean information only) after 5 6 hours' simulation time. The synthetic turbulence structures imposed at the inlet are advected in the domain, and are adjusted by the model dynamics at further downwind distances. After an adjustment distance (about / = 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 inflow turbulence generator can generate realistic well-configured turbulence structures from an adjustment distance downwind of 10 about / = 5-10. For the inflow case without inlet velocity perturbations, there is nearly 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 demonstrated that a very long distance (e.g 100 times 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. Munoz-Esparza et al., 2015) at the inlet. The inflow case 15 without inlet velocity perturbations is not presented in the later sections. Figure 3 shows the development of the y-t averaged local friction velocity, 〈 * 〉( ), for the periodic case and the inflow BASE case (LS1.0), normalised by * , the x-y-t-averaged friction velocity for the periodic case. The variation of the local friction velocity is within ±0.5% * along the streamwise direction for the periodic case and is slightly higher (within 1.5% * ) than 20 that for the inflow case after a downwind distance of / =7. There is a larger variation close to the inlet region ( / < 7)

Development of local friction velocity
for the inflow case. This is because the imposed turbulence on the inflow plane is 'synthetic', which develops in a certain distance in the WRF-LES domain. suggests that the inflow case reproduces successfully the desired mean wind profile. The curves of normalised streamwise velocity variance (〈 ′ 2 〉/ * 2 ) for both cases match well with each other from / = 7-8, although there is a sudden jump close 30 to the inlet and a subsequent decrease until the location of convergence. The horizontal profiles of normalised cross-stream velocity variance (〈 ′ 2 〉/ * 2 ) for the inflow case are in a good agreement after a developing distance of / = 10-12 , compared with those for the periodic case. The development of normalised vertical velocity variance (〈 ′ 2 〉/ * 2 ) is achieved after a developing distance of about / = 5-10. The length scale of the development of shear turbulent stress (〈 ′ ′〉/ * 2 ) is about / = 5 − 15. Since the streamwise velocity variance has a major contribution to TKE, the developing distance for 5 TKE is similar to that for the streamwise velocity variance, i.e. about / = 7-8. The distance needed for different quantities to reach a converged state differs from each other, and it is about / = 5-15. 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, / = 0, 4, 6 and 10, for the inflow case (LS1.0). Inflow cases 10 are not averaged in the streamwise direction so that the development of turbulence at each downwind location ( / ) 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. The normalised mean streamwise velocity component (〈 〉/ * ) profiles of the inflow case match closely those of the periodic case. Although the sampled data 15 are limited, this suggests that the inflow case achieves the desired the mean wind profile.. It is noted that the profiles of the mean velocity and second order moments at the inlet ( / = 0) are overall in a good agreement with these of the periodic case, which further suggests a satisfactory performance of the turbulence generator. The normalised streamwise velocity variance (〈 ′ 2 〉/ * 2 ) converges towards the periodic profile after / = 6 as shown in Fig. 5 (b). Although the vertical profiles of 〈 ′ 2 〉/ * 2 , 〈 ′ 2 〉/ * 2 and 〈 〉/ * 2 for the inflow case show small variations between different locations, they are all in a 20 good agreement with the data of the periodic case. These are consistent with the results shown in Fig. 4. The Reynolds shear stress 〈 ′ ′〉, which is the cross correlation between the streamwise and vertical velocity fluctuations, usually converges slower than the normal Reynolds stresses, e.g. 〈 ′ 2 〉. Overall, the synthetic inflow turbulence generator performs well in terms of the development of the mean flow and the turbulence quantities against the data from the periodic case. Figure 6 illustrates the spectra of the streamwise wind component at a series of downwind locations ( / = 0, 4, 6, and 10) at / = 0.5 for the periodic case and the inflow case (LS1.0). For each x-location, e.g. / = 10, the spectrum for the inflow case was first calculated from the streamwise wind velocity component over a time series of 3600 s with an interval of 5 s for five selected sample locations of ( / = 1.76, 2.16, 2.56, 2.96 and 3.36), namely, ̃( , 2 , , 0.5 ). The spectral data were then averaged over to give the spectra plotted in Fig. 6. 30

Spectral analysis 25
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 . It is shown that the spectrum at the inlet (x/H=0) possesses the most broad range of the -5/3 slope compared to the others. There is an 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 / = 0 to downwind locations, and to recover towards the spectrum of the periodic case. The slight drop suggests a decay of small eddies 5 due to the SGS and molecular viscosities. The spectrum in Munoz-Esparza et al. (2015) drops steeper at high wavenumbers, mainly due to a coarser resolution (noticing that their plots were for with the inertial subrange of -2/3 slope). Our spectrum for has a broad range of the inertial subrange of -5/3 slope, indicated in Fig. 6. This is partially attributed to the fact that our resolution of 20 m in the horizontal direction is much finer than their resolution of 90 m. 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 Munoz-10 Esparza et al. (2015). Munoz-Esparza et al. (2015) also compared the stochastic perturbation method with those obtained using Xie and Castro (2008). These confirm that synthetic turbulence with an inertial subrange in the spectrum generated by using 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. in the order of magnitude 1 meter, similar as that used in the simulations of PALM (PALM, 2017), the inertial subrange in the spectrum is wider. 15

Sensitivity tests of integral length scale in the flow cases
It is not trivial to optimise the integral length scales of the inlet turbulence generator. 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. Length scale ratios from 0.6 to 1.4 to those ( , and ) in the LS1.0 case are tested. For all inflow cases, there is a sudden change near the inlet due to the imposed inflow turbulence. The adjustment distance to well-established 20 turbulence is generally shorter for the case with the smaller integral length scale, i.e. about / = 2-7 for the cases LS0.6-1.4. This suggests that the imposed integral length scales for the inflow turbulence affect slightly the convergence to welldeveloped turbulence. It is also observed that a variation of ±40% in the integral length scales in the cases LS0.6-1.4 yields a variation of less than 3% in the local friction velocity after about / = 5 . This suggests that the sensitivity of the tested integral length scales on the local friction velocity is not significant in the WRF-LES model, which is consistent with that in 25 engineering type CFD solvers in Xie and Castro (2008). Once the inflow turbulence is established (e.g. after / = 10), the local friction velocity is slightly greater for larger integral length scales. Figure 8 shows the effects of integral length scale on the horizontal profiles of the normalised mean streamwise velocity component, normal and shear turbulent stresses, and TKE at / = 0.5. Figure 8 (a) shows that 〈 〉/ * is slightly greater for 30 the length scale ratio less than 1.0. This is likely due to a slightly smaller * , which is common for smaller integral length scale cases (as shown in Fig. 7). Figures 8 (b-d) and (f) show that in general the normal stresses, 〈 ′ 2 〉/ * 2 , 〈 ′ 2 〉/ * 2 , 〈 ′ 2 〉/ * 2 , and 〈 〉/ * 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 the 'accurate' (compared with the periodic case) one, i.e.
the LS 1.0 case, the development distance to converged turbulence is shorter compared to other cases. Figure 9 shows effects of integral length scale at a typical streamwise location ( / = 10) on vertical profiles of the mean 5 velocity, normal and shear turbulent stresses and TKE. These profiles support the conclusions drawn from Fig. 8. For / = 10, both mean and turbulent quantities converge approximately to the periodic case. In general, there are slight differences in 〈 〉/ * between each case. The magnitudes of turbulent quantities for smaller integral length scales are generally smaller than those for larger integral length scales. This suggests that the mean velocity and the turbulent Reynolds stresses are not very sensitive to the integral length scales if they are not too different from the realistic values. 10 Figure 10 shows the effect of the integral length scale on the spectra of the streamwise velocity component at / = 10 and / = 0.5. For all cases in the current study, the spectra with various integral length scales generally match those of the periodic case at a developing distance of / = 10, albeit with slight changes of the spectrum for small wavenumber turbulence. A very small variation of the spectra is within the uncertainty of the calculation of spectrum from the raw data. 15 The spectra show an inertial subrange of -5/3 slope, which are consistent as those in the references, such as Xie and Castro (2008).

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 20 a priori turbulence statistical data for the synthetic inflow turbulence generator. The integral length scales were estimated at appropriate ratios to the boundary layer height as in (Xie and Castro, 2008). The results from the inflow cases were then compared with those from the periodic case. It is important to estimate the integral length scales, which are the key inputs of the inflow turbulence generator. Therefore 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. 25 The inflow case with the baseline integral length scales generates similar turbulence structures to those for the periodic case after an adjustment distance of / = 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 / = 7 is required to recover the local friction force for the inflow case, which is consistent with in the 30 findings of Xie and Castro (2008) and Kim et al (2013). Keating et al. (2004) suggested a development distance of about 20 half-channel depth for modelling a plane channel flow. The difference between this value and our results may be owe to the different synthetic turbulence generation approaches Keating et al. (2004) adopted. 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 case of the atmospheric boundary layer here, the Reynolds number is extremely large. Thus adopting 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 5 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 in the domain were close to the desired profiles, while the turbulence second moment statistics profiles were in reasonable agreement with the 10 desired profiles about / = 5-15 downwind of the inlet. The adjustment distances of second moments are crucial for the assessment of the synthetic inflow turbulence generator. Reducing the integral length scales can shorten the adjustment distance. We found varying the integral length scale does not materially influence the mean velocity profiles, but affects the turbulence second moment statistics more noticeably. The synthetic inflow turbulence generator requires additional computational time compared to periodic boundary conditions. This will be certainly improved by running the synthetic inflow 15 generation subroutine in parallel as a future task. This study is focused on the feasibility of implementing the inflow method (Xie & 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 of codes Star-CD (Xie and Castro, 2009) and OpenFOAM (Kim and Xie, 2016), and the micro-scale meteorology code PALM (PALM, 2017). Although the 20 current implementation in WRF is affordable for a moderate-sized simulation (e.g. tens of meters resolutions), the technical parallelisation of this inflow subroutine in WRF-LES can be the future work for very large simulation domains with high resolutions. Users of our open source subroutine may offer this technical contribution.
In summary, the synthetic inflow turbulence generator is implemented successfully into the idealised WRF-LES model. The 25 generated two-dimensional slices of data are correlated both in space and in time in the exponential form. The spectrum of these data shows a broad inertial subrange of -5/3 slope, and this again suggests the capability of the method to generate high Reynolds number flows. These 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 30 periodic case. The limitation of this method is the requirement of a priori turbulence statistic 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 effect of the integral length scale. We conclude that within a certain range of the integral length scale, the numerical results are not significantly sensitive. 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 km-resolution down to LES domains with metre resolutions.