Large eddy simulation of the neutral atmospheric boundary layer: performance evaluation of three inflow methods for terrains with different roughness

In large eddy simulations (LES) of the atmospheric boundary layer (ABL), the choice of inflow method and the specification of surface roughness are crucial for obtaining accurate results. The consistency between inflow conditions and wall boundary condition has been investigated in depth for the RANS approach, whereas it needs further analysis for LES. This paper investigates the combined effect of inflow method and aerodynamic roughness length for LES of the neutral ABL. Three basic inflow generators in combination with three terrain types are evaluated in terms of streamwise homogeneity of the vertical profiles of mean velocity and turbulence kinetic energy. The results show that the precursor method best preserves the homogeneity of the profiles for all roughness lengths considered, with mean absolute deviations up to 0.7% for mean velocity and 3.5% for turbulence kinetic energy at the outlet of the computational domain. Among the synthetic methods, the Vortex Method performs satisfactorily for mean velocity and, to a lesser extent, turbulence kinetic energy, whereas the random flow generation (RFG) method leads to the largest deviations from the target profiles. Finally, the value of the aerodynamic roughness length is found to only weakly influence the performance of the inflow methods

In the past much attention has been paid to the accurate CFD modeling of the atmospheric boundary layer (ABL), both using Reynoldsaveraged Navier-Stokes (RANS) and LES approaches, with particular regard to the correct implementation of boundary conditions. Specifically, inconsistencies between inlet conditions and wall treatment for RANS modeling of the ABL have been investigated by many authors (Blocken et al., 2007a(Blocken et al., , 2007bGorl e et al., 2009;Hargreaves and Wright, 2007;Parente et al., 2011;Richards and Hoxey, 1993). Such inconsistencies can lead to unintended streamwise gradients in the vertical profiles of mean wind speed and turbulence quantities. As a consequence the incident profiles, defined as the vertical profiles of mean velocity and turbulence quantities detected in an empty computational domain at the position where the obstacle(s) would be placed in a non-empty domain, might differ from those prescribed at the inlet boundary. Since these incident profiles may have a considerable impact on the results of numerical simulations of flow around buildings or other obstacles, it is nowadays considered best practice to verify the homogeneity of the vertical profiles of the ABL in an empty computational domain (Blocken et al., 2007a(Blocken et al., , 2007bBlocken, 2015;Franke et al., 2007) before running the main simulation including the objects of interest.
In case LES is chosen to resolve the flow, a time-varying velocity field has to be specified at the inflow boundary using an appropriate inflow generator; in addition, a suitable model that accounts for surface roughness needs to be employed.
With respect to the inflow methods, the generated unsteady flow should faithfully represent the mean wind speed and turbulence characteristics of the flow field (Tabor and Baba-Ahmadi, 2010;Tamura, 2008;Wu, 2017). Comprehensive reviews of the topic exist in the literature (Keating et al., 2004;Tabor and Baba-Ahmadi, 2010;Wu, 2017). According to Tabor and Baba-Ahmadi (2010) the inflow methods can be classified in two main families: precursor simulation methods and synthetic methods.
With regard to the roughness modeling, in LES three categories of methods exist (Rodi et al., 2013). One approach is to explicitly model the roughness elements, whose walls are treated as no-slip walls in the CFD code. An alternative approach consists of applying momentum forcing terms, which is commonly used to model the effect of vegetation on the flow field. The third approach takes into account the ground roughness through prescription of shear stresses at the wall.
Although comparative studies have been performed in the past in order to assess the separate impact of inflow turbulence generation methods (Sections 2.1-2.3) and aerodynamic roughness length (Section 2.4), to the best of our knowledge no comparative studies that take into account changes of both inflow methods and aerodynamic roughness length have yet been performed. Hence, the main focus and novelty of this paper concerns the combined effect of different inflow generators and different terrain roughness lengths.
In the simulations reported in this manuscript the performance of a precursor method is compared to that of two basic but well-known synthetic methods, i.e. the Vortex Method (Mathey et al., 2003(Mathey et al., , 2006Sergent, 2002) and the RFG method (Kraichnan, 1970;Smirnov et al., 2001), for the LES of the neutral atmospheric boundary layer in an empty domain above three different terrain types: rural, suburban and urban. The above-mentioned synthetic methods are selected for the present study since they are frequently employed for comparative studies (e.g. Bazdidi-Tehrani et al., 2016;García et al., 2015;Poletto et al., 2013;Yan and Li, 2015) and in recent studies concerning ABL flows (e.g. Gousseau et al., 2013;Iousef et al., 2017). In order to model the terrain roughness, a prescribed wall shear stress is applied to the wall, since this approach is commonly adopted for ABL flows (Rodi et al., 2013). In particular, the model proposed by Thomas and Williams (1999) is used due to its effectiveness in modeling roughness effects and ease of implementation. In order to assess the effect of the aerodynamic roughness length on the performance of the inflow methods, the streamwise homogeneity of the vertical profiles of mean velocity and turbulence kinetic energy is evaluated.
In the following section, an overview of inflow methods and roughness treatments is provided. In Section 3, computational domain and numerical settings of the CFD simulations are described, including the chosen inflow methods and wall shear stress boundary condition. In Section 4, the results of the numerical simulations are presented: first a quality check of the flow in the precursor domain is performed, then the profile homogeneity of mean velocity and turbulence kinetic energy is discussed along with an analysis of turbulence statistics. Finally, the main points of this work are discussed and summarized in Section 5.

Precursor methods
An effective method to obtain inflow data is to simulate a precursor domain with periodic boundary conditions driven by a forcing term (e.g. a pressure gradient) and, once the flow is fully developed, sample velocity data in a cross-sectional plane and store them in a database. In this way a library is created that can be used to provide inflow boundary conditions for the main simulation domain (successor domain). Note that in this way a periodicity is introduced in the successor domain due to the limited extent of the precursor domain. Chung and Sung (1997) compared different procedures to map data from a precursor domain including phase and amplitude jittering (Lee et al., 1992), which are two different ways to manipulate the data from a precomputed library in order to remove the periodicity from the stored velocity signal. However, important limitations were found for all methods considered. Data from previously stored databases can be rescaled (e.g. Schlüter et al., 2004) in order to match prescribed statistical properties. Lund et al. (1998) proposed a method through which the inner and outer layer velocities, sampled in a plane downstream of the inlet, are rescaled according to similarity laws and reintroduced at the inlet of the precursor domain. Velocity data are sampled from a cross-sectional plane between the inlet and the recycling plane, and can be either stored to generate a library or directly mapped to the inlet of the successor domain (e.g. Stevens et al., 2014). Nozawa and Tamura (2002) and recently Yang and Meneveau (2016) extended Lund's method for flows over rough surfaces. Another sub-category of precursor methods is the internal mapping (see Tabor and Baba-Ahmadi, 2010), in which precursor and successor domain are combined in a single computational volume. Similar to Lund's method, recycling and rescaling from a sampling plane located downstream of the inlet are applied to generate the new inlet boundary condition; the successor domain is the part of the domain downstream of the sampling plane.

Synthetic methods
The family of synthetic methods can be divided in several groups. The first one is the group of Fourier methods, in which turbulence is generated through a summation of harmonic functions. The random flow generation (RFG) method by Smirnov et al. (2001), who modified the original method by Kraichnan (1970), falls in this category and is divergence-free for homogeneous turbulence and nearly divergence-free for inhomogeneous turbulence. Huang et al. (2010) modified Smirnov's method, so that any assigned energy spectrum could be reproduced; they termed this method the discretizing and synthesizing random flow generation (DSRFG) method. Castro et al. (2011) proposed a further modification of the DSRFG, which was successfully tested in an empty domain and then applied by Li et al. (2015) for LES of flow around a square prism with a focus on pressure and force coefficients. Aboshosha et al. (2015a) also proposed a modification of DSRFG, termed consistent discrete random flow generation (CDRFG), in order to overcome some limitations of the method employed by Huang et al. (2010) and improve reproduction of target spectra and coherence in the flow.
Another relevant group is that of digital filter methods whose aim is to generate velocity fluctuations starting from a set of random data (with zero mean and variance equal to one) using a digital filter based on a correlation function (for instance, Gaussian or exponential). In this method correlations are therefore imposed directly, not through a prescribed energy spectrum (Wu, 2017). Some examples of digital filter methods are described by di Mare et al. (2006), Klein et al. (2003), Veloudis et al. (2007), and Xie and Castro (2008). Kim et al. (2013) proposed a modified version of Xie and Castro's method, which has the important property of being divergence-free. Okaze and Mochida (2017), starting from the work of Xie and Castro (2008), developed a synthetic method based on the Cholesky decompositions of the turbulence flux tensor to generate turbulent fluctuations of both velocity components and a scalar, with prescribed temporal and spatial correlations.
Among the remaining synthetic methods, it is worthwhile to mention the principal orthogonal decomposition (POD) methods, through which inflow conditions are generated from available experimental data (e.g. hotwire or PIV measurements; for further references see Tabor and Baba-Ahmadi, 2010), the Vortex Method and the synthetic eddy method (SEM; Jarrin et al., 2006). The Vortex Method was originally proposed by Sergent (2002) and successively modified (Mathey et al., 2003(Mathey et al., , 2006. With this method fluctuations are generated through a synthesized vorticity field and superimposed on the mean velocity profile. The SEM generates inflow turbulence with prescribed mean velocity, turbulence length scales and Reynolds stresses, through a sum of synthetic eddies that are convected through a virtual box that encloses the inlet plane of the computational domain. Both Vortex Method and SEM are not divergence-free. Nevertheless, a divergence-free modification of the SEM (DFSEM) has been proposed by Poletto et al. (2013).
Since the Vortex Method and the RFG method are the synthetic methods employed in the CFD simulations reported in the present paper, more information about these methods is provided in Section 3.2.2.

Comparative studies
Several comparative studies on the performance of the abovementioned methods have been performed. Keating et al. (2004) presented the results of four inflow generation methods for a turbulent plane channel flow, emphasizing the importance of the spectral content of the generated field at the inlet. Tabor and Baba-Ahmadi (2010) investigated the performance of four inflow methods simulating a channel flow. They concluded that although the precursor methods are more accurate than synthetic ones, the latter have the advantage of being easier to implement and use. Another aspect to consider is that, although no additional calculations are required for synthetic methods, they often require a longer inlet section in which the flow develops towards actual turbulence; this adaptation length may make the computational effort unacceptable and is therefore an important performance indicator. In the work of Yan and Li (2015) the effects of different inflow generators have been evaluated in an investigation of wind loads on a tall building. The four inflow methods investigated were: i) recycling method (Nozawa and Tamura, 2002); ii) Vortex Method; iii) RFG method; iv) DSRFG method. Numerical tests in an empty domain revealed a decay of turbulence intensity for the RFG and recycling methods as the flow travels through the computational domain, for which remedial measures to improve the reproduced turbulence level were proposed, while a better agreement with experiments was obtained with the DSRFG and Vortex Method. Regarding the simulated wind loads the authors reported an adequate performance of the DSRFG method. Lampitella et al. (2012) and Lampitella (2014) compared the performance of a precursor method, the Vortex Method and the RFG method for turbulent duct flow, with square and circular sections. They pointed out the higher accuracy of the precursor method compared to the Vortex Method and, to a larger extent, to the RFG; in addition, they also showed that the Vortex Method requires a development region (in order to match the target profile) of about ten times the hydraulic diameter. Bazdidi-Tehrani et al. (2016) performed a comparative analysis on the prediction of the flow field and contaminant dispersion around an isolated building at wind tunnel scale. In their study four inflow turbulence generation methods were considered: no fluctuations, internal mapping, Vortex Method and RFG. The results showed that using the internal mapping technique leads to a better agreement with experimental data of time-averaged concentrations and concentration fluctuations compared to the other inflow generators. In addition, when internal mapping and the Vortex Method are used, higher values of turbulence intensity are predicted in the wake region compared to the other two methods. García et al. (2015) compared the RFG method to generate fluctuations on top of a uniform mean flow against Kaimal's spectral model (Kaimal et al., 1972) applied to a logarithmic wind profile, and low-turbulence uniform flow for the aerodynamics of a train subjected to cross-wind. They reported large differences in integral aerodynamic coefficients, local pressures on the leeward side of the train and velocity spectra. Poletto et al. (2013) compared DFSEM against other synthetic inflow generators, including the Vortex Method, for turbulent channel flow. In their work it is shown that, when DFSEM is employed, a shorter development region is needed in order to correctly reproduce the target vertical profiles (i.e. those obtained from a periodic LES) of Reynolds stresses, compared to the other methods considered. Tamura and Ono (2003) used a synthetic inflow method (Hoshiya, 1972;Kondo et al., 1997) based on an inverse Fourier transform of a target spectrum (specifically, a Karman-type spectrum) to show that the value of turbulence intensity prescribed at the inflow boundary influences the physical mechanisms of shear layer separation in the case of flow past square and rectangular prisms. Therefore, if streamwise fluctuations are not accurately preserved the turbulence characteristics may differ from the ones prescribed at the inlet, with implications for the accuracy of the numerical simulations.

Roughness modeling
The first method, consisting of explicit modeling of roughness elements (e.g. Nozawa and Tamura, 2002;Yan and Li, 2015;Yang and Meneveau, 2016), requires a fine grid resolution around the roughness elements and therefore a large number of cells. In addition, the geometrical characteristics of the obstacles (height, density and shape) should be properly chosen in order to match the desired aerodynamic roughness length. This information can be derived from wind tunnel experiments, either directly or through empirical correlations (see e.g. Fang and Sill, 1992;Lettau, 1969;Raupach et al., 1980;Wooding et al., 1973).
With respect to the methods belonging to the second category (momentum forcing), some authors (e.g. Kanda and Hino, 1994;Shaw and Schumann, 1992;Watanabe, 2004) used a sink to account for momentum losses due to the presence of vegetation. Scotti (2006) and Stoesser (2010) modeled roughness using ellipsoids and grains (whose mean diameter is prescribed). Here either the momentum equation is corrected taking into account the volume fraction occupied by the roughness elements, or the velocity is forced to zero (applying a source term in the momentum equation) in the cells containing roughness elements. An overview of other techniques belonging to this family of methods is reported in Rodi et al. (2013).
With regard to the methods belonging to the third category (fixed wall shear stress prescribed at the wall), Thomas and Williams (1999) proposed a modified version of the Schumann-Gr€ otzbach (Gr€ otzbach, 1987;Schumann, 1975) boundary condition, which relates the instantaneous shear stress at a certain location of the rough wall to the instantaneous velocity components at the wall-adjacent cell centroid. In their approach the wall shear stress is given by a linear combination of mean and instantaneous velocity. The weight of each part of the velocity is determined by a constant damping coefficient which is problem-dependent. Xie et al. (2004) proposed a modification of the wall model described by Thomas and Williams (1999), where the damping coefficient can be estimated from experimental data or similarity relations. They compared the results obtained using their model with results of other models (Schumann, 1975;Thomas and Williams, 1999). In terms of mean velocity they reported a significant improvement compared with Schumann's model. Little difference was found compared with the model of Thomas and Williams (1999), although some improvement was obtained in the surface layer. Nevertheless, it is important to note that these comparisons were made for only one value of the damping coefficient of the model proposed by Thomas and Williams (1999). A comparison between different approaches is also presented in the work of Lim et al. (2009), where small differences are reported between the different wall shear stress models adopted. Stoll and Port e-Agel (2006) conducted a comparative study between four different wall shear stress models for a neutral ABL: i) Schumann-Gr€ otzbach; ii) shifted Schumann-Gr€ otzbach (Piomelli et al., 1989); iii) local Schumann-Gr€ otzbach (Albertson and Parlange, 1999;Mason and Callen, 1986;Moeng, 1984); iv) Marusic, Kunkel and Port e-Agel's (MKP) method (Marusic et al., 2001). All methods were tested for four aerodynamic roughness lengths, except the last method which was tested for only one value. The results pointed out that increasing ground roughness resulted in increasing variance, skewness and flatness of the surface shear stress, and that the effects on the velocity fluctuations were confined to the surface layer. In addition, the authors reported an unrealistic dependence of the vertical velocity gradients and energy spectra on the roughness length for the first three models, whereas more realistic results were obtained with the MKP model.
In conclusion, in all the studies considered and reported in this section the effects of inflow methods and aerodynamic roughness length are evaluated separately. Conversely, the comparative study presented in this paper tackles these two important aspects together, since they may not be decoupled in the framework of LES, and eventually aims at assessing the combined effect of three inflow generators and three roughness lengths in view of streamwise homogeneity.

Computational domain and grid
All simulations in the present paper are performed at wind tunnel scale, as explained later. The computational domain consists of two sub- The first sub-domain is called the precursor (or auxiliary) domain while the second one is the successor domain. The precursor domain is used to generate: i) the instantaneous velocity profiles, directly introduced in the successor domain when the precursor method is used; ii) the target profiles of mean velocity and turbulence quantities, which are used as input when the synthetic methods are employed to generate inflow turbulence.
The computational grid adopted ( Fig. 2) is the same for both precursor and successor domain. In order to discretize each sub-domain n x Â n y Â n z ¼ 200 Â 75 Â 58 ¼ 870,000 cells are employed; n x , n y and n z represent the number of cells in the x-, y-and z-direction, respectively. Cells are equally spaced along the streamwise (x-) and spanwise (y-) directions. The mesh is made up of cubical cells in the upper part of the domain (z/H > 0.16), while in the lower part of the domain they become rectangular due to grid refinement near the bottom wall (Fig. 2). In addition, a successive cell size ratio of 1.1 is employed in the vertical (z-) direction over the bottom 16% of the domain, with a first cell height of 8 Â 10 À3 m. Here a maximum cell aspect ratio of 4.5 is reached.

