Wake properties and power output of very large wind farms for different meteorological conditions and turbine spacings: A large-eddy simulation case study for the German Bight

Germany’s expansion target for offshore wind power capacity of 40 GW by the year 2040 can only be reached if large portions of the Exclusive Economic Zone in the German Bight are equipped with wind farms. Because these wind farm clusters will be much larger than existing wind farms, it is unknown how they affect the boundary layer flow and how much power they will produce. The objective of this large-eddy-simulation study is to investigate the wake properties and the power output of very large potential wind farms in the German Bight for different turbine spacings, stabilities and boundary layer 5 heights. The results show that very large wind farms cause flow effects that small wind farms do not. These effects include, but are not limited to, inversion layer displacement, counterclockwise flow deflection inside the boundary layer and clockwise flow deflection above the boundary layer. Wakes of very large wind farms are longer for shallower boundary layers and smaller turbine spacings, reaching values of more than 100 km. The wake in terms of turbulence intensity is approximately 20 km long, where longer wakes occur for convective boundary layers and shorter wakes for stable boundary layers. Very large wind 10 farms in a shallow, stable boundary layer can excite gravity waves in the overlying free atmosphere, resulting in significant flow blockage. The power output of very large wind farms is higher for thicker boundary layers, because thick boundary layers contain more kinetic energy than thin boundary layers. The power density of the energy input by the geostrophic pressure gradient limits the power output of very large wind farms. Because this power density is very low (approximately 2 W m−2), the installed power density of very large wind farms should be small to achieve a good wind farm efficiency. 15


Introduction
At present, the global installed wind power capacity from offshore wind farms is increasing rapidly. According to the expansion targets of the current leading offshore wind markets (the United Kingdom, Germany and China), the offshore wind power capacity will be subject to significant growth over the next decades. The German expansion target for offshore wind power capacity is 40 GW by the year 2040, which is more than the global installed offshore wind power capacity of 32.5 GW in the 20 year 2020 (WindSeeG, 2020;Herzig, 2020). The otherwise undisturbed flow at offshore sites will be increasingly modified by wind farms, affecting the wind farm power output but also the meteorological conditions in the wake. For wind farms with a state-of-the-art size of approximately 100 turbines and a length of approximately 5 km, these effects have been extensively investigated experimentally and numerically and are generally well understood. Instead of using an idealized wind farm shape, we investigate a potential future wind farm scenario in the German Bight, which is shown in Fig. 1. The scenario assumes that all priority areas for future wind farms are equipped with 15 MW wind turbines. This results in a total number of up to 2088 wind turbines with a total wind farm capacity of up to 31 GW. More than 7 billion grid points are required to fill the large domain with a turbine wake resolving grid. To our knowledge, this large-eddy 60 simulation case study exceeds other studies in terms of wind farm area and total wind turbine number by at least one order of magnitude.
The numerical model, setup and boundary conditions are described in section 2. The simulation results regarding the wake properties and the power output are shown and discussed in section 3. Section 4 concludes and discusses the results of the study.

Numerical model
The simulations were performed with the PArallelized Large-eddy simulation Model PALM (Maronga et al., 2020), which is developed at the Institute of Meteorology and Climatology of Leibniz Universität Hannover, Germany. Several wind farm flow investigations have been successfully conducted with this code in the past (e.g., Witha et al., 2014;Dörenkämper et al., 70 2015). PALM solves the non-hydrostatic, incompressible Navier-Stokes equations in Boussinesq-approximated form, spatially filtered over a grid volume. The equations for the conservation of mass, momentum and internal energy then read as: where an overbar indicates filtered quantities and a double prime indicates subgrid-scale (SGS) quantities, i, j, k ∈ {1, 2, 3}, u i , u j , u k are the velocity components in the respective directions (x i , x j , x k ), t is time, f i = (0, 2Ω cos(φ), 2Ω sin(φ)) is the Coriolis parameter with the Earth's angular velocity Ω = 0.729×10 −4 rads −1 and the geographical latitude φ. The geostrophic 80 wind speed components are u g,j and the basic state density of dry air is ρ 0 . The modified perturbation pressure is π * = p * + 2 3 ρ 0 e, where p * is the perturbation pressure and e = 1 2 u i u i is the SGS turbulence kinetic energy. The gravitational acceleration is g = 9.81 ms −2 and δ is the Kronecker delta.
The SGS model uses a 1.5-order closure according to Deardorff (1980), modified by Moeng and Wyngaard (1988) and Saiki et al. (2000). Recently, the modified version of Dai et al. (2021) has been implemented in PALM, which allows for coarser grid 85 spacings in stable boundary layers due to reduced grid spacing sensitivity. This modified version is used for the simulation of wind farms in a stable boundary layer.
The following features of PALM are relevant for the performed simulations. It is possible to prescribe a surface heating or cooling rate, instead of prescribing a surface heat flux. Stable boundary layers can also be generated by imitating warm air advection by using a large-scale forcing. Convective boundary layer growth can be compensated by applying a large-scale 90 subsidence to the potential temperature field. A Rayleigh damping layer can be used in order to avoid gravity wave reflection at the top of the domain.
The wind turbines are represented by an advanced actuator disc model (ADM) that acts as an axial momentum sink and an angular momentum source (inducing wake rotation). The ADM is described in detail by Steinfeld et al. (2016). The disc element forces are distributed to the neighbouring grid points by a smearing kernel, which causes a power overestimation of 95 12.5 % for a smearing kernel radius of 1∆x. The wind turbine power output is corrected for the power overestimation before entering the wind farm power output analysis 4 https://doi.org/10.5194/wes-2021-83 Preprint. Discussion started: 23 August 2021 c Author(s) 2021. CC BY 4.0 License.

