{'fraction_non_alphanumeric': 0.026949918799104595, 'fraction_numerical': 0.011455910108414168, 'mean_word_length': 4.734709287691921, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Wind turbine wakes in forest and neutral plane wall boundary layer large-eddy simulations

Wind turbine wake flow characteristics are studied in a strongly sheared and turbulent forest boundary layer and a neutral plane wall boundary layer flow. The reference simulations without wind turbine yield similar results as earlier large-eddy simulations by Shaw and Schumann (1992) and Porte-Agel et al. (2000). To use the fields from the homogeneous turbulent boundary layers on the fly as inflow fields for the wind turbine wake simulations, a new and efficient methodology was developed for the multiscale geophysical flow solver EULAG. With this method fully developed turbulent flow fields can be achieved upstream of the wind turbine which are independent of the wake flow. The large-eddy simulations reproduce known boundary-layer statistics as mean wind profile, momentum flux profile, and eddy dissipation rate of the plane wall and the forest boundary layer. The wake velocity deficit is more asymmetric above the forest and recovers faster downstream compared to the velocity deficit in the plane wall boundary layer. This is due to the inflection point in the mean streamwise velocity profile with corresponding turbulent coherent structures of high turbulence intensity in the strong shear flow above the forest.


Introduction
Efficient atmospheric boundary layer simulations in the presence of the forest are crucial for understanding the wind turbine wake flow above canopies. Realistically turbulent boundary layers with wind turbine wakes are crucial to understand for a number of reasons: as a basis for sound propagation studies, calculations of loads on turbine blades, and for the for wind turbine site assessment. Since the pioneering large-eddy simulation (LES) of Shaw and Schumann [1], the forest stands have been treated as a porous body of horizontally uniform (leaf) area density a(z) with constant drag coefficient. This approach is established to simulate forest canopy turbulence [2] above horizontally homogenous forested terrain, forest edge flows [3], as well as stratified canopy flows [4] and is sometimes called field-scale approach. Finer scale and field campaigns consider the heterogeneity of the canopy at the plant-scale [5] as done before for urban flow [6]. We introduce a newly developed method for wind turbine wake simulations above forests with EULAG (Eulerian/semi-Lagrangian geophysical fluid solver) [7].
A wide variety of turbulent atmospheric boundary layers exists and each boundary layer exhibits a certain type of turbulence structure [8]. For our purpose it is crucial that the LES reproduces these structures as well as the turbulence statistics of the respective boundarylayer flow. The domain averaged wind profile U (z) in the neutral boundary layer follows the logarithmic law. In the forest boundary layer a characteristic mixing layer flow evolves with an inflection point in the mean velocity profile at the canopy top [2]. This inflection point initiates an instability where organized flow structures in the shear flow evolve [9] and transient to a turbulent state. Horizontally averaged momentum flux, velocity variance, skewness and the ratio of sweeps and ejections from the turbulence statistics have been evaluated in large-eddy simulations at the field-scale [1] and in plant-scale simulations [5]. These physical quantities compare favorably between large-eddy simulations and measurements [4]. In field measurements of forest canopy turbulence, the streamwise u-and vertical w-velocity fluctuations follow a characteristic k −5/3 Kolmogorov scaling above the trees [2]. The k −1 outer-flow scaling and a Kolmogorow scaling for larger wave numbers was found experimentally [10] in laboratory measurements as well as in the large-eddy simulations [11] of neutral plane wall boundary layer flows. The velocity spectra in the neutral plane wall boundary layer are known to exhibit a self-similar inverse scaling with height z [10] as the turbulent structures develop from the ground up to larger heights. Streaks of relatively low streamwise momentum are found on the ground [8]. Hairpin vortices form around these low momentum regions [10]. Some of these vortices grow in size due to vortex stretching, tilting and twisting [9]. At larger heights, a cascade of turbulent structures over a larger spectrum of wave numbers is found.
Replicating the mean fields and the turbulence structure at the inlet of a large-eddy simulation independent of the outlet of the simulations is a requirement, e.g., for simulating microscale phenomena as internal boundary layers [12], forest edge flows [3], or for wind-turbine wake flows [13]. One way is to prescribe stochastic fluctuations on the mean inflow velocity profile at the inlet in large-eddy simulations [14]. Such stochastic fluctuations do not necessarily exhibit the characteristics of fully developed turbulence, e.g., when the fluctuations cannot reproduce the shape, size and orientation of eddies present in the atmospheric boundary layer [12]. Therefore, Mayor et al. [12] and Porté-Agel et al. [15] rescale the simulated flow structure from the end of the domain and introduce this flow in the inflow with the so called fringing technique [16]. Besides the non-isotropy and inhomogeneity of the wake turbulence that is then introduced at the inflow, the rescaling also yields to an imbalance of the flow quantities [16]. This requires a transition region [12] up to ten boundary layer thicknesses [16] until the flow reaches a quasi stable equilibrium turbulent state. To circumvent this, Churchfield et al. [13] use the temperature and velocity flow fields from the inlet of a so called precursor simulation that is stored at every time step of the numerical realization. Recently, Stevens et al. [17] applied a methodology to simulate the flow of the precursor simulation with spectral discretization at the same time parallely in another computational domain by transferring the required flow fields at each time step to the inflow of the wind turbine wake flow simulation by processor communication with the message passing interface (MPI). In the inflow, the flow fields need to be adjusted to the fields from the precursor simulation.
We introduce and discuss a similar and efficient computational method for finite-difference discretization with EULAG. The method avoids additional MPI communication or disk input/output (I/O) operations. We use open boundaries and no flow adjustment or fringing is required to provide the inlet of the simulation with a fully turbulent flow which characterizes the respective boundary layer regime.
We address the following questions: • Is the boundary layer turbulence simulated correctly upstream of the wind turbine?
• How does the structure of a wind turbine wake differ in both boundary layers?
• How do turbulence intensity and eddy dissipation rate evolve in the wakes?
After explaining the methodology of the two hydrodynamic solvers in EULAG in section 2, we present an intercomparison of the simulated wind turbine wakes and the boundary layers in section 3 and conclude in section 4.

Method
For the purpose of this experiment, EULAG is extended to accommodate two independent hydrodynamic solvers: a and b with two computational domains A and B (Figure 1). Each solver carries its own set of dependent variables. The two solvers are integrated simultaneously, can run independently or can be coupled as done here. Solver a simulates a fully developed homogeneous turbulent boundary layer, while solver b simulates a heterogeneous boundary layer with a wind turbine wake in the domain. In a homogeneous boundary layer, as a plane layer or a forest layer, solver a computes fully developed continuously turbulent fields of u, v, w, and e with correct turbulence statistics in a temporally evolving LES with cyclic boundary conditions. To simulate a heterogeneous boundary layer (with wind turbine wake) requires independent inflow and outflow boundary conditions in the spatial LES of solver b. All developed fields and moments from the simulation of solver a are copied to solver b as initial conditions when solver b starts. During the run all flow fields from A are exactly reproduced in B at the inlet including, e.g., turbulent kinetic energy e. The exact replication of the fields at the inlet from solver a leads to the correct turbulence statistics of the respective boundary layer in the inflow of the wind turbine wake flow in solver b.

Governing equations
In the governing equations of the large-eddy simulations with EULAG, the Boussinesq approximation (1) is adopted for the neutral boundary layer flow. To achieve optimal computational efficiency, we use the three-dimensional domain decomposition with the MPI [19]. Each solver has full access to the local flow fields ψ(x, y, z, t) of both solvers at each time step t without MPI communication, as the two equation systems a and b are solved locally on each processor. The symbol ψ represents velocity components u, v, w, or turbulent kinetic energy e. Within this experiment, the flow fields at the inlet boundary of solver b, i.e., ∂B| x=0 reproduce the flow fields evaluated within solver a on the leftedge at time t, i.e., The upper indices a,b refer to solver a and b, respectively. Solver b begins to integrate the governing equations when the flow in solver a is fully turbulent. At this state, the initial conditions for the governing equations in solver b are given by the flow fields calculated in solver a.
Here, the equations for the field-scale simulations are presented, where v a,b denotes the velocity vector and d/dt = ∂/∂t + v a,b · ∇ is the total derivative. The density-normalized pressure is π a,b = p a,b /ρ 0 with Boussinesq reference density ρ 0 . Dissipation in the momentum equation (2) is represented with D(v a,b ). The wind turbine is simulated with the actuator disk approach incorporating the drag force of a wind turbine with axial induction factor 0.25 as in previous simulations [18,20]. The force of the wind turbine acting on the atmosphere is the time-filtered flow velocity at the disk and the drag coefficient c T = 4/3 is calculated at the disk accordingly to momentum theory [18,21,22]. The sub-grid scale model for turbulent kinetic energy [23] is represented as S(e a,b ) in equation (3) and applied for a wide range of scales and scenarios. The forest model by Shaw and Schumann for neutral stratification parameterizes the forest as a porous medium with drag forest drag coefficient c f or = 0.15 and leaf area density a(z) in both solvers. The equation for turbulent kinetic energy incorporates a sink −2 e a,b τ −1 over an inverse time scale τ −1 = c f or a(z)V a,b . This sink represents the dissipation of the turbulent kinetic energy by, e.g., smaller scale wake vortices of tree branches.

Flow configuration
The physical processes are modeled as follows: momentum flux on the ground is parameterized with a surface drag coefficient c surf = 0.01 and a resulting surface stress proportional to c surf the drag coefficient and the magnitude of the horizontal velocity (u a,b ) 2 + (v a,b ) 2 at the lowest mesh point [24]. A zero flux condition is applied at the top lid. The domain averaged streamwise velocity is 10 m s −1 . The simulations are carried out for more than 30 through-flows, until the turbulent flow approaches a quasi stable equilibrium. In the simulations with a forest the canopy has a height of 20 m and a leaf area index of 2 and the same leaf area density a(z) is used as in the simulations by Shaw and Schumann [1]. The wind turbine is simulated with a hub height of 140 m and a rotor diameter of 100 m in the neutral plane wall and forest boundary layer large-eddy simulations.
The domain size in the numerical simulations is L x =2 km in streamwise direction, L y = 384 m in lateral direction and L z = 300 m in vertical direction with corresponding mesh points (n, m, l) = (1024, 196, 150). Physically, the domain top corresponds to an inversion in the atmosphere [1]. Additional simulations were carried out by doubling the domain height. All simulations are resolved with an isotropic Cartesian mesh of 2 m mesh width and a temporal resolution of 1/16 s. The domain boundaries are cyclic in the horizontal in solver a. In solver b the boundaries are cyclic in y-direction and reproduce the flow fields from solver a at the inflow boundary with a Dirichlet condition for all fields except for the pressure that is solved with the Neumann boundary condition. The outflow boundary in solver b adjusts to a gradient free flow with respect to the condition ψ b | N +1 i=n = ψ b | N i=n−1 in discretized form, where N + 1 is the next, N the current time step and i is the running index in x-direction.

Results
At the field scale, our new methodology with the two hydrodynamic solvers is applied to a forest boundary layer and to a neutral plane wall boundary layer flow, both without and with a wind turbine wake.  Fully developed turbulence is simulated to a quasi stable equilibrium in solver a. The flow structure becomes evident in instantaneous velocity fields as, e.g., in the streamwise velocity fields u presented for both boundary layers ( Figure 2). The velocity fields u(x, y, z 0 ) one mesh point above the ground at z 0 = 2 m in the plane wall boundary layer and one mesh point above the canopy top at z 0 = 22 m for the forest canopy flow is shown for the last simulated time. Different size, shape, and orientation of the simulated turbulent structures in these boundary layers is evident.
In the neutral plane wall boundary layer, regions of lower momentum are elongated in streamwise direction [8,10]. These structures are the so called streaks [8] and are evident at several locations with a characteristic width of about 50 m and a length of approximately 500 m in the simulated plane wall boundary layer near the ground (Figure 2, top). Directly aloft the forest canopy top, the velocity field u exhibits a different structure. The simulated streaks of relatively low momentum are wider in lateral direction and occur more frequently in the forest boundary layer. The amplitudes of the fluctuations of the streamwise velocity are more pronounced in the turbulence evolving above the forest (Figure 2, bottom). The fluctuations appear over a wider range of scales in this forest canopy flow than in the neutral plane wall boundary layer. It is evident that the structure of the simulated streamwise velocity is multiscale. This multiscale nature of the flow fields above forest canopies was pointed out as a characteristic property in forest canopy large-eddy simulations by Shaw and Schumann [1].

Spectra of the simulated streamwise u-velocity fields in both homogeneous boundary layers
To characterize the scaling properties of the turbulent structures at the last simulated time of both boundary layers, spectra of the streamwise velocity are explicitly calculated and compared with previous results from literature [10,11]. Laboratory experiments and previous large-eddy simulations indicate a scaling of −1 for the largest scales and smallest wave numbers in the neutral plane wall boundary layer [10]. This scaling is simulated for both, the neutral plane wall and forest boundary layer for wave numbers k 1 ≤ 0.3 z (Figure 3). Within the inertial subrange, the velocity variance φ 11 scales with the classical −5/3-slope. This scaling can be  Figure 3. Spectra of the streamwise velocity variance φ 11 as explained in the text below. The spectra are calculated for the plane wall (top) and forest canopy boundary layer (bottom). The color scale ranges from green for a height z/H = 0.05 over yellow to red with a maximum height of z/H = 0.95, where H is the boundary layer height. explained with dimensional analysis and not with the Kolmogorov theory as the flow is highly non-isotropic close to the ground [10], where the vertical velocity ceases. Inside the forest in the crown space, where the leaf area density is highest at z/H = 0.05, higher wave numbers and smaller scales carry relatively lower absolute values of variance φ 11 than aloft the forest canopy. As prognosed by Perry et al. [10] and found in previous large-eddy simulations of neutral plane wall boundary layer turbulence [11], all the spectra collapse to one spectrum following the inverse scaling law outside the dissipative range, when scaled properly with the friction velocity u * and height z.
3.3. Structure of the wake flow above the plane wall and forest boundary layer These boundary layers provide the inflow for the wind turbine wake simulations. Compared to the wind turbine wake in the neutral plane wall boundary layer (Figure 4, top), the wake structure as seen, e.g., in the temporal mean streamwise velocity field u (x, y 0 , z), above the forest (Figure 4, bottom) is strongly asymmetric and recovers faster downstream. A reason for earlier wake recovery is the enhanced inflow turbulence intensity evolving due to the inflection point in the mean velocity profile above the forest. The higher turbulence intensity leads to a stronger mixing and relatively earlier wake recovery. The effect of the domain top at 300 m and 600 m was compared. The significantly shorter and vertically more asymmetric wind turbine wake was found in both forest boundary layer simulations.

Properties of the non-equilibrium turbulence in the wind turbine wakes
To compare the recovery of the wake flow between neutral plane wall and forest boundary layer quantitatively, the velocity deficit is defined similarly to previous work in literature [25] as at hub height z hub and lateral position of the hub y hub at the inlet at x = 0, where a temporal mean velocity U hub occurs. The velocity deficit can be estimated with an empirical formula as ∆U (x) = x −0.57 [25]. A number of measurements were compared to this empirical function in a previous paper [25]. The range over all the measurements is indicated in Figure 5. The simulated wind turbine wake velocity deficit values in the neutral plane wall and forest boundary layer occur within this range. The velocity deficit in the neutral plane wall boundary layer reaches the value of 0.1 after 15 rotor diameters. In the forest boundary layer the velocity deficit reaches this value much earlier after ≈ 5 rotor diameters about three times as fast downstream of the wind turbine. Another crucial quantity to compare is the maximum additional turbulence intensity ∆I [20] due to the wind turbine and calculated as at the position of the top shear layer as done by Jimenez [20], where the subscript m denotes the maximum value of the added variance ∆ u 2 and the resulting standard deviation σ u,m is normalized with the local mean inlet velocity U . The resulting values are below the values of an empirical model developed by Frandsen [26] and closer to the values derived from large-eddy simulations of wind turbine wakes in neutral plane wall boundary layer flows by Jimenez [20]. Overall, the values by Jimenez are higher by a maximum of 0.02 as Jimenez simulated the added turbulence intensity of several wind turbines in a column and the simulations here consider one wind turbine wake specifically. The turbulence can be differentiated in different classes. The additional turbulence in the wake evolves locally due to the presence of a shear layer in the wind turbine near-wake. Due to a stronger shear layer and higher amplitude fluctuations in the inflow forest boundary layer turbulence, the additional turbulence intensity in the wake ∆I is  about twice as high in the near wake when immersed in the forest boundary layer, compared to the near wake in the plane wall boundary layer. The corresponding locally higher turbulence intensity dissipates faster in the forest boundary layer, so that the added turbulence intensities agree in both boundary layers within a range of 0.02 after about four rotor diameters. At the height of the shear layer at the wind turbine top, where the added turbulence intensity is maximum, we calculate the eddy dissipation rate from the subgrid scale turbulent kinetic energy e as Scipion [27] does: where is the local mixing length corresponding to one grid increment in our large-eddy simulations. The eddy dissipation rates in the inflows correspond to the simulated dissipation rates in solver a. The turbulence in solver a is in quasi stable equilibrium and thus the simulated dissipation rates are constant in downstream direction. is one order of magnitude higher in the simulated forest boundary layer compared to the plane wall boundary layer. The values of eddy dissipation rate in the wake and inflow agree with measurements from Lundquist and Bariteau [28] of between 10 −4 and 10 −2 m 2 s −3 in the plane wall and forest boundary layer turbulence, respectively. As in the measurements, the dissipation rate is strongly enhanced in the wake by one to two orders of magnitude to values above 10 −2 m 2 s −3 . The maximum values of agree in both wind turbine wake boundary layers. The dissipation rate reaches an equilibrium state in the simulated forest boundary layer after 10 rotor diameters downstream of the wind turbine. There, the is 0.1 m 2 s −3 higher than in the inflow forest boundary layer. This equilibrium value of is not reached in the plane wall boundary layer, where the dissipation rate in the wake is still one order of magnitude higher than in the inflow. This result emphasizes the importance of modeling the wind turbine wake in a separate numerical domain with non-cyclic boundary conditions.

Conclusions
A novel computational approach is introduced and successfully applied to simulate homogeneous and heterogeneous neutral plane wall and forest boundary layer flows with EULAG. The two hydrodynamic solvers in EULAG allow to exactly reproduce the turbulence structure and statistics from the fully developed horizontally homogenous boundary layer at the inlet of the heterogeneous wind turbine simulation with open boundary conditions. The simulations reproduce the characteristic spectra for the variance of the streamwise velocity in the plane wall. The spectra in the forest boundary layer compare well. We find that the wind turbine wake recovers significantly earlier and exhibits a stronger vertical asymmetry in the forest boundary layer. These conclusions are found to be independent of domain height. The earlier recovery is due to stronger turbulent mixing above the forest and suggests a possible nearer placement of downstream wind turbines. Hence, relative to U hub , more wind energy per unit area can be harvested, however, at the cost of stronger turbulent fluctuations acting on the wind turbine blades above the forest.