Boundary conditions
An overview of the boundary conditions applied at each surface is given in (Fig. 1). For the precursor domain, periodic boundary conditions are employed in the streamwise direction, while zero normal velocity and zero normal gradients of all variables (symmetry in Fig. 1) are applied at the top and lateral boundaries. The latter choice, adopted also by other authors (Bazdidi-Tehrani et al., 2016;Dagnew and Bitsuamlak, 2014;Gousseau et al., 2011), allows for the adoption of larger cell dimensions close to the top and lateral boundaries, compared to the case where a wall boundary condition is adopted. In the second sub-domain, the periodic conditions are replaced by a velocity inlet boundary condition and zero diffusion flux at the outlet (outflow in Fig. 1).

Wall shear stress boundary condition
A rough wall boundary condition is applied in both precursor and successor domain in which stresses are prescribed at the wall along x-and y-direction, τ x and τ y . Specifically, the wall boundary condition proposed by Thomas and Williams (1999) is used, which calculates the stresses as a linear combination of the mean and fluctuating part of the velocity where ρ is the air density, ð〈u〉; 〈v〉Þ and ðu; vÞ are the time-averaged and instantaneous velocities in the x-and y-direction, respectively, u * is the friction velocity, β is the damping coefficient and u log is the velocity at the first grid point evaluated with the logarithmic law where κ is the von K arm an constant, set equal to 0.415, z P is the midpoint of the first cell adjacent to the bottom wall, and z 0 is the aerodynamic roughness length. In Eq.
(1), the mean values ð〈u〉; 〈v〉Þ are obtained by locally averaging in time, not along the bottom wall as done in Thomas and Williams (1999) and in the study of Xie et al. (2004), since the latter will generally not be meaningful in studies employing non-empty domains. In addition, in Thomas and Williams (1999) u log is computed as (2) where a slight modification in the numerator is adopted in order to overcome a limitation concerning the height of the midpoint of the first cell (z P > z 0 ) which would have led to unacceptably high values of z P for the definition employed in Thomas and Williams (1999). Aboshosha et al. (2015b), following the work of Blocken et al. (2007a), commented that the constraint z P > k S,ABL % 30z 0 , where k S,ABL is the equivalent sand-grain roughness height for the ABL, must be respected also for LES when wall functions are employed. Nevertheless, it is worthwhile to clarify that no specifications about the value of z P (in relation to k S,ABL ) are mentioned in the work of Thomas and Williams (1999), since Eq. (1) do not represent a sand-grain-roughness wall function. Furthermore, as pointed out by Blocken et al. (2007a), although the constraint z P > k S,ABL is logical from the physical point of view (flow equations should not be solved inside the solid below the surface roughness), there is no reason why it should not be applied mathematically to enforce consistency. In Eq. (1) the damping coefficient β regulates the relative contribution of the mean and fluctuating part of the velocity to the calculation of the wall shear stresses and influences the magnitude of the RMS values of velocity near the ground. The presence of this damping factor originates from the fact that, for rough surfaces, the contribution of the mean flow to the wall shear stress is higher than the contribution of the fluctuating part of the velocity. Physically this is explained through the observation that velocity fluctuations originating from the presence of roughness elements at the wall may be damped before they interact with the roughness element nearest to the one from which they originate, and therefore contribute less to the wall shear stress calculated by means of Eq. (1) (Xie et al., 2004). If the damping factor β is set equal to 1 the original equation proposed by Schumann (1975) is obtained; in this case the wall shear stresses are proportional to the velocity fluctuations.
An interesting point of discussion is the value to assign to β. In Thomas and Williams (1999) a value of 0.3 is arbitrarily assigned. One option is to tune this parameter, which would then become problem-dependent. In the present work the adopted value for the damping factor is not optimized since it is not an aim of this study to compare different wall models; the same value of 0.3 that was also used by Thomas and Williams is therefore adopted. The wall shear stress boundary condition is implemented through a User Defined Function (UDF) in ANSYS Fluent (ANSYS, 2013), which is the CFD code employed for all calculations in this paper. As a consequence of Eqs. (1) and (2) the midpoint of the first near-wall cell must lie in the logarithmic layer. Correspondingly, the range of values of y þ is 30 < y þ < 500 for all simulations.
At the bottom of the domain different values of the aerodynamic roughness length corresponding to rural, suburban and urban terrain are imposed. As reported in Table 1 three cases are considered, with model scales of 1:292, 1:273 and 1:269, respectively. The values of model scales, friction velocity and terrain roughness used for the neutral ABL simulations correspond to those reported in experimental work by Kozmar (2011).
The flow in the precursor domain is driven by a constant streamwise force (equivalent to a constant streamwise pressure gradient), which is implemented as a source term in the x-momentum equation and whose value is not known a priori. In all the simulations the source term is adjusted using a trial-and-error approach in order to reproduce the theoretical logarithmic mean velocity profile and the friction velocity for the corresponding roughness length as closely as possible.