Case selection
To produce meaningful and relevant results, the simulations should represent the most common meteorological conditions in the German Bight. A climatology with frequency distributions of wind speed, wind direction, boundary layer (BL) height and 100 stability information extracted from the COSMO-REA6 reanalysis dataset can be found in Appendix A. The analysis was provided by Thomas Spangehl (German Weather Service) and it is based on hourly data of a 24-year period  at 54 • 30 N, 6 • 00 E, which is located inside Zone 3 (cf. Fig. 1). Wind speed and direction are evaluated at 178 m height, which is the closest COSMO model level to the hub height of 150 m of the wind turbine used in the simulations.
Due to the high computational cost per simulation, only a limited number of simulations were carried out. This study consists 105 of five simulations with varying stability, turbine spacing and BL height. An overview is given in Table 1. Two cases with a neutral boundary layer (NBL), two cases with a convective boundary layer (CBL) and one case with a stable boundary layer (SBL) are simulated. In the two NBL-cases NBL-700-5D and NBL-700-7D the turbine spacing is set to s = 5D and s = 7D, where D is the rotor diameter of the turbine. The turbine spacing for all other cases is s = 7D. The NBL is capped by an inversion layer with a lapse 110 rate of Γ = +1 K km −1 to achieve a BL height of approximately 700 m, which is a very common BL height in the German Bight, according to the COSMO-REA6 climatology (cf. Fig. A3 and A4). The correct term for such an inversion-capped NBL is conventionally neutral boundary layer (CNBL). However, the cases are named NBL-700-7D and NBL-700-5D to avoid confusion with the CBL-cases.
Because CBLs are more frequent and are generally thicker than SBLs in the German Bight, two CBL-cases CBL-700-7D 115 and CBL-1400-7D with a BL height of h ≈ 700 m and h ≈ 1400 m, respectively, are simulated. This represents the spread of CBL heights in the German Bight (cf. Fig. A3). Note that a CBL is the only BL type for which the BL height can be controlled freely by the initial temperature profile without the need to change other parameters. The (steady-state) BL height of CNBLs and SBLs can not be controlled directly but is rather a function of friction velocity, Coriolis parameter, free atmosphere (FA) stratification and surface buoyancy flux (Zilitinkevich et al., 2007).

120
The BL height of the SBL-case SBL-300-7D is h ≈ 300 m, so that the wind turbines with a rotor top height of 270 m are still within the BL and do not penetrate into the FA. 300 m is a small but still typical value for an SBL in the German Bight (cf. Fig. A4).
The wind speed at hub height is set to 10 ms −1 for all cases. This wind speed is less than the mean wind speed in the German Bight (10.8 ms −1 , cf. Fig. A1) to stay below the rated wind speed of v rated = 10.59 ms −1 of the IEA 15 MW reference wind 125 turbine (Gaertner et al., 2020). Thus the turbine operates at a high thrust coefficient and the turbine power is a function of the wind speed. The surface roughness length in all cases is z 0 = 1 mm. The wind direction is set to 225 • , which is one of the most common wind directions in the German Bight. Because the main axis of the wind farm cluster in Zone 3 has a southwest-northeast orientation, strong wake effects can be expected for this wind direction.

130
The domain and wind farm layout are shown in Fig. 2. The domain length and width are L x = 204.8 km and L y = 163.84 km, respectively. These lengths correspond to n x = 10240 and n y = 8192 grid points in x-and y-direction for isotropic grid spacings of ∆x = ∆y = ∆z = 20 m for all cases. These spacings yield a density of 12 grid points per rotor diameter, which is enough to resolve the most relevant eddies inside the wind turbine wakes. Above the BL, where no turbulence must be resolved, the grid is stretched vertically to a maximum of ∆z max = 50 m to save computational cost. The stretch factor is 135 f stretch = ∆z(k + 1)/∆z(k) = 1.08 and the stretching starts at z s (cf. i.e. 204.8 km/10 ms −1 ≈ 5.7 h). The last 4 h are used for the evaluation, e.g. averaging and flux calculations.
At the crosswise lateral boundaries, cyclic boundary conditions are applied and at the outflow plane, radiation boundary conditions are applied. At the inflow plane, steady-state vertical profiles of a precursor simulation are prescribed (details about the precursor simulations are given in the next section). To have a turbulent and stationary inflow from the beginning of the main simulation, the flow field is initialized by the instantaneous flow field of the last time step of the precursor simulation.