Inflow turbulence generation methods
The commercial software ANSYS Fluent 15.0 (ANSYS, 2013) is used to perform the calculations reported in the present paper. Three inflow methods are investigated: the precursor method (PM) and two synthetic methods: the Vortex Method (VM; Mathey et al., 2003Mathey et al., , 2006Sergent, 2002) and the RFG method (Kraichnan, 1970;Smirnov et al., 2001),    available in ANSYS Fluent with the name Spectral Synthesizer (SS). For the precursor calculation, a domain with periodic boundary conditions is adopted. The velocity components are instantaneously mapped from a central cross-sectional plane of this domain (map in Fig. 1) to the inlet of the successor domain. The mapping between the precursor domain and the concurrent successor simulation is implemented in ANSYS Fluent through a User Defined Function (UDF), described by Lampitella et al. (2012) and Lampitella (2014). From the precursor simulations the vertical profiles of mean streamwise velocity and turbulence quantities, averaged in the map plane along the spanwise direction, are also used as input for the synthetic methods. Specifically, the profiles of turbulence kinetic energy, turbulence dissipation rate and Reynolds stresses (only the diagonal components for the Vortex Method, diagonal and off-diagonal components for the Spectral Synthesizer) are assigned.
The first synthetic method considered is the Vortex Method. In the VM velocity fluctuations are superimposed on the mean velocity field at the surface where the method is applied. This method is based on the solution of the Lagrangian form of the vorticity equation through a particle discretization. These particles, also called "vortex points" (Mathey et al., 2006) contain information about the vorticity field. The local size of the vortices, σ, is defined using the turbulence kinetic energy, k, and the turbulence dissipation rate, ε, assigned at the inlet boundary In Eq. (3), c ¼ 0:16 (ANSYS, 2013). In order to avoid the size of the vortices being smaller than the resolved scale a lower boundary equal to the grid size is set for σ (ANSYS, 2013). In addition, in ANSYS Fluent the linear kinematic model (LKM; Mathey et al., 2006) is implemented to simulate the effect of bidimensional vortices on the streamwise direction. The calculated velocity fluctuations along the three directions, u 0 i , are then eventually rescaled according to the diagonal Reynolds stresses, u 0 The second synthetic method considered is the Spectral Synthesizer. This method, as remarked in Section 2.2, generates a divergence-free velocity field for a homogeneous turbulence and nearly divergence-free for an inhomogeneous turbulence (Smirnov et al., 2001), taking as input the Reynolds stress tensor u i 0 u j 0 and the length and time scales of turbulence, determined as a function of the k and ε profiles. The prescribed correlation tensor u i 0 u j 0 is diagonalized through an orthogonal transformation tensor, a ij , and the scaling factors c i are obtained. Once a ij , c i and the length and time scales are known, a 3D transient velocity field vðx; tÞ is generated based on Fourier harmonics with random coefficients following a normal distribution. The modeled energy spectrum follows a Gaussian spectral model. Eventually, in order to compute the velocity field which satisfies the prescribed properties (correlation tensor, turbulence length and time scales), vðx; tÞ is rescaled and an orthogonal transformation is applied The resulting velocity fluctuation field is then superimposed on the mean field.