145
Because the precursor domain is much smaller than the main domain, the flow field is filled cyclically into the main domain. It is important to note that the width of the main domain is a non-integer multiple of the width of the precursor domain, to trigger the break-up of the unnatural periodicity in y-direction of the flow field, that is introduced by the cyclic fill method.
The turbulent state of the inflow is maintained by a turbulence recycling method, that maps the turbulent fluctuations from the recycling plane at x = x r onto the inflow plane at x = 0 (Lund et al., 1998;Kataoka and Mizuno, 2002). The turbulent 150 fluctuation Ψ (y, z, t) at each time step is defined as the difference between the absolute value Ψ(x r , y, z, t) and the horizontal  line average in y-direction Ψ(x r , z, t) y at that height: where Ψ can be a velocity component, the potential temperature or the subgrid-scale turbulent kinetic energy. The turbulent fluctuation is added to the mean inflow profile Ψ inf low (z) at the inflow plane. Instead of adding it at the same y-location, it 155 can be added at y + y shif t : Ψ(0, y + y shif t , z, t) = Ψ inf low (z) + Ψ (y, z, t) . The application of the y-shift effectively reduces the strength of streamwise elongated streaks in the mean wind speed of NBLs (Munters et al., 2016). 1 The otherwise inhomogeneous inflow with crosswise variations of wind speed of up to 10% would hamper the evaluation of the wind farm power output and wake. A homogeneous inflow wind speed is of utmost importance in 160 wind energy studies, because the wind turbine power is proportional to the third power of the wind speed. The y-shift is chosen in such a way, that the flow is recycled many times before reaching its initial y-position, which is achieved if the least common multiple of the y-shift and the domain width is a large number. The y-shift is also applied to the non-NBL cases, because it reduces crosswise variations of wind speed that are caused by wind farm induced flow blockage. The flow blockage leads to a reduced mean wind speed at some y-locations of the recycling plane, which is "interpreted" as turbulent fluctuation, and thus 165 mapped onto the inflow. Due to the self-reinforcing behaviour of this process, the crosswise variations can build up quickly without y-shift, even if the wind farms have a distance of 15 km to the recycling plane.
The turbulence recycling is limited to a height just above the BL height, so that potential BL growth between inflow and recycling plane will not affect the inflow BL height. The recycling plane is located 10 km downstream of the inflow plane, which gives the turbulent structures enough time to interact and decorrelate before becoming recycled. For the CBL-cases, the 170 absolute value of the potential temperature is recycled instead of its turbulent fluctuation, so that the inflow temperature rises according to the increasing surface temperature. This method is not needed in the SBL-case, because the surface temperature is constant in time (details in the next section).
The priority areas of Fig. 1 are rotated 45 • clockwise, so that the inflow at hub height is parallel to the x-axis for a wind direction of 225 • . The priority areas are filled with a regular array of the IEA 15 MW reference wind turbine, that has a 175 rotor diameter of D = 240 m, a hub height of z hub = 150 m and a rated power of P rated = 15 MW (Gaertner et al., 2020).
The wind turbines are staggered, i.e. every second column is shifted by half a turbine spacing in y-direction (cf. Fig. 2). The staggered configuration represents the real-world variation of wind directions better than the very special case of an aligned configuration. Additionally, power output and wake strength are less sensitive to potential wind direction changes (that might occur further downstream inside the wind farm) for the staggered configuration, as revealed by own test simulations with 180 smaller wind farms. The turbine spacing in the x-and y-direction is the same (s x = s y = s). The total number of wind turbines is n wt = 1063 for s = 7D and n wt = 2088 for s = 5D, resulting in a total installed wind farm capacity of 15.9 GW and 31.3 GW, respectively. With a total wind farm area of 3000 km 2 , the resulting installed power density is P 7D = 5.3 MW km −2 and P 5D = 10.4 MW km −2 . Note that s = 7D and P = 5 MW km −2 are typical values for currently existing wind farms in the German Bight but that even with s = 5D the total installed wind farm capacity stays below the 2040-expansion target 185 of 40 GW. Note also that, for the sake of simplicity, all existing wind turbines in the priority areas are replaced by the much larger 15 MW wind turbine.
1 Elongated, streak-like structures in the instantaneous streamwise wind speed (also called superstructures or Very-Large-Scale Motions) are a natural phenomenon of NBLs. However, these structures can be as large as 20 times the BL height (Fang and Porté-Agel, 2015), so that they can not be captured between inflow and recycling plane. Thus, the same structure is recycled repeatedly without breaking up or moving in y-direction. As a result, streaks of high and low wind speed appear in the averaged velocity field even for very long averaging times. A y-shift does not avoid the appearance of streaks in the instantaneous velocity field but due to the changing y-location of the streaks the strength of the streaks in the mean velocity field is reduced effectively (Munters et al., 2016).

Precursor simulations
Steady-state inflow profiles and a turbulent flow field for each main simulation are obtained by a precursor simulation with cyclic boundary conditions in both lateral directions. In order to save computational time, the precursor domains are much 190 smaller than the main domain (cf. Table 1). The domain sizes are different for the different cases in order to ensure that the largest structures of each BL type are covered several times. The number of vertical grid points, the stretching and Rayleigh damping levels are the same as in the corresponding main simulation. It is important that the turbulence and the mean flow are stationary at the end of the precursor simulation. If the mean flow that is prescribed at the inflow plane is not in steady state it will try to reach it during its passage through the main domain, causing streamwise changes in mean quantities such 195 as wind speed and direction. While steady-state turbulence is reached after only a few hours, achieving a steady-state mean flow can take several days due to the slow decay of the inertial oscillation, which has a period of 14.6 h at a latitude of 55 • N.
The physical simulation times of the precursor simulations are 96 h for the cases NBL-700-7D, NBL-700-5D and CBL-1400-7D, 48 h for the case CBL-700-7D and 24 h for the case SBL-300-7D. The initial velocity and potential temperature field is horizontally homogeneous. Horizontal velocity components u and v are set to the geostrophic wind components u g and v g at is used. Further setup details vary significantly between the different cases and hence are described separately in the following sections.

NBL
The initial potential temperature profile of the NBL-cases is linear and has a vertical temperature gradient (lapse rate) of Γ = +1 K km −1 from the surface to the domain top (cf. Fig. 3). At the surface, a Neumann-condition for the potential temperature 210 is applied and the surface heat flux is set to zero. Shear-driven turbulence production leads to the formation of a neutrally stratified BL that grows until it reaches a steady BL height of 780 m. The BL height is defined as the height at which the shear stress reaches 5 % of its surface value. The conventionally neutral boundary layer is separated from the FA by a capping inversion that has a stronger stratification than the FA.

215
The initial temperature profile of the CBL-cases consists of a constant potential temperature between the surface and the desired BL height h = 700 m or h = 1400 m for the cases CBL-700-7D and CBL-1400-7D, respectively. Above that height the potential temperature has a constant lapse rate of Γ = +3.5 K km −1 , which corresponds to the International Standard Atmosphere. A Dirichlet-condition is applied for the surface temperature and a constant surface heating rate ofθ 0 = +0.050 K h −1 andθ 0 = +0.025 K h −1 is used to drive the CBL of the cases CBL-700-7D and CBL-1400-7D, respectively. The heating rates 220 differ by a factor of 2 to achieve approximately the same surface heat flux Q 0 and Monin-Obukhov length L (cf. Table 1), so that only the effect of a changing BL height is seen in the results.
Boundary layer growth is avoided by applying a large-scale subsidence that acts only on the potential temperature field. The subsidence velocity is zero at the surface and increases linearly to its maximum value w sub at the height h and is constant above. The subsidence velocity is chosen in such a way that the temperature increase in the FA exactly matches the surface Large-eddy simulations of CBLs are usually driven by a constant heat flux, i.e. a Neumann-condition for the surface temperature. However, we decided to use a Dirichlet-condition, because of two reasons: -It allows for spatial variations in the surface heat flux which may be caused by enhanced mixing inside the wind farms.

230
In reality, the resulting change in sea surface temperature (on the scale of hours) would be very small due to the good turbulent mixing inside the ocean mixed layer during strong winds and due to the high heat capacity of water in contrast to that of air. Thus, it is more realistic to prescribe a horizontal homogeneous surface temperature than a horizontal homogeneous heat flux.
known in advance and thus the subsidence velocity required for obtaining a constant BL height is also known in advance and does not have to be found iteratively.

SBL
The initial potential temperature profile of the SBL-case is linear and has a vertical temperature gradient (lapse rate) of Γ = +3.5 K km −1 from the surface up to the domain top. A Dirichlet-condition is applied for the surface temperature, as this is 240 the correct surface forcing method for SBLs (Basu et al., 2008). Generating a steady-state SBL is not as simple as it is for the CBL. A straight-forward method would be to use a surface cooling rate. However, due to the long simulation time, required for the decay of the inertial oscillation, the elevated inversion at the top of the SBL would become unrealistically strong (Kosović and Curry, 2000). We developed a method to generate a steady-state SBL in which the potential temperature profile is constant in time and the strength of the elevated inversion can be freely adjusted.

Data analysis
Statistical data that are presented in the results section are obtained in the last 4 h of the 10 h main simulations. Temporal 260 averages are denoted by an overbar (e.g. v h ) and horizontal averages by angled brackets (e.g. θ ). The temporal averaged horizontal wind speed v h is calculated as the average of the absolute values of the wind vector: Resolved turbulent fluxes of momentum are calculated with the temporal eddy-correlation (EC) method. The correlation of two turbulent quantities (e.g. u = u−u and w = w−w) can not be calculated directly during the simulation, because the respective mean quantities are not known in advance. However, the resolved turbulent flux can be calculated after the simulation if the correlation of the absolute quantities is calculated during the simulation: The presentation and discussion of the results is divided into the two sections wake properties and power output. In the first 270 section, the wake properties of very large wind farms and their effect on the BL flow is discussed. In the second section it is discussed how the power output of very large and small wind farms is affected by the variation of the turbine spacing and the meteorological conditions. which is 1 % of the inflow wind speed.
For the CBL-cases, the wind speed is reduced to 6.5 m s −1 for h = 700 m and 8 m s −1 for h = 1400 m inside the large wind farms in Zone 3. Also, the wake length of the large wind farms is much longer for the shallow BL than for the thick BL. This BL height dependency occurs because a thicker BL contains more kinetic energy that can be transported down to the wind 295 turbine level by turbulent vertical mixing than a shallow BL. The wake length and speed deficit of small wind farms (e.g. N-1 and N-2) is relatively unaffected by the BL height because the wind farm induced internal boundary layer does not reach the top of the BL.
The wake is deflected counterclockwise in the CBL-cases, as well. The deflection angle is approximately 5 • for the case with the shallow BL and 1 − 2 • for the case with the thick BL. The higher deflection angle for the case with the shallow BL is 300 caused by a greater speed deficit compared to the case with the thick BL. A comparison between the cases NBL-700-7D and CBL-700-7D shows that the speed deficit inside Zone 3 is greater for the CBL-case. Also, the wakes of the small wind farms are longer for the CBL-case. This is contradictory to the well known fact that wind turbine and wind farm wakes are generally shorter in CBLs than in NBLs and SBLs (Porté-Agel et al., 2020, sec. 2.3 and 3.4.2). To achieve the same hub height wind speed for both cases, the geostrophic wind speed is 6 % greater for 305 the NBL-case (cf. Table 1) than for the case CBL-700-7D. Additionally, the wind speed is supergeostrophic in the upper half of the NBL and thus the mean BL wind speed is approximately 10 % greater in the NBL-case than in the CBL-case. Stability likely has little-to-no affect, because the stratification of the CBL-case is only weakly unstable (L = −420 m, cf. Table 1).
In the stable case SBL-300-7D, the wind speed is reduced to below 7 m s −1 in the first 20 km of the large wind farms in Zone 3. The wind speed deficit is greater and the wake is more than 20 km longer for the small wind farms, compared to the 310 other cases with s = 7 D. The wake of the large wind farms in Zone 3, however, is not longer than in the cases NBL-700-7D and CBL-700-7D. This occurs because the speed recovery in the wake of this large wind farms is not driven by momentum flux divergence (which is stability dependent) but rather by a favorable pressure gradient (details are given in the next section).
The Because the wind farms in this study also have a finite size also in the crosswise direction, it can be seen that the pressure 335 perturbation also significantly affects the wind direction. Due to the streamwise reduction in wind speed, the flow diverges in the crosswise direction inside the wind farms. In the wake, where the flow accelerates, horizontal convergence can be observed.