Numerical settings
For pressure-velocity coupling the non-iterative pressure implicit scheme with splitting of operators (PISO; Issa, 1986) is employed. The selected time step size is 2.5 Â 10 À4 s for all simulations, yielding a In addition, it was reported that variations of the C S influence the numerical prediction of the flow field (e.g. Gousseau et al., 2013;Murakami et al., 1987) and that the value of C S can be adjusted empirically (Basu and Port e-Agel, 2006). For the simulations reported in the present work, the value of C S ¼ 0.2 is chosen after preliminary numerical tests where different values of C S were used. The calculations where the PM is used are initialized with results from steady-state RANS calculations employing the realizable k-ε turbulence model (Shih et al., 1995), a logarithmic inlet velocity profile, and turbulence kinetic energy and   (Kozmar, 2011) and computed values of the friction velocity, ðu * Þ theor and ðu * Þ CFD respectively, for the precursor domain (indicated as A, auxiliary) and successor domain (PM, Precursor Method; VM, Vortex Method; SS, Spectral Synthesizer turbulence dissipation rate profiles according to Richards and Hoxey (1993). The calculations employing the two synthetic methods are initialized with an instantaneous velocity field obtained by perturbing a horizontally homogeneous solution given by the spanwise-averaged velocity profile in the map plane, obtained from a previous LES simulation where the PM is employed. As already mentioned above, these profiles (target profiles) are also used as input profiles for the synthetic methods to generate synthetic turbulence. In order to ensure statistical convergence great care is taken that all simulations are run for a sufficiently long time. For rural (case 1), suburban (case 2) and urban (case 3) terrain (Table 1) the simulations are initialized for a period of respectively 18.5, 18.3 and 35.8 large-eddy turnover times (H/u * ), after which statistics are sampled for respectively 36.7, 48.9 and 63.6 large-eddy turnover times. The moving average of the streamwise velocity and the indicator e P conv are used to quantify the convergence. The latter is defined as (e.g. Gousseau et al., 2013 where I is the number of intervals, taken equal to nine, in which the averaging period is divided, P is the monitored point and T avg the averaging time. The values of the indicator and the moving average of the streamwise velocity are monitored in four points: P 1 , P 2 , P 3 , and P 4 , located on a vertical line (y/H ¼ 0) in the map plane with z/H ¼ 0.01, 0.1, 0.3 and 0.5, respectively. In Fig. 3 the moving average for the four points is reported for the case 1. After some initial variation the moving average does not change significantly for the four points, indicating that an acceptable convergence of the simulation is reached. This is in line with the corresponding values of the indicator e P conv , which are initially high (37.5% for the point P 1 ), then they decrease and eventually drop below 2% for the point P 1 and 1% for the points P 2 , P 3 and P 4 .
Eventually, in order to perform statistical analyses, velocities are sampled every 10 time steps, thus with a sampling frequency of 400 Hz. In this way a reasonable amount of data is stored without excessively slowing down the calculations with I/O operations; at the same time, only information about high frequencies (falling in the dissipation range) is overlooked.

Mean velocity and turbulence kinetic energy
Results along horizontal lines for different values of the height z/H (0.02, 0.05, 0.1, 0.3 and 0.7) are presented in Fig. 4 in terms of deviations of the local mean velocity from spanwise-averaged velocity along the line considered.
It is shown that, although the mean velocity is not completely symmetrical in the cross-sectional plane, the mean absolute deviations from the spanwise average (~3-5%) are acceptable for all cases considered. Locally larger deviations are found close to the lateral boundaries of the domain (up to~12-15%) due to the presence of the symmetry boundary condition at the sides, which forces the orthogonal velocity component and the normal gradients of other variables to zero. In preliminary simulations where periodic boundary conditions were applied at the lateral boundaries a large inhomogeneity in the lateral direction was observed due to a lateral instability, although the mean values of the deviation were comparable with the ones presented in Fig. 4. A cyclic instability in the spanwise direction is reported in the work of Thomas and Williams (1999) who, for their case, overcame the problem by applying a filter in the LES code. In the present work the adoption of zero shear stress boundary conditions at the lateral boundaries effectively limited the spanwise instability.
In Fig. 5 the results for four vertical lines in the map plane at y/H ¼ 0, 0.25, 0.5 and 0.75 are presented for non-dimensionalized mean velocity and turbulence kinetic energy. The deviation from the spanwiseaveraged profiles is again limited for both quantities, with larger local differences at the lateral boundaries. For instance, considering the mean streamwise velocity, for the three vertical lines at y/H ¼ 0, 0.25 and 0.5, the mean absolute deviations from spanwise-averaged profiles are lower than 3.5% for all cases, while they increase up to 10.8% (case 3) at y/H ¼ 0.75. In general an acceptable level of spanwise homogeneity is present, especially in the central part of the domain which is usually the region of greatest interest in the numerical simulations.
Finally, the comparison between computed, ðu * Þ CFD , and theoretical, ðu * Þ theor , values of the friction velocity are shown in Table 2. The resulting friction velocity is lower than the theoretical one for all terrain types but is most prominent for case 3. This indicates that, especially for the case with higher roughness (urban terrain), the computed wall shear stresses are underestimated when using the wall boundary condition proposed by Thomas and Williams (1999), and in particular when the Spectral Synthesizer is used due to a possible laminarization of the flow, as explained in Section 4.2.1.

Zero-lag correlations, autocorrelations and integral length scale
An analysis of flow statistics is also reported, with particular regard to probability density functions, skewness and kurtosis as a function of height, zero-lag correlations, autocorrelations with variable time lag, computed according to the definitions reported in ESDU 74030 (1976), and longitudinal integral length scales. Velocity components are sampled along three orthogonal lines, oriented according to the three main directions in the domain, at z/H ¼ 0.23 and y/H ¼ 0.
In Fig. 6 the probability density function (PDF) of the three instantaneous velocity components in the plane z/H ¼ 0.23 is shown at the final simulation time for all the cases. For all velocity components the PDF is satisfactorily close to a normal distribution, which is commonly assumed to be proper for the neutral ABL (ESDU 85020, 2001). In order to quantify the symmetry of the distribution and its flatness compared to the normal distribution, the values of skewness (γ 1 ) and kurtosis (β 2 ) are reported as a function of height in Fig. 7. For a perfectly Gaussian distribution their values are respectively 0 and 3. For all cases there is a similar trend: the values of skewness of all components denote a weakly skewed distribution (jγ 1 j < 1.5), and in particular the deviation from the ideal value is lower for the spanwise velocity than for the other velocity components. As also observed by Lu and Port e-Agel (2010), a low value of the skewness for the spanwise velocity component is necessary for the symmetry of the flow in the spanwise direction. In addition, the negative values of skewness for the streamwise velocity indicate that the prevailing sign of the turbulent fluctuations of this quantity is negative. This is also reported in previous LES studies (Lu and Port e-Agel, 2010;Piomelli, 1993). Finally, the values of kurtosis agree moderately (2 < β < 4) with the ideal value of 3 for all the velocity components at a wide range of heights, likely because of the influence of the symmetry condition applied to the lateral boundaries.
In Fig. 8 zero-lag two-point correlations, ρ ii (x/H) with i ¼ u, v and w are reported. In case 1 and case 3, as expected, the values of ρ quickly drop to zero, becoming negligible in a short distance for all components. Along the streamwise direction this distance is about one-fourth of the domain length for the streamwise velocity, and much less for the other components. In case 2 the correlation for the streamwise velocity along the x-direction reaches a minimum value, higher than zero, which depends on the interval of time samples considered. Eventually, the value of ρ ii (x/H) obtained by considering the whole sampling interval is reported and used for the calculation of the longitudinal integral length scale of turbulence. It is worthwhile to observe that the finite length of the precursor domain is a limitation which is clearly shown by the zero-lag correlations (Fig. 8) and autocorrelations (Fig. 9). As a matter of fact the value of zero-lag correlations has a secondary peak equal to one at x/ H ¼ 4, which is where periodic conditions are applied. In addition, with respect to the autocorrelations, the value of ρ ii (τ) quickly decreases towards zero, then an oscillating behavior with decreasing peaks is observed or all the velocity components. These secondary peaks are a consequence of the periodicity of the domain: the peaks are almost in correspondence with the flow through time, equal to 0.40 s for case 1, 0.43 s for case 2 and 0.47 s for case 3, and they are present for 5-6 flowthrough times before the function is confined to values close to zero.
The longitudinal integral length scale of turbulence, x L u , is calculated and compared with international standards ESDU 85020 (2001), Eurocode (2005) and AIJ guidelines (AIJ, 2004) for the three aerodynamic roughness lengths considered (Table 3). For all the cases there is a large underestimation of x L u compared to the ESDU standard, whereas a moderate agreement is obtained with the values recommended by Eurocode and AIJ for the reference height considered (z/H ¼ 0.23). With respect to the influence of aerodynamic roughness length it can be observed from the predicted values that changes in roughness do not significantly influence the accuracy of prediction of x L u .

Mean velocity and turbulence kinetic energy
Vertical profiles of mean velocity and turbulence kinetic energy are evaluated at different positions in the successor domain and compared to the target profiles, generated through the precursor domain (as reported in Section 3.1), in order to quantify the performance of the inflow methods examined. This comparison is performed for the three aerodynamic roughness lengths considered. All quantities reported are averaged in time and along the spanwise direction. Fig. 10 presents the results of mean streamwise velocity for the crosssectional planes at x/H ¼ 0, x/H ¼ 2 and x/H ¼ 4. A very good homogeneity is observed for the PM, for which the mean absolute deviation from the target profile is very limited compared to the synthetic methods.    The mean absolute deviations from target profiles (Table 4) obtained with the PM at the outlet are 0.2%, 0.3% and 0.7% for cases 1, 2 and 3, respectively. The VM shows satisfactory agreement with prescribed profiles in terms of mean absolute deviation (< 2% at the outlet for all cases), which is nevertheless higher than the PM. Finally, when the SS is used the mean absolute deviations at the outlet section are in the range 6.4-7.5%, depending on the case considered and are therefore considerably higher than for the other methods.  Consequently, the implemented version of this method in ANSYS Fluent turns out to be unable to produce homogeneous profiles of mean streamwise velocity in an empty domain, possibly due to laminarization of the flow (Lampitella, 2014;Lampitella et al., 2012), as confirmed by the lower values of friction velocity (Table 2) and mean velocity in the lower part of the domain compared to the target profiles.
In terms of streamwise homogeneity of turbulence kinetic energy (Fig. 11), the PM more clearly outperforms the other methods for all terrain types considered. The mean absolute deviation from target profiles (Table 5) at the outlet boundary is between 2.8% and 3.5% depending on the roughness length, which is considerably lower than what is obtained with the other methods. When the VM is used, the values of turbulence kinetic energy increase throughout the domain until they come very close to the target profile. The mean absolute deviation along the vertical direction decreases along this development length from 25.0% at the inlet plane (x/H ¼ 0) up to 18.5% at the outlet (x/H ¼ 4) for case 1, from 31.6% to 22.3% for case 2 and from 36.4% to 23.5% for case 3. Hence, the PM displays the least deviation from the target profile compared to the synthetic methods. This can be attributed to the fact that fully developed turbulence is immediately achieved for the PM, whereas a longer distance is needed for the VM. As a matter of fact, for the VM the deviations from the target profile become smaller along the streamwise direction. Finally, the SS shows by far the worst performance among the three methods, with turbulence that is quickly damped by the solver once it enters the domain. According to Lampitella (2014), this is attributable to incorrect spatial correlations of the velocity fluctuations generated at different positions of the inlet plane, possibly caused by an incorrect parallel implementation. The result is a sudden decay of turbulence kinetic energy, with an mean absolute deviation from the target profile which increases up to almost 80% at the outlet for all terrain types considered (Table 5).
In addition to mean absolute deviations from target profiles, in Figs. 12 and 13 the local deviations of mean streamwise velocity and turbulence kinetic energy as a function of height are reported. The aim is to quantify the local impact of the inflow methods, especially in the region of interest for urban simulations. Given an urban area with a tallest building whose height is H b and considering the typical size of the building wake, the region of interest in the vertical direction can tentatively be defined as 2H b . This height (2H b ) corresponds to H/3 for the present domain if best practice guidelines (Franke et al., 2007) are applied, according to which a minimum distance of 5H b has to be observed between the building roof and the upper boundary of the computational domain. In terms of mean velocity (Fig. 12), for all roughness lengths locally smaller deviations from target profiles are obtained when the PM is applied. For the VM the deviations are locally higher in the lower part of the domain (increasing, for instance, up to 13.4% at the outlet for case 3; this value is lower than 9% for the PM). In addition, the SS locally shows a very large decay of mean velocity compared to the other methods. The same trend is observed for turbulence kinetic energy (Fig. 13). For example, if case 1 is considered, the maximum absolute deviations for the PM, VM and SS are respectively 3.9%, 12.4% and 86.5% for z/H < 0.33. It is also interesting to note that in the upper part of the domain, far from the region of interest, the deviations become very large not only for the SS but also for the VM.

Zero-lag correlations, autocorrelations and power spectral density
The zero-lag correlations (Fig. 14) and autocorrelations (Fig. 15) show a significant variation between the three methods. In particular, the curve for the PM is very similar to the one for the precursor domain. This is determined by the periodicity of the velocity field due to the limited extension of the precursor domain, which is the main limitation of the PM. As a consequence, the resulting periodicity of the velocity signal at the inlet of the successor domain yields an increase of ρ ii (x/H) in a region close to the outlet section, and secondary peaks of ρ ii (τ) for lag-times τ, which are multiples of the flow-through time. Differently from the PM, the VM and the SS are not affected by the periodicity since turbulence is not directly injected from the precursor domain but artificially synthesized. For the VM the correlations do not result in an exponential decay, showing instead an oscillating behavior which is almost completely damped at the end of the domain. This aspect, together with the fact that the deviation from target profiles decreases along the streamwise direction, indicates that compared to the PM a longer distance is needed to obtain uncorrelated profiles in space and time, and probably longer still (> 4H) to achieve fully developed turbulence. In contrast, when the SS is used the zero-lag correlations and autocorrelations show an immediate decay for all the roughness lengths considered, which in this case is due to the fact that the turbulence is instantly destroyed by the solver in a short distance due to the incorrect correlations of the synthetic velocity field.
Finally, the power spectral density (PSD) of the streamwise velocity S uu (f) at three different locations in the successor domain along the streamwise direction (x/H ¼ 0, x/H ¼ 2 and x/H ¼ 4, at z/H ¼ 0.23) is reported for the three methods and the three cases (Fig. 16). It is clear that for the PM the spectrum is not subjected to significant variations along the streamwise direction, which is a further indication of homogeneity of the fluctuations in x-velocity since the integral of the PSD along the frequency range represents the variance of the monitored variable. The evolution of the spectra is different for the synthetic methods. Indeed, the need for an adaptation length is visible for the VM and the SS, whose PSDs at x/H ¼ 2 and x/H ¼ 4 differ significantly from those at the inlet. It is also evident that the energy content of the spectrum for the SS downstream of the inlet is one order of magnitude lower than for the other two cases. This is also reflected in the lower amplitude of fluctuations in the velocity time-traces. As an example, the streamwise velocity time histories for a point at x/H ¼ 2, obtained with the three inflow methods considered, are shown in Fig. 17. Given the lower amplitude of the power spectral density S uu when the SS is used, the variance of the signal is consequently lower than the other two cases. As a matter of fact, the PSDs related to the PM and the VM have comparable amplitude. As a result, the variance of the velocity signal when the PM and the VM are employed is higher than for the SS. The typical -5/3 slope hypothesized by Kolmogorov in the inertial subrange is observed up to f 40 Hz for the PM. For the synthetic methods, in the same frequency range the slope is steeper than -5/3, in line with what has been observed by Tabor and Baba-Ahmadi (2010). For higher values of the frequency a faster decay of the spectrum is present. This is also reported in previous LES studies ( Thomas and Williams, 1999;Xie et al., 2004), and is hypothesized by these authors to be due to the dissipation of the discretization schemes and the finite length of the top-hat filter, respectively.

Influence of aerodynamic roughness length
For the cases under examination it is found that the value of the aerodynamic roughness length does not play an important role for the inflow methods considered, since the homogeneity of the physical quantitites is almost unaffected by this parameter. Only a modest impact is reported for the VM in terms of homogeneity of the turbulence kinetic energy, for which the deviation from the target profile increases from 18.5% (case 1) to 22.3% (case 2) up to 23.5% (case 3) at x/H ¼ 4. When the cross-sectional plane at x/H ¼ 0 is considered the same trend is observed, i.e. the deviation from target profiles of turbulence kinetic energy increases with increasing roughness length (from 25% to 36.4%). This is not observed when the plane at x/H ¼ 2 is considered, where an increase in the mean absolute deviation of turbulence kinetic energy from target profiles occurs when roughness is increased between case 1 and case 2 (21.6-27.2%), but not between case 2 and case 3 (26.7%). It can be concluded that the aerodynamic roughness length, which in the reported simulations is implicitly taken into account by means of a  surface shear-stress boundary condition, has only a weak influence on the performance of the three inflow methods considered.

Discussion and conclusions
Large eddy simulations of the neutral atmospheric boundary layer are performed at wind tunnel scale using three inflow generation methods: a precursor method and two basic synthetic methods, i.e. the Vortex Method and the Spectral Synthesizer. The performance of each method is assessed for three aerodynamic roughness lengths, corresponding at full scale to rural, suburban and urban terrain. The main results are summarized below: The analysis of the flow field in the precursor domain shows that in the central part of the domain, the most relevant for CFD simulations, an acceptable level of spanwise homogeneity is present. As a matter of fact, the deviations of the mean streamwise velocity from the spanwise-averaged values are limited in that region, whereas larger mean absolute deviations are found at the sides of the domain, with a maximum for case 3 (10.8% at y/H ¼ 0.75). In agreement with findings from previous studies, the PM is the most suitable method to preserve the homogeneity of the flow in the streamwise direction, with a discrepancy at the outlet section with respect to target profiles lower than 1% for mean velocity, and 2.8%, 3.1% and 3.5% for the turbulence kinetic energy for cases 1, 2 and 3 respectively. The zero-lag correlations and the autocorrelations show a realistic decay over a relatively short distance. Among the synthetic methods, the VM outperforms the SS. The turbulence kinetic energy decay from the inlet to the outlet is 18.5%, 22.3% and 23.5% for the former and 79.0%, 79.1% and 79.6% for the latter for cases 1, 2 and 3 respectively. It is worthwhile to note that since the application of the VM does not necessarily require a precursor simulation, it is more flexible and computationally less expensive than the PM, though less accurate. Eventually, the VM may represent a good trade-off between accuracy and computational requirements, provided that: i) a sufficiently long development section in the computational domain is present along which fluctuations can evolve towards the prescribed profiles; ii) deviations of mean velocity and turbulence kinetic energy in the lower part of the domain are considered to be acceptable for the case under examination. Among the three methods tested, the implemented version of the SS in the solver does not allow for an accurate reproduction of the turbulence fluctuations in an empty domain compared to the PM and the VM. The application of modified versions of the SS (e.g. Aboshosha et al., 2015a;Castro et al., 2011;Huang et al., 2010) that were developed for improved performance (i.e. higher quality of the synthesized turbulence, which satisfies prescribed statistics and spectra) might therefore also improve the streamwise homogeneity of the vertical profiles of mean velocity and turbulence kinetic energy. Nevertheless, testing improved inflow methods is not the purpose of this research, since this has already been done by other authors (see Section 2.3). Roughness plays a minor role in the performance of the inflow methods. In particular, differences in the decay of turbulence kinetic energy between different aerodynamic roughness lengths are very small when the PM and SS are applied. The impact on the performance of the VM is slightly higher than for the other two methods, with larger deviations from target profiles being observed with increasing roughness. As a consequence, for the cases under examination, two important aspects of the LES of the atmospheric boundary layer, i.e. the inflow turbulence generation method and the value of aerodynamic roughness length, are found to be only weakly coupled and can essentially be dealt with separately.
In conclusion it is also important to highlight some limitations of this study, in particular: Only one value for the damping factor of the wall shear stress boundary condition is considered. Nevertheless, optimizing this coefficient is not the objective of this paper. For the same reason only one subgrid-scale model is employed. Further simulations including different values of the coefficient β and different subgrid-scale models could clarify the influence of these two aspects on the performance of the inflow methods. The effect of the streamwise dimension of the successor domain on the performance of the three inflow methods is not assessed. The role of this parameter may be understood from the zero-lag correlations and the evolution of turbulence parameters in the successor domain. For the PM and the SS it plays a minor role, since in both cases (for different reasons) the autocorrelations drop to zero before the outlet section and, in addition, in both cases the discrepancy with target profiles of velocity and turbulence kinetic energy slightly increases at the outlet section compared to the middle plane (x/H ¼ 2; see Tables 4 and 5). In contrast, for the VM: i) a longer distance than for the PM is required to damp the correlation coefficients; ii) the discrepancy from target profiles decreases along the streamwise direction. Therefore, it is reasonable to assume that increasing the longitudinal dimension of the domain would lead to lower deviations of the velocity fluctuations from target profiles when the VM is employed. Nevertheless, in order to correctly assess the effects of this parameter, further numerical tests, in which different lengths of the domain are considered, need to be performed.