Reasons for wake deflection and slow speed recovery
What is the reason for the slow speed recovery and the wake deflection inside and behind the large wind farms? In order to answer that question, Fig. 5 shows streamwise (parallel to streamlines) and crosswise (perpendicular to streamlines) components 340 of the pressure gradient force 3 , the Coriolis force F c and the resolved vertical turbulent momentum flux divergence, also called frictional force F f , at z = 150 m and y = 120 km. The pressure gradient force can be divided into the geostrophic pressure gradient force F gp , which is constant and is defined by the geostrophic wind, and the perturbation pressure gradient force F pp , that can vary horizontally due to wind farm induced pressure perturbations. The forces are averaged over 4 turbine spacings along x and 2 turbine spacings along y in order to eliminate peaks in F pp that are caused by single turbines. Thrust forces of 345 the turbines are not included. The analysis is made from a Lagrangian frame of reference, examining the forces on an air parcel during its passage through the wind farms. From an Eulerian frame of reference, the sum of all forces, including the advection tendencies, would sum to zero, because the flow is stationary. Streamwise force components in Fig. 5a show that the accelerating geostrophic pressure gradient force and the decelerating momentum flux divergence are in balance and sum to zero at the inflow. The streamwise component of the Coriolis force is 350 zero because this force acts always perpendicular to the flow. Inside the wind farms, the momentum flux divergence is positive and thus is an accelerating component. It is the dominant driving force because it is more than seven times greater than the geostrophic pressure gradient force.
An increasing perturbation pressure in front of the wind farm leads to a negative perturbation pressure gradient force and thus flow deceleration (often called blockage effect). However, inside the wind farms the perturbation pressure gradient force 355 is positive due to a favorable pressure gradient (decreasing pressure). In the near wake the momentum flux divergence is high and leads to a fast speed recovery. The momentum flux divergence decreases fast until it becomes negative in the far wake, so that the speed recovers slowly in the far wake. The only force that remains for driving the flow is the geostrophic pressure gradient force. At the inflow, this force is in balance with the momentum flux divergence, but in the wake this is not the case due to two reasons: first, the negative momentum flux divergence is weaker than at the inflow due to a smaller wind speed 360 and thus a reduced near-surface momentum flux. Second, the streamwise component of the geostrophic pressure gradient force has increased because the wake flow is deflected counterclockwise (i.e. to lower pressure). These results show that the wake deflection is an elementary feature of the wake that supports the wind speed recovery. They also show that mixing of momentum from the BL to the wind turbine level is not the dominant process that drives the speed recovery in the far wake of very large wind farms.

365
The wake deflection can be explained by examining the crosswise force components that are shown in Figure 5b perturbation pressure gradient force are negative inside the wind farm and inside the wake and are therefore opposing the wake deflection. The negative perturbation pressure gradient force is a result of the pressure distribution around the wind farms that is caused by the wind farm shape (cf. Fig.4d). The reason for the negative momentum flux divergence is the enhanced downward mixing of negative y-momentum of the overlying flow, which is veered to the right (cf. Fig. 3). For small wind farms this process can be dominant and may result in clockwise wake deflection (Van Der Laan and Nørmark Sørensen, 2017). However, 375 for very large wind farms, as in this study, the effect of the reduced Coriolis force is dominant. An appropriate parameter for estimating the importance of Coriolis effects is the Rossby number. Coriolis effects become dominant for Rossby numbers close to or below 1. For the large wind farms in Zone 3 with a length of L wf ≈ 100 km at mid-latitudes (Coriolis parameter f ≈ 10 −4 ) and a wind speed of U ≈ 10 m s −1 the Rossby number becomes:

380
indicating that Coriolis effects play an important role for flows in wind farms of this size.

Turbulence intensity at hub height
The turbulence intensity TI is defined as in Porté-Agel et al. (2013): where TKE is the resolved turbulence kinetic energy defined as: where u 2 , v 2 and w 2 are the resolved-scale variances of u, v and w, respectively.
The TI at hub height is shown in Fig. 6. Inside the wind farms, the TI reaches a fully developed state after approximately 4 rows and is constant farther downstream. A smaller turbine spacing leads to a greater TI inside the wind farms. For the case NBL-700-7D, a TI of 10 % is reached inside the wind farms, but more than 14 % is reached in the case NBL-700-5D. In 390 the CBL-cases, the TI inside the wind farms reaches 10 % in case CBL-700-7D and approximately 12 % in CBL-1400-7D.
Although the ambient TI is only approximately 3 % for the case SBL-300-7D, the TI inside the wind farms reaches values similar to NBL-700-7D (approximately 10 %).
The wake in terms of TI is generally shorter than the wake in terms of wind speed (cf. Fig. 4). The shortest wakes occur in the shallow SBL (SBL-300-7D) and the longest wakes occur in the thick CBL (CBL-1400-7D). This is the opposite behaviour 395 than that of the wake in terms of wind speed. The wake length in terms of TI weakly depends on wind farm size with slightly longer wakes for larger wind farms. However, this effect is caused by the definition of the TI, in which the wind speed variances are normalized by the mean horizontal wind speed. The mean horizontal wind speed is smaller in the wake of the large wind farms than in the wake of the small wind farms, resulting in a higher TI in the wake of the large wind farms. The wind farm size dependency of the wake length vanishes if the TKE is used for measuring the wake length instead of the TI (not shown).
In the NBL-cases and especially in the SBL-case, the TI in the far wake drops below the ambient TI at the inflow. This effect is caused by the reduced wind speed in the far wake which leads to a reduction in the shear-driven turbulence production. In the CBL-cases, there is also buoyancy-driven turbulence production, which is unaffected by the reduced wind speed in the wake and thus maintains the TI level. That buoyancy-driven turbulence production has a large impact on hub height TI is also verified by the fact that the ambient TI is greater in the CBL-case with the thick BL (CBL-1400-7D) than in the case with the 405 shallow BL (CBL-700-7D): The convective velocity scale w * is greater in the case CBL-1400-7D than in case CBL-700-7D and hence also the buoyancy-generated velocity variances are greater. Because buoyancy acts as a TKE-sink in the SBL-case, it can not compensate for the reduction in shear-driven turbulence production and thus the TI in the wake drops below 2 %.
This effect is amplified by the entrainment of warm air into the BL, that leads to a stabilization (increased lapse rate) at hub height and therefore stronger turbulence damping (cf. Fig. 6f and Fig. 8). The entrainment of warm air results in a temperature 410 increase at hub height of approximately 0.3 K for the large wind farms and approximately 0.1 K for the small wind farms.  Figure 7 shows vertical cross sections of the horizontal mean wind speed for all cases. The cross sections are located at y = 120 km and thus cross the large wind farms in Zone 3. The inversion layer height z i is marked by lines at which the maximum vertical potential temperature gradient occurs.

415
The inversion layer (IL) height is affected by the presence of the wind farms in all five cases. In the NBL-cases the IL is displaced upwards by 200 − 300 m, whereas a larger displacement occurs for the smaller turbine spacing. Because the mean wind speed inside the BL decreases in the streamwise direction, the IL must be displaced upwards in order to maintain a constant mass flux inside the BL. The increase in IL height is not caused by entrainment of warm air into the BL (as can also be seen in the profiles of potential temperature in Fig. 8). This phenomenon has also been observed by Allaerts and Meyers 420 (2017), who also stated that the mass flux conservation is the reason for the IL displacement. Abkar and Porté-Agel (2014) stated that a smaller turbine spacing results in a larger BL height for an infinite wind farm in a CNBL. The IL displacement causes an acceleration of the flow in the FA for the NBL and CBL cases-Details about this effect are described in the next section.
In the CBL-cases, the IL height increases above the wind farms and decrease above the wake, reaching its initial value at 425 approximately 70 km downstream of the last wind farm trailing edge. The IL displacement is larger for the shallower BL, i.e.
The IL displacement is most significant for the case SBL-300-7D. The IL height increases from 300 m to 500 m. Allaerts and Meyers (2016) also reported that larger IL displacements occur for shallower BLs (+60 % for h = 250 m). In the case SBL-300-7D, the IL height increase is caused by vertical displacement due to mass conservation and also by entrainment of 430 warm air into the BL (see Fig. 7f). The entrainment of warm air into the BL leads to a warming of the lower part of the BL.
However, the temperature at the height of the original IL is reduced because the warm air in the IL is replaced by relatively cold air from the BL.
Because the laminar flow in the FA is adiabatic, the isotherms in Fig. 7f can be interpreted as streamlines. They show that gravity waves are excited by the wind farms. There are small scale gravity waves with a wave length that corresponds to the 435 turbine spacing and a large scale gravity wave with a wave length that approximately corresponds to the wind farm length. The negative and positive temperature deviations in the wave crest and trough, respectively, cause a positive and negative deviation in the perturbation pressure at the surface, as it is shown in Fig. 4f. A detailed analysis of the wind farm induced gravity waves goes beyond the scope of this study. However, it is noted that the qualitative pressure and temperature distributions correspond to the findings of Allaerts and Meyers (2017) and Wu and Porté-Agel (2017). of Γ = 1 K km −1 (resulting in L z = 0.22λ z ) and Γ = 5 K km −1 (resulting in L z = 0.49λ z ). It is not clear whether the vertical wave length is the only relevant parameter for choosing the correct domain height or whether the wind farm length also has to be considered. Further research is needed to find setup guidelines that ensure that wind farm induced gravity waves are covered as realistically as possible.
3.1.5 Profiles of wind speed, wind direction and potential temperature in the wake

450
To examine the effect of the wind farms on the BL in more detail, profiles of wind speed, wind direction and potential temperature are shown in Fig. 8. The profiles are evaluated at the inflow (x = 0 km), in the near wake (x = 120 km) and in the far wake (x = 180 km) of the large wind farms in Zone 3 (y = 120 km).
The wind speed profiles show that the wind farm induced wind speed deficit spreads over the entire height of the BL. The effective vertical mixing in the CBL-cases results in an approximately height-constant wind speed at the inflow and in the 455 wake. In the case CBL-1400-7D, the wind speed in the upper part of the BL is even lower than in the lower part of the BL at x = 120 km. In the NBL-cases, the vertical mixing is not as effective and thus a significant wind shear exists over the entire BL in the wake. The wind speed profiles of the case SBL-300-7D show that the BL has grown from 300 m to 500 m and that the super-geostrophic maximum is eliminated completely. The IL displacement causes an increase in wind speed in the FA above the BL. The maximum increases (approximately 1 m s −1 ) are observed for the cases with the greatest IL displacements. That 460 suggests that the wind speed excess above the BL is also caused by the continuity constraint, i.e. the wind speed has to increase in order to maintain a constant mass flux between the IL and the domain top. For the CBL-cases, where the IL height decreases again behind the wind farms, the wind speed above the BL decreases to below-geostrophic in the far wake (x = 180 km). Note that these effects could be overestimated, because of the artificial boundary that is introduced by the Rayleigh damping layer that starts several hundred meters above the BL. The sensitivity of this effect on the Rayleigh damping height has not been 465 investigated because the scope of this study is on BL-internal effects.
The wake deflection shown in the horizontal cross sections can also be seen in the wind direction profiles. The wind farm induced wind direction change is approximately constant over the entire height of the BL. The largest deflection angles of up to −10 • are observed for the cases with the greatest speed deficit (SBL-300-7D and NBL-700-5D) because the Coriolis force reduction is greatest in these cases. The smallest deflection angle is observed in the case CBL-1400-7D with no deflection in 470 the upper half of the BL. All cases have in common that, in contrast to the counter clockwise deflection in the BL, the flow in the FA is deflected clockwise. As a result, the wind veer in the inversion layer increases to approximately 10 • . The flow deflection in the FA is also a Coriolis effect. Because the wind speed in the FA is supergeostrophic, the Coriolis force is greater than the geostrophic pressure gradient force and therefore the flow is deflected clockwise. The largest deflection angle of more than 10 • is observed for the case NBL-700-5D. For this case the wind speed excess in the FA is greatest and stays positive 475 until x = 180 km. Thus the Coriolis force is greatest and acts on the flow for a long time so that a high deflection angle is achieved. Note that this effect might also be overestimated due to the potentially overestimated wind speed excess. In the case SBL-300-7D, the combination of clockwise deflection in the FA and counterclockwise deflection in the BL leads results in a total wind veer of approximately 40 • between the surface and the FA. The effect of the wind farms on the potential temperature profiles is largest for shallow BLs (SBL-300-7D) and negligible 480 small for thick BLs (CBL-1400-7D). The potential temperature profiles inside the well-mixed BLs of the NBL-and CBL-cases are nearly unaffected by the wind farms. The greatest changes take place in the inversion layer, which is displaced upwards in order to maintain a constant mass flux in the BL, as already described in the previous section. The profiles show that the potential temperature inside the BL is unchanged and thus BL-warming due to entrainment of warm air from the FA is not the reason for the increased IL height. On the contrary, the temperature at the height of the original IL decreases by approximately 485 0.5 K because it is replaced by colder air from the underlying BL. The potential temperature profile of the SBL-case is heavily modified by the wind farms. The temperature in the BL increases by approximately 0.5 K due to entrainment of warm air from the FA into the BL. The IL rises from 300 m to 500 m due to the combined effect of BL warming and IL displacement. Because the surface temperature is constant, a new SBL forms in the far wake. This new SBL is shallower and more stably stratified than the original SBL at the inflow.
The reference power is obtained by an additional simulation of a single turbine using the same inflow profiles as for the case 505 NBL-700-7D. The reference power is P ref = 12.56 MW. The wind farm efficiency can also be interpreted as the wind turbine efficiency averaged over all wind turbines of the wind farm. The wind turbine efficiency of a wind turbine generating P wt is defined as: The wind farm efficiencies of the small and the large wind farm are listed in Table 2  height has also been observed by Allaerts and Meyers (2016), who reported a 17.6 % increase in power deficit for a BL height reduction from 1000 m to 250 m.

530
A comparison between the cases NBL-700-7D and CBL-700-7D shows that greater wind farm efficiencies are obtained for the NBL, although better efficiencies are expected for the CBL due to the better vertical mixing. The wind farm efficiencies for the NBL-case are greater because the inflow wind speed in the bulk of the BL is higher for the NBL-case than for the CBL-case (cf. Fig. 3). This result shows that it is important to consider not only the wind speed at hub height, but also the wind profile inside the entire BL to make accurate wind farm performance predictions.

Energy flux analysis
To examine the dependency of the wind farm efficiency on the turbine spacing and the BL height in more detail, an energy flux analysis is made in this section. The analysis is a simplified version of the analyses made by Abkar and Porté-Agel (2014) and Allaerts and Meyers (2017). The extraction of kinetic energy by the wind turbines can be compensated for by two sources of energy: The resolved downward turbulent flux of mean kinetic energy at rotor top level, averaged between y = 120 km − s y and y = 120 km + s y , is calculated by multiplying the shear stress with the corresponding wind velocity component at that height: 545 and the power density of the energy input by the geostrophic pressure gradient on the flow below rotor top level z t = 270 m is calculated as: The work done by the geostrophic pressure gradient on the rest of the BL (between z t and z i ) is calculated as:

550
In Fig. 10 these power densities are compared with the power densities of the wind turbines: For the case NBL-700-7D it can be seen that the first-row wind turbines operate at the reference power, so that a high power density of W wt = 12.56 MW/(7D) 2 = 4.45 W m −2 is achieved. The power density of downstream wind turbines is determined by the kinetic energy flux. Because the kinetic energy flux decays from 3 W m −2 at the beginning of the first wind farm 555 to 2 W m −2 at the end of the last wind farm, the power density of the wind turbines also decays to below 2 W m −2 . The good correlation between the wind turbine power density and the kinetic energy flux has also been found by Stevens et al. (2016) for the fully developed regime in a 9 km long wind farm. The work done by the pressure gradient on the flow below the rotor top level achieves a power density of approximately 0.6 W m −2 . It is thus not the dominating energy source inside the wind farms but it still contributes approximately 20 % to the total energy input W vkef + W gpg,wt . The wind turbines extract approximately 560 70 % of the total energy input W kef + W gpg,wt , which is a relatively large value. Johnstone and Coleman (2012)   shows (Fig. 10b), a reduction of the turbine spacing from s = 7 D to s = 5 D results in a doubling of the power density of the first-row wind turbines, but the power density of the last row wind turbines is as low as for s = 7 D. This result indicates that the turbine spacing for very large wind farms should be chosen much larger than for small wind farms to achieve a good wind farm efficiency. That the wind farm power output is limited by the vertical kinetic energy flux has also been found by Badger  Fig. 10f). Consequently, the wind farm induced flow effects result in a significant increase of the energy input by the pressure gradient, as it can be seen in Fig.10a and b. This effect occurs only in the wake although the BL height and wind direction already change inside the wind 585 farms. The reason is the decrease in the absolute wind speed that tends to reduce the ageostrophic wind speed component and thus compensates the increasing ratio of ageostrophic to geostrophic wind speed (counterclockwise wind direction change).
In the wake, the wind speed recovers and thus the ageostrophic wind speed component becomes larger than upstream of the wind farms. The described effect is largest for case with the small turbine spacing NBL-700-5D, because the BL growth and the wake deflection angle is largest for this case.

590
Figures 10c and d show that a doubling of the BL height has approximately no effect on the energy input by the pressure gradient on the undisturbed inflow. The effect of the thicker BL is compensated by a smaller ageostrophic wind speed component inside the BL. This is indicated by a much smaller angle between hub height wind and the geostrophic wind of α = 3.4 • for the case CBL-1400-7D than α = 9.5 • for the case CBL-700-7D (cf. Table 1 and Fig. 3). Consequently, the ageostrophic wind speed component inside the BL of the stationary inflow adjusts in such a way that the resulting energy input by the pressure 595 gradient balances the energy extraction by TKE production near the surface. As stated earlier, the power output of infinitely large wind farms is determined by the energy input of the pressure gradient. Hence, the power output of infinitely large wind farms does not depend on the BL height, at least for this idealized setups with a stationary CBL inflow. However, for very large, but finite-size wind farms, as in this study, a significant amount of the energy extracted by the wind turbines comes from the horizontal influx of mean kinetic energy inside the BL. As this influx is proportional to the BL height, much more energy is 600 available for the wind turbines in the case CBL-1400-7D than in the case CBL-700-7D. This results in a higher kinetic energy flux and also a slower decay of the vertical kinetic energy flux for the case CBL-1400-7D than for the case CBL-700-7D, as shown in Fig. 10c and d. The same holds for the wind turbine power densities, because they are directly related to the vertical kinetic energy flux. Therefore, higher wind farm efficiencies are achieved for the thicker BL.
The case SBL-300-7D is very special, because the rotor top level matches the BL height of the inflow. Thus, the energy 605 input by the pressure gradient above the rotor top level, as well as the vertical kinetic energy flux at the rotor top level, is zero which results in a smaller power density of the first-row wind turbines compared to the other 7D-cases and in a constant power density from approximately row 10. This redistribution is done by a favorable perturbation pressure gradient inside the wind farms and reaches power densities of approximately 1 W m −2 (not shown in Fig. 10). In the wake, the vertical kinetic energy flux at rotor top level drops to zero again, which is consistent with the very low TI in the wake (cf. Fig. 6).
These results show that the power output and the wake of very large wind farms behave very differently compared to small 615 wind farms. The main findings and their implications are summarized in the next section.

Conclusions
This study investigates wake properties and power output of very large wind farms with different turbine spacings in boundary layers (BL) of different stabilities and heights. Very large wind farms do not only change wind speed and turbulence intensity (TI) at wind turbine level but rather affect several flow quantities inside the entire BL and even above the BL. BL growth, 620 counterclockwise flow deflection inside the BL and clockwise flow deflection above the BL are the main effects that distinguish large from small wind farm flows. Wake lengths of very large wind farms are longer for shallower BLs and smaller turbine spacings, reaching values of more than 100 km. Thus, very large wind farms in the German Bight have the potential to affect the wind farm performance of neighbouring states such as Denmark or the Netherlands. The wake length in terms of TI is relatively independent of the wind farm size and is in general much smaller (approximately 20 km) than the wake length in 625 terms of speed deficit. Longer TI-wakes occur for convective BLs and shorter wakes for stable BLs due to the buoyancy-driven turbulence production or destruction.
For shallow, stable BLs very large wind farms trigger large scale gravity waves in the free atmosphere that cause significant flow blockage, affecting also smaller wind farms that are nearby. Some tuning of the domain height and the boundary conditions was necessary to capture this phenomenon correctly. Because shallow BLs occur quite frequently in the German Bight, it is an 630 important task to find best practice rules for simulation setups that capture this phenomenon as realistically as possible.
The wind speed recovery inside the wind farms is mainly driven by the turbulent vertical momentum flux but the wind speed recovery in the wake of very large wind farms is mainly driven by the geostrophic pressure gradient force. Thus, it is expected that the wake recovery of very large wind farms rather depends on the ageostrophic wind speed component than on parameters that affect the turbulent momentum flux such as stability or TI. Further investigations are needed to proof this hypothesis.

635
The power output of very large wind farms is limited by the available kinetic energy inside the BL and the energy input by the geostrophic pressure gradient. The achieved power density of turbines in the upstream part if the wind farms is significantly affected by the BL height, whereas the power density of the downstream turbines approaches the power density given by the energy input of the geostrophic pressure gradient. Because this power density is only as small as 2 W m −2 , the turbine spacing wake deflection towards lower pressure tend to increase the power input by the geostrophic pressure gradient, which could have a positive effect on the power output of downstream wind farms.
Overall, the results show that very large wind farms trigger much more complex flow effects than small wind farms do. It will be necessary to consider at least some of these effects in simple wake models in order to accurately predict the power output of very large wind farms. One of the next research tasks could be to derive empirical rules for predicting the power 645 output of very large wind farms by performing a more systematic and idealized set of simulations.
Code and data availability. The PALM code can be accessed under https://palm.muk.uni-hannover.de. PALM input files and plot scripts are available at https://doi.org/10.25835/0004522. Output data are available on request.