The actuator farm model for large eddy simulation (LES) of wind-farm-induced atmospheric gravity waves and farm–farm interaction

This study introduces the actuator farm model (AFM), a novel parameterization for simulating wind turbines within large eddy simulations (LESs) of wind farms. Unlike conventional models like the actuator disk (AD) or actuator line (AL), the AFM utilizes a single actuator point at the rotor center and only requires two to three mesh cells across the rotor diameter. Turbine force is distributed to the surrounding cells using a new projection function characterized by an axisymmetric spatial support in the rotor plane and Gaussian decay in the streamwise direction. The spatial support's size is controlled by three parameters: the half-decay radius $r_{\mathrm{1}/\mathrm{2}}$, smoothness s, and streamwise standard deviation σ. Numerical experiments on an isolated National Renewable Energy Laboratory (NREL) 5MW wind turbine demonstrate that selecting $r_{\mathrm{1}/\mathrm{2}}=R$ (where R is the turbine radius), s between 6 and 10, and $\mathit{\sigma }\approx \mathrm{\Delta }x/\mathrm{1.6}$ (where Δx is the grid size in the streamwise direction) yields wake deficit profiles, turbine thrust, and power predictions similar to those obtained using the actuator disk model (ADM), irrespective of horizontal grid spacing down to the order of the rotor radius.

Using these parameters, LESs of a small cluster of 25 turbines in both staggered and aligned layouts are conducted at different horizontal grid resolutions using the AFM. Results are compared against ADM simulations employing a spatial resolution that places at least 10 grid points across the rotor diameter. The wind farm is placed in a neutral atmospheric boundary layer (ABL) with turbulent inflow conditions interpolated from a previous simulation without turbines. At horizontal resolutions finer than or equal to $R/\mathrm{2}$, the AFM yields similar velocity, shear stress, turbine thrust, and power as the ADM. Coarser resolutions reveal the AFM's ability to accurately capture power at the non-waked wind farm rows, although it underestimates the power of waked turbines. However, the far wake of the cluster can be predicted well even when the cell size is of the order of the turbine radius.

Finally, combining the AFM with a domain nesting method allows us to conduct simulations of two aligned wind farms in a fully neutral ABL and of wind-farm-induced atmospheric gravity waves under a conventionally neutral ABL, obtaining excellent agreement with ADM simulations but with much lower computational cost. The simulations highlight the AFM's ability to investigate the mutual interactions between large turbine arrays and the thermally stratified atmosphere.

1 Introduction

The global offshore wind capacity has been expanding at a rapid rate in the past few years, with the global installed capacity estimated to reach 500 GW by the end of 2030 according to the 1.5 °C global temperature rise scenario (IRENA, 2023). This signifies a 14-fold increase compared to 2020 levels. As the number and size of offshore wind farms increase, they are often clustered to maximize the use of the available wind energy resources as well as to minimize the infrastructure costs (Akhtar et al., 2021; Junqueira et al., 2021). Wind farm clusters can lead to reduced wind farm power production due to the impact of wakes shed by upstream wind farms, which may persist for several kilometers downstream as reported from airborne (Platis et al., 2018; Lampert et al., 2020), lidar (Bodini et al., 2021; Schneemann et al., 2020), and satellite synthetic aperture radar (SAR) measurements (Hasager et al., 2015; Ahsbahs et al., 2020), as well as mesoscale numerical studies (Pryor and Barthelmie, 2024; Fischereit et al., 2022; Wang et al., 2023). The extent of wind farm wake propagation depends on factors such as atmospheric stability, turbulence intensity, boundary layer height, and wind farm size (Schneemann et al., 2020; Cañadillas et al., 2020; Pryor et al., 2021). Moreover, recent numerical studies (Stipa et al., 2024b; Maas and Raasch, 2022) have shown how wind-farm-induced atmospheric gravity waves can be triggered within a shallow boundary layer under stable free-atmosphere stratification, leading to horizontal pressure gradients that influence the flow in the region upstream of the farm – leading to what is commonly referred to as blockage – as well as the wind farm wake. This underlines the importance of studying the mutual interactions between neighbouring farms as well as the interaction of farm clusters with the atmospheric boundary layer (ABL) and the thermally stratified free atmosphere aloft.

Wind farm clusters have been studied using various numerical models differing in the breadth of the resolved temporal and spatial scales. These include engineering models; computational fluid dynamics (CFD) models, such as large eddy simulation (LES) or Reynolds-averaged Navier–Stokes (RANS); or high-fidelity mesoscale numerical models such as the Weather Research and Forecasting (WRF) model (Skamarock et al., 2019).

Engineering models usually combine different sub-models to capture various physical processes such as individual turbine wakes (Bastankhah and Porté-Agel, 2014; Blondel and Cathelain, 2020), their interaction and merging (Niayifar and Porté-Agel, 2016), and blockage effects (Branlard and Meyer Forsting, 2020). On the one hand, turbine wake models tend to underestimate wind farm wake loss (Fischereit et al., 2022; van der Laan et al., 2023b) due to the assumptions regarding the wake profile, wake expansion, and superposition of wakes. Additionally, Gayle Nygaard et al. (2020) showed that individual wake models strongly overestimate wind farm wake recovery when applied to large clusters. On the other hand, individual blockage models underestimate the full extent of the blockage effect as they do not consider the wind farm interaction with the atmosphere. In this regard, Stipa et al. (2024b) and Devesse et al. (2023) showed that the accuracy of these engineering models improves substantially when they are coupled with a reduced-order mesoscale model, making them more suitable when looking at large wind farm clusters.

High-fidelity mesoscale models like WRF utilize a coarse grid resolution such that the flow around individual wind turbines is not resolved and hence the wind farm drag force must be parameterized. The Fitch et al. (2012) wind farm parametrization model represents the effect of wind turbines as a momentum sink. Specifically, a portion of the flow kinetic energy is assumed to produce electricity, while the rest is converted into turbulent kinetic energy (TKE). Eriksson et al. (2015) simulated the Lillgrund wind farm, located off the coast of southern Sweden, and compared WRF's wind farm parametrizations against LES data, where wind turbines were modeled using the actuator disk method (Mikkelsen, 2004). In this analysis, the WRF model overestimated the wind farm power and predicted faster wake recovery compared to the LES. Studies by Vanderwende et al. (2016) and Peña et al. (2022) also indicate that adopting WRF's turbine parameterizations yields lower velocity deficits and higher TKE in the wind turbine wake than LES. This suggests that, when possible, mesoscale models require validation against microscale models, which are capable of resolving finer temporal and spatial scales.

In recent years, there has been a significant increase in studies employing LES to examine the flow around large wind farms and the evolution of the cluster wake (see for example Maas, 2023; Maas and Raasch, 2022; Cheung et al., 2023; Stipa et al., 2024a; Lanzilao and Meyers, 2022b; Stieren and Stevens, 2022, among others). This trend is supported by the growing availability of computational resources, allowing the temporal fluctuations and large-scale flow features of the ABL to be captured alongside the dynamics of wind turbine wakes, which are both equally important to accurately simulate the flow around large wind farm clusters. Nevertheless, this type of simulation still requires significant computational resources due to the size of the flow domain and grid resolution constraints. For instance, Stieren and Stevens (2022) used LES to study the impact of a rectangular wind farm wake on another wind farm located downstream in fully neutral ABL conditions, with a horizontal domain size of about 390 km² and 1800 × 480 × 480 mesh cells. To simulate real-world wind farm clusters, an even larger domain is usually required. Maas and Raasch (2022) performed a series of LESs of the wind farm clusters in the German Bight under different ABL stability conditions using a horizontal domain size of approximately 33 620 km², around 8.4×10⁹ mesh cells, and a finest grid spacing of 20 × 20 × 20 m. Cheung et al. (2023) simulated the 10 000 km² AWAKEN wind farm site in Oklahoma (USA) using 21.4×10⁹ mesh cells, a background mesh size of 20 × 20 × 20 m, and a finest grid spacing of 2.5 × 2.5 × 2.5 m. Additionally, the vertical extent of the domain should also be increased when simulating the effects of wind-farm-induced gravity waves within the free atmosphere due to the large vertical wavelength of these waves (Allaerts and Meyers, 2017).

Wind farm LES studies typically employ a precursor–successor method, where the precursor simulation is used to develop the ABL turbulent inflow which is subsequently utilized in the successor simulation, which includes the wind turbines. In the context of ABL LES, the grid resolution is dictated by the requirement to correctly capture the theoretical law-of-the-wall (LOTW) scaling, which imposes a specific cell aspect ratio at the wall, as well as a minimum number of cells within the boundary layer (Brasseur and Wei, 2010). Notably, the inclusion of wind turbines in the successor simulation introduces an additional grid resolution constraint, typically tighter than that required to capture LOTW scaling, which depends on how wind turbines are represented within the numerical domain. For example, the actuator line model (ALM; Sørensen and Shen, 2002) requires several mesh cells across the rotor radius (Martínez-Tossas et al., 2015a), and it is usually employed in simulations of isolated wind turbines or small clusters of two to three machines due to its additional requirement that the rotating blade tip cannot cross more than one mesh cell per time step, substantially limiting the time step size. For this reason, LES studies of large wind farms usually employ the actuator disk model (ADM; Mikkelsen, 2004). In fact, according to Wu and Porté-Agel (2013), the ADM requires at least seven grid points across the turbine diameter for sufficient spatial resolution, while it leaves the time step to be determined according to the flow solution. Other authors fit more grid points within one diameter; for example Cheung et al. (2023) used around 40, Maas and Raasch (2022) used 12, and Lanzilao and Meyers (2023) used 9. This restricts the grid size to the range of 10 to 20 m for most LES studies, depending on this specific choice and on the wind turbine diameter. Notably, those studies that aim at investigating the effect of wind direction changes (see for example Stieren et al., 2021) require maintaining this condition not only in the spanwise and vertical directions, but also in the streamwise direction, as wind turbines or the wind rotate dynamically.

Within this landscape, it is evident that alternate wind turbine models that overcome the grid resolution constraint introduced by existing actuator models will be beneficial in reducing the computational cost of conducting LESs of large wind farms, especially when looking at farm–farm interactions. Recently, van der Laan et al. (2023a) developed a RANS-based wind farm parametrization model similar to a forest canopy model that applies a wind farm drag obtained by filtering each wind turbine location using a two-dimensional Gaussian kernel. While the model compares well against RANS ADM simulations when looking at the entire cluster, it requires multiple RANS ADM simulations for every wind farm represented with the new parameterization in order to compute the wind-drag coefficient relation, required as a model input.

In the present work, a new actuator model, referred to as the actuator farm model (AFM), is developed and validated. In contrast to conventional actuator models, the AFM requires a single actuator point positioned at the rotor center and only three to four mesh cells across the rotor diameter. This essentially reduces the grid constraint only to that imposed by the simulation's ability to capture LOTW scaling. The turbine force is projected from the actuator point to the surrounding grid cells using a new projection function characterized by axisymmetric spatial support in the rotor plane and Gaussian decay in the streamwise direction. Although the AFM solution of a single turbine simulation approaches the ADM solution when a similar grid size is used, the application domain of the AFM is tailored to problems requiring large domains that would otherwise be too computationally intense to be carried out using the ADM, such as studies of cluster wake evolution, as well as farm–farm and farm–ABL interactions. In terms of model fidelity at the turbine scale, the AFM LES lies in between the more detailed ADM LES and the parameterizations employed in numerical weather prediction codes (e.g. Fitch et al., 2012, or simple canopy models). Two illustrative applications of the AFM are the investigation of wind-farm-induced atmospheric gravity waves and farm–farm interaction via a hybrid AFM–ADM approach, whereby an upstream wind farm is modeled using the AFM with a coarser grid, while a downwind farm of interest is represented using the ADM within a nested finer grid.

The present work is organized as follows. Section 2 introduces the classic non-rotating ADM, the novel AFM formulation, and the grid nesting method used in the present study. Section 3 presents the parametric analysis conducted on an isolated wind turbine to choose the best set of AFM parameters, as well as to investigate the model's sensitivity to the grid size. Section 4 describes the AFM simulations performed on both an aligned and a staggered wind farm layout, showing their comparison against ADM results. In Sects. 5 and 6, the AFM combined with grid nesting is leveraged to study the interaction between two aligned wind farms and to simulate wind-farm-induced gravity waves. Moreover, challenges posed by coarsening the grid in the context of wind farm LES are discussed in Appendix A, and a wall model correction is therefore proposed. In Appendix B, the effect of such correction is analyzed, and precursor simulations used to provide a time-resolved inlet boundary condition to the wind farm analyses of Sects. 4 and 5 are described. Finally, Sect. 7 highlights the conclusions of the present study.

2 Methodology

For the LESs presented in this paper, we use the open-source finite-volume code TOSCA (Toolbox fOr Stratified Convective Atmospheres) developed at the University of British Columbia, Canada. The governing equations correspond to mass and momentum conservation, while sub-grid-scale stresses are calculated with the model proposed by Meneveau et al. (1996), where the dynamic Smagorinsky model coefficient is averaged along the flow pathlines in a Lagrangian sense. A potential temperature transport equation can also be solved to account for stability effects inside the ABL and in the free atmosphere. As specific details about TOSCA's numerical method are provided in Stipa et al. (2024b), together with an exhaustive validation, the present section solely focuses on new features that are relevant for the current paper. Specifically, the conventional non-rotating uniform ADM is described in Sect. 2.1, while Sect. 2.2 presents the developed AFM. Lastly, Sect. 2.3 describes the domain nesting technique used in TOSCA, referred to as overset mesh, which will be used for the simulations described in Sects. 5 and 6.

2.1 Uniform actuator disk model

The ADM represents each individual wind turbine within the computational domain by discretizing each rotor disk with a certain number of actuator points in the radial and azimuthal directions. In principle, the ADM allows a radially varying blade force to be modeled upon knowing the blade chord, twist, and type of airfoil at each radial location and by providing lift and drag look-up tables for each airfoil (Martínez-Tossas et al., 2015a). In this case, both wind turbine thrust and torque are modeled and the wind turbine can be additionally equipped with angular velocity and pitch controllers. For this class of the ADM, thrust and power coefficients are computed variables that vary in time, making it more difficult to compare their results against other classes of the ADM where the turbine thrust coefficient is a model input (the so-called “uniform” ADM; Porté-Agel et al., 2010; Calaf et al., 2010; Jimenez et al., 2007, 2008). As a consequence, in order to exactly match the turbine thrust coefficient applied in the ADM to that of the AFM when comparing the two models, we employ the uniform ADM throughout the entire paper. This class of ADM only applies a thrust force to the flow, making it unable to capture the tangential force exerted on the flow by the wind turbine. As a consequence, uniform ADMs are only expected to produce accurate results in the far region of the wind turbine wake. For brevity, the modifier uniform will be omitted throughout the remainder of the paper.

In this type of model, wind turbine thrust at each actuator point is calculated as

$$\begin{array}{}\text{(1)}& {\mathbf{T}}_p=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left|\right|{\mathbf{U}}_{\mathrm{d}}\left|\right|{\mathbf{U}}_{\mathrm{d}}{C}_{\mathrm{T}}^{\prime }{\mathrm{dA}}_{\mathrm{p}},\end{array}$$

where dA_p is the portion of rotor disk area associated with each actuator point and U_d is the disk velocity. The latter is obtained by first interpolating the flow velocity from the background mesh at the actuator points and then averaging among all actuator points. The disk-based thrust coefficient $C_{\mathrm{T}}^{\prime }$ can be determined from the thrust coefficient C_T using the relation

$$\begin{array}{}\text{(2)}& C_{\mathrm{T}}^{\prime }={\displaystyle \frac{C_{\mathrm{T}}}{(\mathrm{1}-a{)}^{\mathrm{2}}}},\end{array}$$

where a is the turbine axial induction factor (Calaf et al., 2010). The point force at each actuator point given by Eq. (1) should be projected to the mesh cells and transformed into a body force field that represents the effect of the wind turbine on the incoming flow. This operation is performed by means of the projection kernel g_AD(x), defined as

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{g}_{\mathrm{AD}}\left(\mathit{x}\right)& ={\displaystyle \frac{\mathrm{1}}{{\mathit{ϵ}}^{\mathrm{3}}{\mathit{\pi }}^{\mathrm{3}/\mathrm{2}}}}\exp (-{\displaystyle \frac{(x_c-x_p{)}^{\mathrm{2}}}{{\mathit{ϵ}}^{\mathrm{2}}}}-{\displaystyle \frac{(y_c-y_p{)}^{\mathrm{2}}}{{\mathit{ϵ}}^{\mathrm{2}}}}\\ & -{\displaystyle \frac{(z_c-z_p{)}^{\mathrm{2}}}{{\mathit{ϵ}}^{\mathrm{2}}}}),\end{array}\end{array}$$

where subscripts c and p refer to the mesh cell center and actuator point, respectively. The quantity ϵ is the projection width, and it should be set to 1.5–2 times the mesh size along the rotor plane (Calaf et al., 2010; Martínez-Tossas et al., 2015b). Using Eq. (3), the wind turbine body force at each mesh cell can be evaluated as

$$\begin{array}{}\text{(4)}& {\mathit{b}}_c=\sum _{p=\mathrm{1}}^{N_p}{g}_{\mathrm{AD}}{\mathbf{T}}_p,\end{array}$$

where N_p is the total number of actuator points. Note that the body force contribution at any given cell c may come from different actuator points. Moreover, the sum of T_p among all actuator points yields a force that is equivalent to the total wind turbine thrust but directed in the opposite direction.

2.2 Actuator farm model

The AFM is based on a similar concept to that of the ADM, but, instead of representing each wind turbine as a distribution of actuator points, a single point is used, located at the rotor center. Hence, the thrust force at this single actuator point coincides with the total wind turbine thrust, and it is evaluated as

$$\begin{array}{}\text{(5)}& \mathbf{T}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left|\right|{\mathbf{U}}_{\mathrm{d}}\left|\right|{\mathbf{U}}_{\mathrm{d}}{C}_{\mathrm{T}}^{\prime }\mathit{\pi }{R}^{\mathrm{2}},\end{array}$$

where R is the turbine radius. The main implication of using a single actuator point is that the spatial support of the body force field, given by the projection kernel's convolution with the location of the actuator points in the ADM, is equivalent to the projection kernel itself in the AFM. For this reason, the spatial support of the projection function should be of the order of the rotor disk. One option is to use the Gaussian kernel expressed by Eq. (3), with σ≥R, to distribute the turbine thrust to the neighbouring mesh cells. However, this approach produces a body force that is too concentrated close to the rotor center (not shown here), resulting in an artificial shedding of the wind turbine wake for high thrust coefficients. For this reason, a new projection function has been developed which emulates the one resulting from the convolution of Eq. (3) with the location of the ADM points. In addition, any kernel used within a generic actuator model should integrate to unity over the volume so that the total wind turbine thrust is recovered upon integrating the body force.

First, we define a Cartesian coordinate system 𝒞, having its origin at the rotor center, the x axis aligned with the wind direction, the z axis pointing in the vertical direction, and the y axis forming a right-handed coordinate frame. Within 𝒞, we can define a function f(x) as

$$\begin{array}{}\text{(6)}& f\left(\mathit{x}\right)={\displaystyle \frac{\mathrm{1}}{\exp \left(\frac{\sqrt{y^{\mathrm{2}}+z^{\mathrm{2}}}-r_{\mathrm{1}/\mathrm{2}}}{s}\right)+\mathrm{1}}}\exp \left(-{\displaystyle \frac{x^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}\right),\end{array}$$

where $r_{\mathrm{1}/\mathrm{2}}$ is the radius in the (y,z) plane where the function's magnitude decays by $\mathrm{1}/\mathrm{2}$, s is a smoothing parameter, and σ is the standard deviation of the Gaussian decay along x. Equation (6) is axisymmetric on the rotor plane with respect to the rotor center, and it coincides with a Gaussian function in the x direction. To be used as a projection kernel, it must be divided by its integral over the volume so that the integral of the resulting function equals unity. In order to more easily perform the integral, Eq. (6) is expressed in cylindrical coordinates. As a result, expressing the differential volume given by dxdydz as rdrdθdx, the definite integral of Eq. (6) over the entire domain can be written as

$$\begin{array}{}\text{(7)}& \begin{array}{rl}\mathcal{I}& =\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\mathrm{d}\mathit{\theta }\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\exp \left(-{\displaystyle \frac{x^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}\right)\mathrm{d}x\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{r}{\exp \left(\frac{r-r_{\mathrm{1}/\mathrm{2}}}{s}\right)+\mathrm{1}}}\mathrm{d}r\\ & =\mathrm{2}\mathit{\sigma }{\mathit{\pi }}^{\mathrm{3}/\mathrm{2}}\left[s^{\mathrm{2}}{\mathrm{Li}}_{\mathrm{2}}\left(-\exp \left({\displaystyle \frac{r_{\mathrm{1}/\mathrm{2}}-r}{s}}\right)\right)\right.\\ & {\left.-sr\log \left(\exp \left({\displaystyle \frac{r_{\mathrm{1}/\mathrm{2}}-r}{s}}\right)+\mathrm{1}\right)\right]}_{\mathrm{0}}^{\mathrm{\infty }},\end{array}\end{array}$$

where Li₂ is the poly-logarithmic function of the second kind and the expression between square brackets is the integral of the portion of Eq. (6) which depends on r. Its values at zero and infinity are evaluated as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle }\begin{array}{rl}& {\left[s^{\mathrm{2}}{\mathrm{Li}}_{\mathrm{2}}\left(-\exp \left({\displaystyle \frac{r_{\mathrm{1}/\mathrm{2}}-r}{s}}\right)\right)-sr\log \left(\exp \left({\displaystyle \frac{r_{\mathrm{1}/\mathrm{2}}-r}{s}}\right)+\mathrm{1}\right)\right]}_{\mathrm{0}}\\ & =s^{\mathrm{2}}{\mathrm{Li}}_{\mathrm{2}}\left(-\exp \left({\displaystyle \frac{r_{\mathrm{1}/\mathrm{2}}}{s}}\right)\right),\end{array}\text{(9)}& {\displaystyle }\begin{array}{rl}& \underset{r\to \mathrm{\infty }}{lim}\left[s^{\mathrm{2}}{\mathrm{Li}}_{\mathrm{2}}\left(-\exp \left({\displaystyle \frac{r_{\mathrm{1}/\mathrm{2}}-r}{s}}\right)\right)-sr\log \left(\exp \left({\displaystyle \frac{r_{\mathrm{1}/\mathrm{2}}-r}{s}}\right)+\mathrm{1}\right)\right]\\ & =\mathrm{0}.\end{array}\end{array}$$

Hence, the functional form of the projection function which integrates to unity is finally given by

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{g}_{\mathrm{AF}}(x,r)& =-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\sigma }{\mathit{\pi }}^{\mathrm{3}/\mathrm{2}}{s}^{\mathrm{2}}{\mathrm{Li}}_{\mathrm{2}}\left[-\exp \left(\frac{r_{\mathrm{1}/\mathrm{2}}}{s}\right)\right]}}\\ & \cdot \left({\displaystyle \frac{\exp \left(-\frac{x^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}\right)}{\exp \left(\frac{r-r_{\mathrm{1}/\mathrm{2}}}{s}\right)+\mathrm{1}}}\right),\end{array}\end{array}$$

where a minus sign has been applied to Eq. (8) due to the integration in Eq. (7). Using Eq. (10), the wind turbine body force at each mesh cell can be evaluated as b_c=g_AFT. Notably, Eq. (10) has two free parameters, namely the half-decay radius $r_{\mathrm{1}/\mathrm{2}}$ and the smoothing s, while the streamwise standard deviation σ can be chosen following the same approach as the ADM.

As an illustrative example, the wind turbine force calculated using U_d=9 m s⁻¹, R=60 m, and $C_{\mathrm{T}}^{\prime }=\mathrm{0.7}$ is shown in Fig. 1 for both the ADM and AFM on a vertical plane through the turbine rotor center. The grid size along y and z is set to 5×5 m; σ is set to 20 m for both Eqs. (3) and (10); and $r_{\mathrm{1}/\mathrm{2}}$ and s are set to 60 m and 6, respectively. As can be appreciated, the definition of g_AF allows for the recovery of a body force field that is similar to that obtained from the ADM. However, this is achieved in the AFM using a single function that projects from the wind turbine rotor center instead of the summation of many different Gaussian functions centered at each actuator point. This property allows the AFM to require fewer mesh cells along the wind turbine radius to properly resolve the projection function in space.

The developed AFM projection function is limited by the wind turbine hub height to avoid errors in the body force owing to some of the force being projected outside of the domain. In practice, the $r_{\mathrm{1}/\mathrm{2}}$ parameter should be of the order of the wind turbine radius. Figure 2 shows the radially dependent component of Eq. (10) for two values of $r_{\mathrm{1}/\mathrm{2}}$ and different values of smoothing parameter s. The vertical dashed line indicates the wall, located at $h_{\mathrm{hub}}/R=\mathrm{1.5}$ from the rotor center. In order not to incur any projection error, Eq. (10) should reach a value close to zero before reaching the edge of the domain. For this reason, high values of s are not ideal, as they would yield a defect in the recovered wind turbine force following integration over the domain.

Figure 1 Comparison of two sections located at x=0 of the body force projection function for the ADM (a) and AFM (b). In the ADM, the force associated with each actuator point is equal to the total wind turbine force scaled by the ratio between the actuator element area and the rotor swept area. For the AFM, the force is projected from a single actuator point, located at the rotor center. The continuous black circle indicates the radial coordinate corresponding to the wind turbine radius. The parameter σ is set to 20 m for both Eqs. (3) and (10), while $r_{\mathrm{1}/\mathrm{2}}$ and s are set to 60 m and 6, respectively.

Figure 2 Dependency of Eq. (10) on the parameters $r_{\mathrm{1}/\mathrm{2}}$ (a: $r_{\mathrm{1}/\mathrm{2}}/R=\mathrm{1}$; b: $r_{\mathrm{1}/\mathrm{2}}/R=\mathrm{1.3}$) and on the smoothing s. The vertical dashed line indicates the wall, located at a distance h_hub from the rotor center ($h_{\text{hub}}/R=\mathrm{1.5}$ in the figure, but this ratio depends on the specific wind turbine under study).

Regarding the calculation of the disk velocity U_d for the AFM, two strategies have been tested within the present paper. A first method, referred to as the rotor disk sampling, consists of tri-linearly interpolating the wind speed at the rotor center from the eight cells surrounding the actuator point. Notably, the same approach is used in the ADM to find the wind speed at each actuator point before averaging. In a second strategy, inspired by Churchfield et al. (2017) and referred to as the integral sampling, the disk velocity is calculated as

$$\begin{array}{}\text{(11)}& {\mathbf{U}}_{\mathrm{d}}=\sum _{c=\mathrm{1}}^{N_{\mathrm{c}}}{g}_{\mathrm{AF}}{\mathit{u}}_c,\end{array}$$

where u_c is the wind speed at cell c and N_c is the number of cells contained in a sphere of radius $\mathrm{2}{r}_{\mathrm{1}/\mathrm{2}}$ from the rotor center. Notably, this radius is only defined for implementation purposes, and increasing it further has no effect, as g_AF decays to zero well before $\mathrm{2}{r}_{\mathrm{1}/\mathrm{2}}$.

One of the main benefits of the AFM is to relax the requirement imposed by the ADM, which dictates that around 10 mesh cells should be used along the rotor diameter. This allows computational resources to be saved by reducing the number of cells in the domain, especially for those wind farm simulations characterized by a domain that extends for tens of kilometers in each direction. However, while these large simulations are the main target for the AFM, it is crucial to understand the limits and implications of reducing the LES spatial resolution below what is commonly employed by the research community. In fact, grid coarsening may cause problems that are closely related to how the inlet boundary condition is prescribed within ABL LESs and to the simulation's ability to satisfy the theoretical law of the wall (LOTW) in the background atmospheric flow when a wall model is used (Brasseur and Wei, 2010). Due to its relevance in relation to using the AFM, this topic is expanded in detail in Appendix A. Here we provide a correction to commonly adopted wall models based on the classic Monin and Obukhov (1954) similarity theory, together with best-practices for the interpolation of time-resolved inlet boundary conditions in the context of precursor–successor simulations.

2.3 Overset mesh

The overset mesh technique (Benek et al., 1983) facilitates grid generation for flow around complex geometries as well as for bodies under relative motion (Meakin, 1983). It involves decomposition of the overall domain into a set of overlapping subdomains, such that the governing equations are independently solved in each of these domains and the information is transferred from one domain to the other at the subdomain interface using an interpolation method. Within wind energy applications, the overset mesh method has been used for blade-resolved simulations of wind turbines (Kirby et al., 2019) as well as simulations capturing synoptic (≈2000 km), meso-scale (≈100 km), and micro-scale (≈100 m) effects simultaneously via nested static grids with both one-way and two-way coupling (Liu et al., 2011; Mirocha et al., 2013). TOSCA applies a one-way coupling strategy where the information from the background (coarser) grid is transferred to the enclosed overset (finer) grid at the overset mesh boundaries, but no feedback is provided from the overset grid to the background grid. The AFM combined with an overset mesh method enables the application of a finer mesh grid in a region of interest, such as a downstream wind farm in which the wind turbines are modeled using the ADM.

The interpolation from the background mesh to the overset mesh begins with the identification of the face center points of the overset mesh cells along the interface, referred to as acceptor points. Then, for every acceptor point, the closest background mesh cell, referred to as donor cell, is identified. Using the relative position of the closest donor cell and the acceptor point, the eight donor cells from the background mesh that enclose an acceptor point are found and a tri-linear interpolation method is used to compute the velocity at the acceptor point. As the governing equations in TOSCA are formulated in generalized curvilinear coordinates, the numerical method solves for the contravariant fluxes instead of the Cartesian velocity (see Stipa et al., 2024b, for further details); hence, from the interpolated velocity at the acceptor points, the contravariant boundary fluxes are calculated at each iteration to obtain the boundary information for the nested inner domain.

The tri-linear interpolation scheme is non-conservative; hence, a correction of the local flux at the interface based on the mass residual is necessary to ensure global mass conservation. Here, the interpolated flux correction is made proportional to the flux, similar to Zang and Street (1995), for global conservation of mass. A local flux-proportional correction ${\stackrel{\mathrm{‾}}{U}}^r$ for the local flux U^r is given as

$$\begin{array}{}\text{(12)}& {\displaystyle }{\stackrel{\mathrm{‾}}{U}}^r=U^r-{\displaystyle \frac{{\mathit{ϵ}}_v\left|U^r\right|}{|U^r{|}_{\text{sum}}}}{\displaystyle \frac{\mathit{n}\cdot {\mathit{\zeta }}_r}{|\mathit{n}\cdot {\mathit{\zeta }}_r|}},\end{array}$$

where ϵ_v is the global mass flux imbalance, $|U^r{|}_{\text{sum}}$ is the sum of the flux magnitude, n is the outward pointing unit normal to the overset boundary and ζ_r is the unit normal to the curvilinear co-ordinate line. Global mass flux imbalance ϵ_v summed over the interface cell faces is given by

$$\begin{array}{}\text{(13)}& {\displaystyle }{\mathit{ϵ}}_v=\sum \left(U^r{\displaystyle \frac{\mathit{n}\cdot {\mathit{\zeta }}_r}{|\mathit{n}\cdot {\mathit{\zeta }}_r|}}\right).\end{array}$$

With respect to the wind farm simulations presented in Sects. 5 and 6, the above interpolation method is used for the overset domain at the streamwise inlet, spanwise lateral boundaries, and the upper boundary. The streamwise outlet employs a zero normal gradient on velocity, while the wall model defined by Eqs. (A6), (A7), and (A8) is used at the bottom wall.

3 Isolated wind turbine

In order to gain crucial insight regarding the optimal choice of the AFM settings, a parametric analysis is conducted by varying $r_{\mathrm{1}/\mathrm{2}}$, s, the velocity sampling strategy, and the horizontal mesh resolution for the uniform flow around an isolated wind turbine. Results are compared against two uniform ADM simulations characterized by a fine and a coarse mesh resolution. It should be kept in mind that the AFM, as the name suggests, is developed to model wind farm clusters rather than isolated turbines. However, the useful knowledge obtained by conducting these idealized isolated wind turbine simulations will be later applied to the wind farm studies presented in Sects. 4, 5, and 6.

Regarding the two ADM simulations used for comparison, the coarser one employs a grid resolution of 30 × 12.5 × 10 m in the streamwise, spanwise, and vertical directions, respectively. In the finer ADM simulation, this initial mesh is gradually refined in all directions to reach a uniform resolution of 2.1 m around the wind turbine. This fine region where the mesh is uniform extends one rotor diameter upstream of the turbine and five rotor diameters downstream. In the vertical direction, it extends ±h_hub above and below the hub height. Notably, the resolution of 12.5 m along y in the coarse ADM case is motivated by the fact that this is close to the largest lateral cell size that allows the National Renewable Energy Laboratory (NREL) 5MW wind turbine to be modeled using the ADM, as 10 mesh cells are placed along the rotor diameter. The streamwise and vertical grid spacings are similar to previous wind farm ABL studies (see Stipa et al., 2024b; Wu and Porté-Agel, 2017, among others). Following Calaf et al. (2010), the ratio $\mathit{\sigma }/\text{max}(\mathrm{\Delta }y,\mathrm{\Delta }z)$ in the ADM is set to 1.5 in order to avoid numerical oscillations when projecting the force from the actuator points. This leads to $\mathit{\sigma }/\mathrm{\Delta }x=\mathrm{0.625}$, which is then extended to all AFM simulations since Eq. (10) has a Gaussian shape in the streamwise direction. For the fine ADM case, where the mesh is uniform around the wind turbine, σ is chosen such that $\mathit{\sigma }/\mathrm{\Delta }=\mathrm{2}$. The projection error for the two ADM cases, evaluated as the relative difference between the cell-integrated body force after projection and the force sum from all actuator points, is equal to 2.5 % and 0.32 % for the coarse and fine cases, respectively. Throughout the paper the wind turbine corresponds to the NREL 5MW reference turbine (Jonkman et al., 2009), characterized by a radius R of 63 m and a hub height h_hub of 90 m.

Regarding the choice of the AFM parameters, the ratio $r_{\mathrm{1}/\mathrm{2}}/R$ is set to 0.8, 1, and 1.2, while the smoothing s is set to 2, 6, and 10. For each combination of these parameters, two types of velocity sampling methods are tested, namely the rotor disk and the integral sampling methods. Finally, four different mesh resolutions are chosen in the spanwise direction, namely 12.5, 20, 40, and 60 m. In the streamwise direction, previous studies showed evidence that an accurate solution can be obtained with a mesh resolution as large as 30 m, so the latter is used in conjunction with the values of 12.5 and 20 m along y. For the 40 and 60 m grids, mesh cells are rendered equal also in the streamwise direction. In the vertical direction, all cases feature a mesh resolution of 10 m, which corresponds to a representative value for wind farm ABL LESs. For all cases, periodic boundary conditions are applied at the spanwise boundaries, while a slip condition is enforced at the upper and lower boundaries. At the outlet, a zero normal gradient on velocity outflow is specified, while the inlet is set to a uniform velocity of 9 m s⁻¹. The turbine rotor is 600 m away from all boundaries except from the lower one, which is located at a distance equal to h_hub. This forces representative values of $r_{\mathrm{1}/\mathrm{2}}$ and f to be used to ensure that the projection function decays to zero before reaching the ground. As all four mesh configurations are uniform, they lead to the number of cells for each simulation reported in Table 1. All cases are advanced in time for 500 s, and data are averaged over the last 250 s. In total, this isolated turbine parametric study involves 72 simulations.

Figure 3 Velocity magnitude on a horizontal plane passing through the hub height for (a) ADM case with $\mathrm{2.1}\times \mathrm{2.1}\times \mathrm{2.1}$ m grid resolution, (b) ADM case with $\mathrm{30}\times \mathrm{12.5}\times \mathrm{10}$ m grid resolution, (c) AFM case with $\mathrm{30}\times \mathrm{12.5}\times \mathrm{10}$ m grid resolution, (d) AFM case with $\mathrm{30}\times \mathrm{20}\times \mathrm{10}$ m grid resolution, (e) AFM case with $\mathrm{40}\times \mathrm{40}\times \mathrm{10}$ m grid resolution, and (f) AFM case with $\mathrm{60}\times \mathrm{60}\times \mathrm{10}$ m grid resolution.

The velocity field resulting from the two ADM cases and from four AFM cases characterized by the same AFM settings ($r_{\mathrm{1}/\mathrm{2}}=R$ and s=6) and different mesh resolution is qualitatively shown in Fig. 3. As can be noticed, the ADM and AFM models predict a very similar velocity field when the same mesh is employed. The ADM model predicts a slightly higher velocity deficit than the AFM model (see the remainder of this section for a quantitative comparison), which is due to a spuriously increased turbine radius when accumulating the body force over all actuator points and to an increased body force towards the rotor center due to the body force being accumulated also from neighbouring points. This does not occur for the finer ADM case, where the standard deviation of the Gaussian projection function is reduced from 18.75 to 4.1 m. In fact, body force accumulation from neighbouring points is minimal in this case, and the actual wind turbine diameter is well represented. The coarsest AFM case is shown to visualize the flow field when the mesh is drastically coarsened. For this case, we also investigated the dependency of the AFM results to the relative position of the rotor disk with respect to the surrounding mesh cells (see Appendix C) and found this never produces a variation in thrust and power above 5 %, a number that is expected to further decrease for finer values of grid resolution.

Table 1 Number of mesh cells and total number of degrees of freedom (N_dofs) for each mesh configuration of the isolated wind turbine cases. The first line of the table corresponds to the fine ADM case, where the mesh is graded in each direction to reach a resolution of 2.1 m around the wind turbine. In all other cases the cell size is constant in each direction.

Figure 4 In (a)–(d), AFM to ADM thrust ratio $T_{\mathrm{AFM}}/T_{\mathrm{ADM}}$ for all AFM cases performed in the present section. Blue and black colours indicate rotor disk (RD) and integral (I) sampling, respectively, while different symbols indicate different values of $r_{\mathrm{1}/\mathrm{2}}/R$. The smoothness parameter is reported on the x axis, and data in the same panel are obtained using the same mesh resolution. Panels (e) and (f) report AFM to ADM thrust ratio $T_{\mathrm{AFM}}/T_{\mathrm{ADM}}$ for all cases performed in the present section, where (e) and (f) correspond to the integral and rotor disk sampling methods, respectively. The line colour identifies the mesh resolution, while each symbol corresponds to a different value of the smoothness parameter. The $r_{\mathrm{1}/\mathrm{2}}/R$ ratio is shown on the x axis.

Figure 4 shows the metric $T_{\mathrm{AFM}}/T_{\mathrm{ADM}}$, where T_ADM and T_AFM are the turbine thrust obtained with the ADM and AFM models, respectively, for all the cases conducted in this section. The ratio $T_{\mathrm{AFM}}/T_{\mathrm{ADM}}$ is not very sensitive to both the smoothing parameter and the mesh resolution, especially for the rotor disk sampling method. Conversely, the results seems to be greatly affected by the $r_{\mathrm{1}/\mathrm{2}}/R$ ratio, with ratios lower and higher than unity underestimating and overestimating turbine thrust, respectively. This behaviour is confirmed by looking at the two rightmost panels, where each panel contains all cases characterized by the same sampling method. For each value of $r_{\mathrm{1}/\mathrm{2}}/R$, the integral sampling method is more sensitive to the smoothing parameter and mesh resolution than the rotor disk sampling, which is expected because the sampled velocity is directly related to the spatial support of the projection function. The rotor disk sampling shows very little spread for any given value of $r_{\mathrm{1}/\mathrm{2}}/R$. For both sampling methods, $r_{\mathrm{1}/\mathrm{2}}/R=\mathrm{1}$ appears to provide the least error for wind turbine thrust.

Figure 5 shows the vertical profiles of velocity magnitude at different streamwise locations and at a spanwise coordinate coincident with the rotor center. The mesh resolution is characterized by different symbols, while the velocity sampling method is identified by their colour. All data correspond to s=6, and each panel refers to a different value of $r_{\mathrm{1}/\mathrm{2}}/R$. As noticed previously, AFM results are very sensitive to the value of $r_{\mathrm{1}/\mathrm{2}}/R$. When $r_{\mathrm{1}/\mathrm{2}}/R$ is low, the same turbine force has to be distributed over a smaller volume, thus increasing the body force locally. Conversely, when $r_{\mathrm{1}/\mathrm{2}}/R$ is large, the body force decreases as the force is projected over a larger volume. As a result, setting $r_{\mathrm{1}/\mathrm{2}}/R=\mathrm{1}$ represents the best choice to capture the velocity field around the wind turbine for all the investigated values of grid spacing and velocity sampling methods.

In Fig. 6, the same analysis is performed by fixing $r_{\mathrm{1}/\mathrm{2}}/R=\mathrm{1}$ and studying the dependence of the velocity profile on the smoothing parameter s. For $\left|y\right|<R$, varying the smoothness leads to a slight overestimation of the wake deficit for low s (more so when using the rotor disk sampling method) and an underestimation for high values of s. The opposite behaviour can be observed for $\left|y\right|>R$. This behaviour is expected, as increasing s increases the spatial support of the projection function, thus increasing the apparent radius of the wind turbine. Although the variation with s is small, s=6 seems to produce the best match in terms of velocity deficit when this is compared against ADM simulations, regardless of the mesh resolution. Moreover, some differences can be observed between the fine and coarse ADM simulations, where the smoothing generated by increasing the standard deviation of the Gaussian projection function leads to a wake deficit overestimation due to an increased body force towards the rotor center and to a smearing of the velocity profile when transitioning from the wake to the outer flow.

Figure 5 Vertical velocity profile at y=0 for s=6 and (a) $r_{\mathrm{1}/\mathrm{2}}/R=\mathrm{0.8}$, (b) $r_{\mathrm{1}/\mathrm{2}}/R=\mathrm{1}$, and (c) $r_{\mathrm{1}/\mathrm{2}}/R=\mathrm{1.2}$. Symbols indicate different mesh resolutions, and colours refer to the velocity sampling strategy. Continuous red and green lines indicate the fine and coarse ADM simulations, respectively. Each sub-panel corresponds to a different streamwise location, marked in the figure as a multiple of the turbine rotor diameter.

Figure 6 Vertical velocity profile at y=0 for $r_{\mathrm{1}/\mathrm{2}}/R=\mathrm{1}$ and (a) s=2, (b) s=6, and (c) s=10. Symbols indicate different mesh resolutions, and colours refer to the velocity sampling strategy. Continuous red and green lines indicate the fine and coarse ADM simulations, respectively. Each sub-panel corresponds to a different streamwise location, marked in the figure as a multiple of the turbine rotor diameter.

4 Isolated wind farm

This section describes the wind farm simulations' setup and results. Both an aligned wind farm and a staggered wind farm consisting of 25 wind turbines, organized in five rows and five columns, are investigated. When the wind turbines are modeled using the AFM, the four different mesh resolutions employed for the isolated wind turbine simulations described in Sect. 3 are used for each wind farm configuration. This allows the sensitivity of both the AFM and the LES to the grid resolution to be studied. For each wind farm, a baseline case employing the ADM is conducted, characterized by a mesh resolution of $\mathrm{30}\times \mathrm{12.5}\times \mathrm{10}$ m in the streamwise, spanwise, and vertical directions, respectively, similar to previous numerical setups by Stipa et al. (2024c) and Lanzilao and Meyers (2023). In all cases, wind turbines are immersed in the same neutral ABL described in Appendix B. In particular, while further advancing in time both precursor cases characterized by a mesh resolution of $\mathrm{15}\times \mathrm{15}\times \mathrm{10}$ and $\mathrm{50}\times \mathrm{50}\times \mathrm{10}$ m for 20 000 additional seconds, y−z slices of velocity are saved at each time step in what is referred to as the inflow database. The coarser precursor employs the wall model correction described in Appendix A, where the value of u^* is obtained from the finer precursor. The generated inflow databases are used to prescribe the inlet velocity field for the wind farm simulations by linearly interpolating in time from the two closest available time samples, as well as bi-linearly interpolating in space from the precursor to the successor two-dimensional boundary meshes. In order to avoid the mismatch in the shear stress profile when the ratio between the target and source mesh size is greater than 2, the inflow data from the $\mathrm{15}\times \mathrm{15}\times \mathrm{10}$ m precursor are only used for the successor cases characterized by a grid spacing of $\mathrm{30}\times \mathrm{12.5}\times \mathrm{10}$ and $\mathrm{30}\times \mathrm{20}\times \mathrm{10}$ m, whereas the wind farm analyses involving a cell size of $\mathrm{40}\times \mathrm{40}\times \mathrm{10}$ and $\mathrm{60}\times \mathrm{60}\times \mathrm{10}$ m use the inflow data obtained from the $\mathrm{50}\times \mathrm{50}\times \mathrm{10}$ m precursor with wall model correction. At the outlet, all wind farm simulations employ a zero gradient condition, while the remaining boundaries are treated similarly to their respective precursor simulation, which are detailed in Appendix B.

The wind farm is characterized by a streamwise spacing S_x of 630 m (five rotor diameters) and a spanwise spacing S_y of 600 m (4.76 rotor diameters). For the aligned case, this yields a total size $L_x^f\times L_y^f$ of 2.52×2.4 km in the streamwise and spanwise directions, respectively, while the same values of S_x and S_y determine a total wind farm size of 2.52×2.7 km for the staggered case, as rows 2 and 4 are shifted by $S_y/\mathrm{2}$ in the spanwise direction. The successor domain is $\mathrm{15.6}\times \mathrm{8.4}\times \mathrm{0.7}$ km, arranged such that 3 km are left on each side of the wind farm (for the staggered case they reduce to 2.7 km on the right side) and between the domain inlet and the first wind farm row. This leads to 10.08 km between the last wind farm row and the domain outlet, which are used to track the wind farm wake evolution. All four mesh configurations are uniform, leading to the number of cells for each simulation reported in Table 2. The AFM settings are based on the results from Sect. 3: hence, $r_{\mathrm{1}/\mathrm{2}}/R=\mathrm{1}$ and s=6. The streamwise standard deviation σ is set to be consistent with the ADM simulations, where the isotropic standard deviation is set to 1.5Δ_y=18.75 m. This corresponds to $\mathit{\sigma }/\mathrm{\Delta }x=\mathrm{0.625}$ for the finer AFM case, which is maintained for all mesh resolutions. In total, we run eight AFM simulations and two ADM simulations. All analyses are advanced in time for 20 000 s, and flow statistics are averaged over the last 15 000 s.

Table 2 Number of mesh cells and total number of degrees of freedom for each mesh configuration of the wind farm cases.

Figure 7 Contours of instantaneous hub-height velocity field from the ADM and AFM simulations. Top and bottom rows correspond to the staggered and aligned wind farm layouts, respectively.

Figure 8 Contours of time-averaged hub-height velocity field from the ADM and AFM simulations. Top and bottom rows correspond to the staggered and aligned wind farm layouts, respectively.

In Fig. 7, the contours of the instantaneous velocity field at the hub height are reported for all simulations apart from the $\mathrm{30}\times \mathrm{12.5}\times \mathrm{10}$ m AFM cases. Individual wind turbine wakes show a lower tendency to meander for the AFM simulations characterized by a lower horizontal grid resolution (40×40 and 60×60 m). Although this may result in a slower individual wake recovery than in the higher-mesh-resolution cases, it is expected, as reducing the grid resolution increasingly filters both the incoming ABL turbulence and the fine flow features that characterize each individual turbine wake. However, as can be noticed by looking at the mean hub-height velocity field reported in Fig. 8, once turbine wakes have merged, the wake of the entire cluster is less sensitive to the grid's ability to capture the evolution of each individual turbine wake, and cases characterized by different turbine model and mesh resolution are in very good agreement.

From the row-averaged thrust and power reported for all cases in Fig. 9, it can be noticed that the ADM and the AFM yield very similar results for the $\mathrm{30}\times \mathrm{12.5}\times \mathrm{10}$ m resolution for both the aligned and staggered cases. At the waked rows, the AFM model predicts a slightly lower values of thrust and power. This can be attributed to the employed velocity sampling strategy, according to which the wind speed is sampled at a single location. Notably, when an upstream-aligned turbine is present, the sampling location coincides with the wake centerline, leading to a lower sampled velocity. At the non-waked rows, AFM predictions are fairly independent of the grid spacing. Conversely, at the waked rows, the AFM predicts lower thrust and power as the grid resolution is reduced. The reason for such underestimation follows from the lower tendency for individual turbine wakes to meander and hence to mix when the grid size is increased.

Figure 9 Row- and time-averaged thrust (a, c) and power (b, d) distributions obtained for the aligned (a, b) and staggered (c, d) wind farm layouts. The black line refers to the ADM case, while the blue, orange, green, and red lines correspond to the results obtained using the AFM with rotor disk sampling on the 30×12.5, 30×20, 40×40, and 60×60 m horizontal mesh sizes.

Figure 10 shows the time-averaged spanwise velocity profiles at the hub height at different streamwise locations inside the wind farm and in the wake for both the aligned and staggered wind farm layouts. As can be observed, the spanwise velocity profiles predicted using the AFM are in good agreement with the ADM results both inside and downstream of the wind farm for both the staggered and the aligned layouts.

Figure 10 Time-averaged spanwise velocity profiles at the hub height, sampled at different streamwise locations inside the wind farm and in the wake, for the aligned (a) and staggered (b) layouts. Wind farm locations are identified with the row ID, while wake locations are identified by their distance in rotor diameters from the last wind farm row. The black line refers to the ADM case, while the blue, orange, green, and red lines correspond to the results obtained using the AFM with rotor disk sampling on the 30×12.5, 30×20, 40×40, and 60×60 m horizontal mesh sizes.

Figure 11 reports the time-averaged hub-height velocity as a function of the streamwise coordinate, further averaged over the wind farm width. As can be noticed, except from the AFM results obtained with the 60×60 m horizontal mesh resolution, the velocity evolution agrees well with that predicted by using the ADM model, especially upstream of the wind farm and in the wake. In fact, it can be argued that the wind speed in the wind farm wake is fairly independent of the grid spacing after the individual turbine wakes have merged.

Figure 11 Streamwise evolution of time-averaged hub-height velocity, further averaged over $y=\pm \mathrm{2.5}$ km, for the aligned (a) and staggered (b) layouts. The black line refers to the ADM case, while the blue, orange, green, and red lines correspond to the results obtained using the AFM with 30×12.5, 30×20, 40×40, and 60×60 m horizontal mesh sizes. Vertical dashed lines correspond to the first and last wind farm rows.

Figure 12 Streamwise derivative of the time-averaged hub-height velocity, further averaged over the wind farm width, for the aligned (a, b) and staggered (c, d) layouts. The black line refers to the $\mathrm{30}\times \mathrm{12.5}\times \mathrm{10}$ m ADM case, while the blue, orange, green, and red lines correspond to the results obtained using the AFM with rotor disk sampling on the 30×12.5, 30×20, 40×40, and 60×60 m horizontal mesh sizes. Panels (a) and (c) refer to the wind farm region, and (b) and (d) refer to the wake region.

Figure 13 Time-averaged vertical shear stress profiles, further averaged over the wind farm width, for the aligned (a) and staggered (b) layouts. Wind farm locations are identified with the row ID, while wake locations are identified by their distance in rotor diameters from the last wind farm row. The black line refers to the ADM case, while the blue, orange, green, and red lines correspond to the results obtained using the AFM with rotor disk sampling on the 30×12.5, 30×20, 40×40, and 60×60 m horizontal mesh sizes. Horizontal dashed lines refer to h_hub±R.

To further expand on this and to assess the differences in wake recovery predicted by the AFM using different grid sizes, we report in Fig. 12 the streamwise derivative of the mean velocity previously shown in Fig. 11. As can be noticed, while the $\mathrm{40}\times \mathrm{40}\times \mathrm{10}$ and the $\mathrm{60}\times \mathrm{60}\times \mathrm{10}$ m mesh resolutions underpredict wake recovery inside the wind farm with respect to the finer meshes, wake recovery is well captured by all mesh resolutions after ≈1 wind farm length downstream of the last turbine row, with the largest deviations observed for the 60×60 m AFM case with the aligned wind farm layout.

One of the effects of the wind farm on the ABL flow is to increase vertical turbulent mixing by enhancing the level of shear stress. This enhances momentum entrainment from above the wind farm, which plays an important role in the wake recovery of the entire wind farm. Figure 13 shows the time-averaged vertical shear stress profiles, further averaged over the wind farm width, at different streamwise locations inside the wind farm and in the wake. Notably, the ADM and AFM results are generally in good agreement except for the AFM case characterized by the largest grid spacing, which underpredicts the shear stress profile evolution inside the wind farm. In the wind farm wake, all cases are in very good agreement for both the aligned and staggered layouts. The same conclusions can be drawn from the time-averaged vertical velocity profile at y=0 m, reported in Fig. 14. In addition, it is evident from this figure how the first cell velocity strongly depends on the employed grid spacing when the LOTW scaling is not captured according to Brasseur and Wei (2010) criteria. This applies to the $\mathrm{40}\times \mathrm{40}\times \mathrm{10}$ and the $\mathrm{60}\times \mathrm{60}\times \mathrm{10}$ m grids, where the velocity at the first cell decreases as the horizontal grid size increases. However, as discussed in Appendix A and Appendix B and confirmed here, this does not impair the results of the wind farm simulations, especially when the wall shear stress experienced away from the wind farm matches that of a simulation complying with the LOTW scaling criteria. Moreover, even the coarsest grid seems to capture the shear stress perturbation generated by the wind farm, which explains why wind farm wake recovery is also well captured by all values of mesh resolution.

Figure 14 Time-averaged vertical velocity profiles at y=0 m for the aligned (a) and staggered (b) layouts. Wind farm locations are identified with the row ID, while wake locations are identified by their distance in rotor diameters from the last wind farm row. The black line refers to the ADM case, while the blue, orange, green, and red lines correspond to the results obtained using the AFM with rotor disk sampling on the 30×12.5, 30×20, 40×40, and 60×60 m horizontal mesh sizes. Horizontal dashed lines indicate h_hub±R.

In summary, regarding the accuracy of the AFM with respect to the ADM, the two approaches are practically equivalent if the same mesh resolution is employed. For the $\mathrm{40}\times \mathrm{40}\times \mathrm{10}$ and the $\mathrm{60}\times \mathrm{60}\times \mathrm{10}$ m grid sizes, the AFM captures the wind farm power at the non-waked rows, while power is underpredicted at the waked turbines. We argue that this is not an issue of the AFM, but it is rather attributable to the inability to properly capture individual wake meandering by these coarser grids, leading to a slower recovery of individual turbine wakes. Nevertheless, all values of grid resolution can accurately capture the velocity distribution both upstream and downstream of the wind farm. Inside the wind farm, velocity profiles agree reasonably well with those predicted by the ADM except for the AFM case characterized by a resolution of $\mathrm{60}\times \mathrm{60}\times \mathrm{10}$ m, which overpredicts wake deficit. As a consequence, a mesh characterized by a horizontal resolution of 40–60 m may be employed for those problems where an array of interest is waked by an upwind wind farm to discretize the flow region upstream of the array of interest. Conversely, the wake of the target turbine array should be discretized with no more than a 30 m horizontal resolution in order to properly capture individual wake interactions. Notably, a 30 m horizontal resolution corresponds to four cells along the rotor diameter for the wind turbine employed in the present study. This is still too coarse to use the ADM, highlighting the cost-saving potential of the AFM.

5 Farm–farm interaction

In this section, we conduct simulations of two interacting wind farm clusters with the objective of understanding if the AFM can be used to model the upstream wind farm cluster at a low computational cost, when the main focus is on a downstream waked wind farm. We choose an idealized case where two aligned wind farms corresponding to the staggered layout of Sect. 4 are separated by a distance of 5 km. The domain extends for $\mathrm{23}\times \mathrm{8.4}\times \mathrm{0.7}$ km in the streamwise, spanwise, and vertical directions, respectively. In a first case, all turbines are modeled using the ADM and the mesh resolution is set to $\mathrm{30}\times \mathrm{12.5}\times \mathrm{10}$ m. All boundaries are treated similarly to the isolated wind farm simulations described in Sect. 4. The inflow data correspond to the fully neutral ABL described in Appendix B, where the grid spacing is set to $\mathrm{15}\times \mathrm{15}\times \mathrm{10}$ m. A second simulation employs one-way-coupled nested domains using the technique described in Sect. 2.3. The size of the outer domain coincides with that of the ADM case, but it is discretized using a $\mathrm{50}\times \mathrm{50}\times \mathrm{10}$ m mesh resolution instead. Moreover, both wind farms are modeled using the AFM. The inner domain, characterized by a grid resolution of $\mathrm{30}\times \mathrm{12.5}\times \mathrm{10}$ m, extends for $\mathrm{16}\times \mathrm{8}\times \mathrm{0.7}$ km, and its inlet boundary is located 4 km after the start of the first wind farm. Here, wind turbines are modeled using the ADM, and velocity is interpolated at the inlet, top, and side boundaries from the outer domain. At the wall, a wall model based on the classic Monin and Obukhov (1954) similarity theory is employed, while the outlet is treated similarly to the outer domain.

Notably, since we employ a one-way nested domain, the downwind wind farm is modeled in both the outer and inner domains to capture its effect on the upstream cluster. In order to highlight the effects of mapping the inflow database between different precursor and successor grid sizes, further detailed in Appendix A, the simulation employing the AFM is carried out twice, using both the same inflow database as the ADM simulation and the inflow database obtained employing the precursor mesh characterized by a resolution of $\mathrm{50}\times \mathrm{50}\times \mathrm{10}$ m with wall model correction. Both precursor simulations are described in Appendix B and correspond to those used in Sect. 4. The boundary conditions in the outer domain are the same as the ADM case except for the wall. Specifically, when the inflow database generated from the coarse precursor is used, the wall shear stress is applied by fixing the friction velocity as described in Appendix A, calculated from the finer precursor. Conversely, when the inflow database generated from the finer precursor is used to prescribe the inflow to the outer domain, friction velocity is calculated locally by means of classic Monin and Obukhov (1954) similarity theory, also described in Appendix A. Table 3 summarizes the main features of the three simulations, such as the employed turbine model, the total number of mesh cells, and the inflow data used to prescribe the ABL flow at the inlet. Throughout the remainder of this section, simulations 1, 2, and 3 in Table 3 will be referred to as ADM, AFM, and AFM with coarse inflow. All simulations are carried out for 20 000 s, and flow statistics are gathered for the last 15 000 s.

Table 3 Summary of turbine model, mesh size used in the precursor and successor simulations, and number of degrees of freedom for the wind farm simulations conducted in the present section. Farm A and farm B are the upwind and downwind wind farms, respectively.

n/a – not applicable.

Figure 15 Instantaneous (a, c) and time-averaged (b, d) hub-height wind speed obtained from cases 1 (a, b) and 3 (c, d) of Table 3. The dashed black line shows the horizontal size of the inner domain in the AFM simulation.

Figure 15 shows the contours of instantaneous and time-averaged hub-height velocity for the AFM with coarse inflow and ADM cases. Although the coarser mesh used in the outer domain of the AFM case filters out the fine turbulence structures that are instead resolved in the ADM case, a very good agreement exists between the two on a qualitative level. A more quantitative comparison between all cases is reported in Fig. 16 by showing the mean vertical velocity profile at y=0 and the mean shear stress profiles, further averaged over $y=\pm \mathrm{2.5}$ km.

Figure 16 (a) Time-averaged vertical profiles of velocity magnitude at y=0 and at different streamwise locations. (b) Time-averaged vertical shear stress profile, further averaged between $y=\pm \mathrm{2.5}$ km, at different streamwise locations. In both panels, streamwise locations are identified by the row number and wind farm ID (A for upstream and B for downstream) when inside a wind farm or by their distance in rotor diameters from the last row of the closest upstream wind farm. The dashed line in (a) corresponds to the freestream velocity obtained from the precursor simulations characterized by the $\mathrm{15}\times \mathrm{15}\times \mathrm{10}$ m grid resolution. Horizontal dashed lines refer to h_hub±R.

First, the vertical velocity profiles predicted by the ADM and AFM simulation that employs the coarse inflow agree well at every streamwise location. Only a slightly higher velocity is observed up to 20 rotor diameters downstream of the first wind farm. Conversely, the AFM simulation that uses the finer inflow consistently underestimates the wake deficit by the same amount in both wind farm wakes. Looking at the shear stress profile, the largest difference can be observed around the first wind farm, while all profiles collapse after 30 rotor diameters downstream of the first cluster. In particular, it can be noticed how the shear stress profile is strongly underpredicted at the first wind farm row of the upwind wind farm by the AFM case where the finer inflow data are used. This issue, analyzed in Appendix A, reaches all the way to the inlet of the outer domain and it is attributed to the velocity mapping from the inflow database. Specifically, when this is interpolated from the finer precursor mesh to the coarser successor grid, non-solenoidal velocity fluctuations are produced which are subsequently altered by the pressure iteration of our solver when it corrects the velocity field in order to make it divergence-free. Conversely, when the coarser inflow data and the wall model correction are used, the inlet shear stress profile agrees with that resulting from the ADM case. Moreover, both AFM cases underestimate – to a lesser extent when the coarser inflow is used – the perturbation in shear stress inside the first wind farm. Notably, this is an expected mechanism when the mesh size is increased, since also the LES filter size grows and more eddies are modeled by an increase in the eddy viscosity. Further downstream, as the flow enters the inner domain, the shear stress profiles from the different cases are in good agreement with each other.

Figure 17 Streamwise evolution of time-averaged hub-height wind speed, further averaged between $y=\pm \mathrm{2.5}$ km from the three cases conducted in the present sections. Vertical dashed lines indicate the two wind farms. A magnification of (a) in the induction region of the first wind farm is depicted in (b).

The effect on velocity of the reduction in shear stress at the domain inlet operated by the mapping procedure when the source and target grids are very different from each other can be clearly visualized in Fig. 17. In fact, looking at the blockage region of the AFM case employing the inflow database generated with the finer grid, one can see that the wind speed increases after the inlet before reaching the wind farm. This phenomenon, explained in Appendix A, is due to the fact that a reduction in shear stress causes a deformation of the velocity profile in order to satisfy the momentum budget inside the boundary layer. This may potentially induce also a spanwise velocity component, as the momentum source terms representing the constant driving pressure gradient in the successor simulation are averaged from the precursor. For this reason, they require the same shear stress profile to fulfil the horizontal momentum balance. The most important consequence of such an issue is that the blockage region is completely misrepresented and the freestream velocity experienced by the upstream wind farm is increased. However, by applying the inflow data calculated using the coarse mesh and the correction to the wall model, the shear stress profile is not altered and the velocity in the induction region is correctly captured. In general, using this approach results in a much better agreement with results from the ADM case. The largest differences are observed in the wake of the first wind farm and are likely due to the wind farm thrust underprediction by the AFM when the horizontal resolution exceeds ≈30 m. To some extent, also the fact that the AFM in general predicts a slightly lower freestream wind speed – and thus a lower thrust – due to the employed sampling method may play a minor role.

Figure 18 Row-averaged thrust (a, b) and power (c, d) distributions from the three cases described in this section for the upstream (a, c) and downwind wind farm (b, d).

Finally, Fig. 18 shows the row-averaged wind turbine thrust and power for the two wind farms obtained from the three different cases. For the two AFM cases, data correspond to the AFM and ADM models for the upwind and downwind wind farms, respectively. In general, the AFM results point to the same conclusions as the isolated wind farm cases presented in Sect. 4, i.e. that a horizontal resolution of 50×50 m is not sufficient to capture the absolute thrust and power. However, trends are reasonably captured, and the error reduces when the coarser inflow data are used, as the shear stress profile and thus the turbulence intensity level are in better agreement with the ADM case. In addition, thrust and power from the downwind wind farm are only slightly overestimated, again with the error decreasing if the coarse inflow data and wall model correction are used. These results suggest that the AFM is a good candidate model for problems involving one or more wind farms waking a downstream cluster of interest. In particular, the less stringent requirement on mesh resolution imposed by the AFM reduces the overall computational cost while still capturing the cluster wake evolution with reasonable accuracy.

6 Wind-farm-induced atmospheric gravity waves

In this section, we conduct LESs of wind-farm-induced atmospheric gravity waves (AGWs), investigating the ability of the AFM to capture the AGW evolution in the free atmosphere. In particular, we simulate a wind farm immersed in a conventionally neutral boundary layer (CNBL), i.e. a boundary layer developing against a potential temperature stratification characterized by a neutral region, followed by a capping inversion layer with strength Δθ and thickness Δh, centered at H, and a linear stable lapse rate γ aloft. For the CNBL, we chose Δθ=5 K, Δh=100 m, γ=4 K km⁻¹, and H=500 m, which also coincides with the height of the boundary layer, as its growth is limited by the capping inversion. These values are frequently observed offshore in the North Sea, as pointed out by Lanzilao and Meyers (2023). The equivalent roughness height z₀ is set to 0.0001 m, and the reference potential temperature θ_ref is equal to 300 K. The precursor simulation uses a velocity controller which aims at maintaining a reference velocity of 9 m s⁻¹ at the wind turbine hub height of 90 m, as turbines correspond to the NREL 5MW reference turbine (Jonkman et al., 2009). The precursor simulation is advanced for 100 000 s in order to spin up turbulence, and the horizontally averaged potential temperature profile is kept constant by the temperature controller described in Stipa et al. (2024b). The flow is initialized with a uniform velocity of 9 m s⁻¹, while temperature follows the model developed by Rampanelli and Zardi (2004). Geostrophic damping using the same settings employed by Stipa et al. (2024b) is applied to remove inertial oscillations that may arise when the initial geostrophic speed is not in geostrophic balance, a condition that cannot be avoided when forcing the wind speed at a height located inside the boundary layer. The precursor domain extends for $\mathrm{6}\times \mathrm{6}\times \mathrm{1}$ km in the streamwise, spanwise, and vertical directions, respectively, and it is discretized using a grid size of 15×15 m in the horizontal directions. Below the start of the inversion layer, the vertical grid size is set to 10 m. From 450 to 500 m the grid is reduced to 5 m and then increased again to 10 m at 550 m to capture the Ellison scale within the inversion layer (Allaerts and Meyers, 2017). The 10 m resolution is then maintained until the upper boundary. After the first 100 000 s of simulation, statistics are averaged for 20 000 s and y−z flow sections are saved at each time step to form the inflow database. The profiles of wind speed and direction, shear stress, and potential temperature from the precursor phase are reported in Fig. 19, while quantitative data are summarized in Table 4.

Figure 19 Vertical profiles of wind speed magnitude (a), wind veer (b), potential temperature (c), and shear stress (d) for the CNBL precursor. Data are averaged from 100 000 to 120 000 s.

Table 4 Values of the geostrophic wind, friction velocity, minimum heat flux q_min within the boundary layer, and geostrophic wind angle obtained from the CNBL precursor.

Regarding the wind farm simulations, their setup is sketched in Fig. 20. For the ADM case, this consists of a domain that extends for 22.62×12 km in the streamwise and spanwise directions, respectively. The wind farm corresponds to the staggered layout described in Sect. 4, where the first row is located at x=0 and separated by 10 km from the inlet boundary. Moreover, the wind turbines at the first-row sides are located at a distance of 4.8 km from the lateral boundaries. Regarding the vertical domain size, this is dictated by the simulation's ability to resolve AGWs. In fact, the total domain height should be at least twice as high as the expected gravity wavelength λ_z (this parameter can be estimated as ${\mathit{\lambda }}_z=\mathrm{2}\mathit{\pi }G/N$, where N is the Brunt–Väisälä frequency and G is the geostrophic wind). Moreover, a Rayleigh damping layer should be used at the upper boundary to avoid AGW reflection, characterized by a layer depth ≥λ_z (Lanzilao and Meyers, 2022a). With reference to the CNBL parameters used in the present study, λ_z≈5.2 km. Hence, the domain height has been set to 14 km, while the start of the Rayleigh damping region has been placed at 7 km (blue box in Fig. 20). The mesh resolution in the vertical direction follows the precursor simulations below 1 km, while it is stretched up to 200 m in the Rayleigh damping region. Specific details are provided in Table 5. In order to also avoid AGW reflections from the inlet boundary, a fringe region and an advection damping region characterized by the same activation functions adopted by Lanzilao and Meyers (2022a) are used (magenta and orange boxes in Fig. 20, respectively), and their parameters are reported in Table 6. Following the same authors, the Rayleigh and fringe region damping coefficients have been set to ν_RDL=0.035 s⁻¹ and ν_FR=0.03 s⁻¹, respectively. The successor simulations employ lateral periodic boundary conditions, a slip wall at the upper boundary, and a wall model based on the classic Monin and Obukhov (1954) similarity theory at the wall.

Table 5 Vertical discretization for successor and concurrent precursor simulations. The parameter N indicates the number of cells in each mesh layer. These extend from z_s to z_e and are characterized by a cell size Δz.

Table 6 Fringe and advection damping region information.

Table 7 Summary of turbine model, mesh size used in the precursor and successor simulations, and number of degrees of freedom for the wind farm simulations conducted in the present section.

As spatially and temporally resolved velocity and temperature fields that are unperturbed by the wind farm are required in the fringe region to compute the damping source terms, a concurrent precursor simulation characterized by a domain size coincident with the fringe region ($\mathrm{3}\times \mathrm{6}\times \mathrm{14}$ km) and having the same mesh resolution as the wind farm domain is carried on in sync with the latter. The concurrent precursor uses inflow slices saved from the CNBL simulation described above to enforce a time-resolved inflow condition. At the outlet we apply a zero gradient boundary condition, while all remaining boundaries are treated similarly to the wind farm simulation. Since the flow slices available from the pre-computed inflow database are 6×1 km large in the spanwise and vertical directions, their data are tiled two times along y and extrapolated along z in order to be mapped at the concurrent precursor inlet.

Figure 20 Location of the Rayleigh damping layer (blue), fringe region (magenta), and advection damping region (orange) relative to the domain. The ADM case (a) only uses a single domain, while the AFM case (b) employs the same outer domain, characterized by a coarser grid resolution and an inner domain (black) featuring the same grid spacing as the ADM case. Wind turbines are modeled using the AFM in the outer domain, while the ADM is employed in the inner domain.

Figure 21 Time-averaged hub-height velocity (a, c) and pressure perturbation (b, d), calculated as $p_p(x,y,z)=p(x,y,z)-p(x_{\mathrm{\infty }},y,z)$, where $x_{\mathrm{\infty }}=-\mathrm{7}$ km and $p(x,y,z)$ is the time-averaged pressure. Panels (a, b) and (c, d) correspond to the ADM and AFM cases, respectively. In the latter, the horizontal inner domain size is identified by the dashed black line. The end of the fringe and advection damping regions are identified by the dashed blue and orange lines, respectively.

Figure 22 Time-averaged vertical velocity (a, c) and pressure perturbation (b, d) at y=0. Panels (a, b) and (c, d) correspond to the ADM and AFM cases, respectively. In the latter, the horizontal inner domain size is identified by the dashed black line. The end of the fringe and advection damping regions are identified by the dashed blue and orange lines, respectively, while the start of the Rayleigh damping region corresponds to the dashed green line.

Regarding the two AFM simulations, they employ two one-way nested domains. The outer domain setup coincides with the ADM simulation described above, with the only differences being that the AFM is used in place of the ADM and that the horizontal mesh resolution is coarsened to 50×50 m. The inner domain (black box in the left panel of Fig. 20) extends for $\mathrm{15.6}\times \mathrm{8.4}\times \mathrm{0.7}$ km, with the inlet boundary placed at $x=-\mathrm{3}$ km, and it is discretized using a mesh resolution of $\mathrm{30}\times \mathrm{12.5}\times \mathrm{10}$ m. Velocity and potential temperature are interpolated at the lateral, upper, and inlet boundaries from the outer domain, while the outlet and bottom boundaries use a zero gradient and a wall model based on classic Monin and Obukhov (1954) similarity theory, respectively. In the inner domain, wind turbines are modeled using the ADM. This arrangement allows AGWs to be captured in the outer domain, where they are forced by the AFM, and their effects to be transferred to the inner domain by interpolating velocity and potential temperature from the outer grid. Notably, in the outer domain, the same inflow data used for the ADM simulation are used to provide an inlet boundary condition to the concurrent precursor. This means that data are mapped from a y−z grid resolution of 15×10 m, employed for the CNBL simulation, to a mesh spacing of 50×10 m. The issue described in Appendix A related to the modification of the shear stress profile by the mapping interpolation is mitigated by imposing the value of u^* reported in Table 4 for both the concurrent precursor and outer domains to correctly match the freestream shear stress profile. A summary is given in Table 7 regarding the turbine model, precursor and successor grid size (that of the concurrent precursor coincides with the successor), and number of degrees of freedom employed for both the ADM and AFM simulations. Both the ADM and AFM cases are advanced in time for 20 000 s, and statistics are gathered over the last 15 000 s of simulation.

Figure 21 shows contours of time-averaged velocity at the hub height and the pressure perturbations produced at the same height by the internal and interface waves triggered by the wind farm in the free atmosphere and within the inversion layer, respectively. As can be noticed, the developed setup combining grid nesting with the AFM agrees well with the ADM results obtained with a more conventional design of the numerical simulation. Differences are only observed inside the fringe region, i.e. in a non-physical portion of the domain, where pressure perturbations are reduced for the AFM case. The reasons for such difference are presently unknown to the authors, but it seems that the fringe region performs better when employing a coarser grid resolution. Regarding the AGW patterns produced in the vertical velocity field and in the pressure perturbation field, they can be visualized in Fig. 22 on a x−z plane located at y=0. Also in this case, they agree extremely well between the two simulations. For instance, the perturbations in pressure resulting from the AFM case seem to be slightly lower than those predicted using the ADM, but, as will be shown in the following analysis, they do not lead to a visible alteration of the results, both in terms of wind speed and turbine quantities, when these are compared against the data extracted from the ADM case.

In this regard, Fig. 23 reports the vertical flow perturbation, magnified 10 times, at different heights. The blue lines correspond to the inversion layer displacement; dotted and dashed lines have been obtained using data from the ADM case and from the outer domain in the AFM simulation, respectively. As can be noticed, the perturbations obtained in the free atmosphere using the AFM are almost identical to those obtained when wind turbines are modeled with the ADM, suggesting that AGWs are not sensitive to how accurately the simulation captures the turbulent flow inside the boundary layer. These AGW-induced vertical flow perturbations are transferred to the inner domain when interpolating the velocity and temperature fields at its boundaries in the AFM simulation, allowing AGW effects on both the wind farm and the boundary layer flow to be modeled.

Figure 23 Vertical flow perturbation at different heights for the ADM (dotted line) and AFM (dashed line) magnified 10 times.

Figure 24 plots the vertical velocity profiles at y=0 and the vertical shear stress profiles, averaged over $y=\pm \mathrm{2.5}$ km, at different streamwise location inside and downstream of the wind farm. Data from the AFM simulation are entirely contained within the inner domain, where wind turbines are modeled using the ADM. For instance, the coarse resolution in the outer domain in the AFM simulation does not alter the incoming boundary layer flow when data are compared against the ADM case. Only small differences in shear stress and velocity exist at the first wind farm row, but these are soon removed by the higher grid resolution adopted in the inner domain and are not propagated downstream. In general, the two simulations are in very good agreement despite the AFM case involving only 30 % of the number of degrees of freedom used in the ADM case. The effect of AGWs on the boundary layer flow can be appreciated in Fig. 25, where the time-averaged and hub-height velocity and pressure fields, further averaged over the wind farm width, are displayed. First, the pressure oscillations in the wind farm wake induced by the lee waves shown in Fig. 21 produce oscillations in velocity, leading to an intermittent recovery of the wind farm wake. Moreover, if compared to previous results obtained with the same wind farm in Sect. 4 and a boundary layer height of 0.7 km, wind farm blockage is greatly increased for the atmospheric conditions analyzed in this section, and important reductions in velocity can be observed up to ≈0.5 km upstream of the first wind farm row. Note that, according to Smith (2024), a fully neutral boundary layer where the domain height coincides with the ABL height corresponds to the rigid-lid approximation of very high stratification above the boundary layer. In particular, as shown in Stipa et al. (2024a), allowing the inversion layer to displace according to the AGW solution in the free atmosphere modifies the pressure field around the wind farm, yielding different values of blockage and individual wake recovery inside the wind farm. Finally, Fig. 26 reports the row-averaged turbine power and thrust from the two simulations, showing an excellent agreement between the two.

Figure 24 (a) Time-averaged vertical profiles of velocity magnitude at y=0 and at different streamwise locations. (b) Time-averaged vertical shear stress profile, further averaged between $y=\pm \mathrm{2.5}$ km, at different streamwise locations. In both panels, streamwise locations are identified by the row number when inside a wind farm or by their distance in rotor diameters from the last wind farm row. The dashed line in (a) corresponds to the freestream velocity obtained from the precursor simulation. Horizontal dashed lines refer to h_hub±R.

Figure 25 Streamwise evolution of time-averaged hub-height velocity (a) and hub-height pressure perturbation (b), further averaged between $y=\pm \mathrm{2.5}$ km. Vertical dashed lines refer to the first and lasts wind farm rows.

Figure 26 Row- and time-averaged turbine thrust (a) and power (b) for the two simulations presented in this section.

Overall, these results demonstrate that, if our focus is only on the atmospheric flow solution, the AFM alone (i.e. without the need of an inner domain) is sufficient to accurately capture AGWs in the free atmosphere and their effects on the ABL flow at a reduced computational cost. However, if accurate wind turbine information is required, an inner domain characterized by a higher grid resolution can be placed around the wind farm, and the outer flow – which contains the AGW solution – can be interpolated at the boundaries to model AGW effects on the wind farm performance. Notably, this is similar to the model proposed by Stipa et al. (2024a) to account for AGW effects on the wind farm flow without extending the domain into the free atmosphere. However, while in Stipa et al. (2024a) the top boundary coincides with the inversion layer and is physically displaced in order to enforce a slip boundary condition, here the slip boundary condition – which corresponds to no penetration – is replaced by interpolating the velocity from the outer domain, thus allowing some degree of permeability at the top boundary.

7 Conclusions

In this study, we introduced the actuator farm model (AFM), a new parametrization that allows the aerodynamics of wind turbines in the context of wind farm LESs to be captured. Unlike similar models such as the actuator disk model (ADM) or the actuator line model (ALM), in the AFM wind turbines are represented using a single actuator point, located at the rotor center, and only two to three mesh cells are required along the rotor diameter. The turbine force is distributed to the surrounding cells by means of a new projection function whose spatial support is axisymmetric in the rotor plane and characterized by a Gaussian decay in the streamwise direction. The size of the spatial support is controlled by means of three free parameters, namely the half-decay radius on the rotor plane $r_{\mathrm{1}/\mathrm{2}}$, the smoothness s, and the streamwise standard deviation σ.

To find the best set of parameters that allow us to obtain similar wake deficit profiles and turbine thrust and power to those predicted using the ADM, we conducted simulations of an isolated NREL 5MW wind turbine in uniform inflow using different values of the horizontal grid spacing. In particular, while σ is chosen using existing best practices from the ADM, our results show that $r_{\mathrm{1}/\mathrm{2}}$ should be approximately the size of the turbine radius R, while values of s should lie between 6 and 10. With this choice of AFM projection parameters, results are fairly independent of the horizontal grid spacing up to a resolution of 60×60 m, i.e. $\mathrm{\Delta }x/R\approx \mathrm{2}$.

The optimal set of parameters ($r_{\mathrm{1}/\mathrm{2}}=R$ and s=6) was used to investigate the AFM performance in predicting the flow around a wind farm with 25 NREL 5MW turbines organized in five rows and five columns in both an aligned and a staggered layout, using horizontal grid spacings of 30×12.5, 30×20, 40×40, and 60×60 m. As the wind farm is immersed in a fully neutral ABL, the time-resolved inflow condition is mapped from a previously conducted precursor simulation, where no turbines are present, using linear interpolation both in space and time.

To avoid alteration in the velocity fluctuations when the target grid is more than twice as coarse as the source grid, the precursor used to prescribe the inlet flow to the wind farm simulations with horizontal grid spacings of 40×40 and 60×60 m has been conducted on a grid characterized by a $\mathrm{50}\times \mathrm{50}\times \mathrm{10}$ m resolution. Notably, this does not satisfy the Brasseur and Wei (2010) criteria and also leads to a reduction in the predicted shear stress magnitude. To correct this issue, we proposed a modified wall model that allows us to recover the shear stress profile obtained when the precursor simulation complies with the Brasseur and Wei (2010) criteria. This is achieved by prescribing the friction velocity used to compute the wall shear stress instead of calculating it based on the velocity at the first cell.

Results obtained using the AFM have been compared against ADM predictions made on the finer grid, which satisfies the ADM requirement of having at least 10 grid cells along the rotor diameters. Specifically, when the same or the 30×20 m horizontal grid spacing is employed, the AFM and the ADM essentially predict identical velocity and shear stress profiles around the wind farm. Moreover, row-averaged turbine thrust and power are in excellent agreement. For the 40×40 and the 60×60 m grid spacing, the AFM captures the wind farm power at the non-waked rows, while power is underpredicted at the waked turbines. Nevertheless, all values of grid resolution allow the mean velocity distribution to be captured both upstream and downstream of the wind farm with good accuracy. Therefore, for those problems where a turbine array of interest is waked by an upwind wind farm, the upwind farm may be modeled using the AFM together with a mesh resolution that places two to three cells across the rotor diameter, while the array of interest should be discretized with no less than four cells across the rotor diameter in order to properly capture individual wake interactions and turbine power. Notably, this resolution is still too poor to use the ADM.

Lastly, the AFM is combined with the nested domain technique and used in two wind farm LES applications to demonstrate its ability to drastically reduce the computational cost whilst predicting similar results in terms of flow field and turbine variables. In particular, we conduct simulations of two aligned wind farms immersed in a truly neutral ABL and of a single wind farm that interacts with a conventionally neutral boundary layer by triggering both internal gravity waves in the free atmosphere and surface waves in the capping inversion layer. Conventionally, these applications are rendered computationally intense by the large domain size required to capture the flow physics and by the fine grid resolution imposed by the ADM model. The proposed AFM allows the grid spacing to be increased, leading to a reduced cell count. In particular, both analyses employ a one-way coupling between an outer domain characterized by a horizontal resolution of 50×50 m and a nested inner domain with a 30×12.5 m grid. Notably, only the solution in the inner domain is influenced by the outer domain. Hence, while the outer domain should contain all the relevant physics, the inner domain only provides a refined solution for the region of interest. As a consequence, turbines are modeled using the AFM and ADM in the outer and inner domain, respectively. In both applications, the combined use of the AFM and grid nesting yields velocity, shear stress, and turbine quantities that are in excellent agreement with those obtained using a finer grid and the ADM throughout. Finally, we also highlight that flow perturbations induced in the free atmosphere and within the boundary layer by wind-farm-induced atmospheric gravity waves obtained using the AFM and a coarser grid size agree almost exactly with ADM simulations conducted on a finer grid.

Future studies will involve using the AFM to study the wind farm response to more realistic atmospheric inflow conditions, introduced within the LES using profile assimilation techniques, as well as the mutual interaction of real-world wind farms with neighbouring clusters (offshore) and with complex terrain features (onshore).

Appendix A: Effect of spatial resolution

Inlet boundary conditions in wind farm LESs are often calculated by means of a different simulation, referred to as the precursor. The precursor does not contain any wind turbine and generally employs periodic boundary conditions in the lateral directions. This allows the flow to be recycled for several flow-turnover times until a fully developed ABL characterized by stationary turbulence statistics is reached. For non-idealized or time-varying atmospheric states corresponding to a specific realization of the planetary boundary layer, the precursor simulation can be forced using profile assimilation techniques or two-dimensional boundary data derived from weather models (see Haupt et al., 2023, for a review). The precursor simulation is then used to derive boundary conditions that characterize the incoming flow for the wind farm simulation, referred to as the successor. In its most simple form, two-dimensional sections at a given streamwise coordinate are saved during the precursor at each time step and then used as the inlet boundary condition for the successor simulation. If free-atmosphere stratification is present and atmospheric gravity waves should be resolved, the precursor can be synchronized with the successor and used to prescribe the inlet flow through momentum and temperature source terms applied throughout a fringe region (see Allaerts and Meyers, 2017; Stipa et al., 2024b; Lanzilao and Meyers, 2022a, among others). More complex ways of forcing a successor simulation involve the use of one- or two-way-coupled nested domains such as in WRF LESs (Sanchez Gomez et al., 2022, 2023). In all these cases, the inflow data used to set the inlet boundary condition for the successor should exhibit the correct law-of-the-wall (LOTW) scaling in the velocity profile, as well as the expected shear stress profile and friction velocity for the specific conditions under investigation. In addition, the inflow data should be mapped to the wind farm simulation without being altered by the mapping procedure. In the reminder of this section, the effects of grid coarsening (in both the precursor and successor) on these aspects will be addressed in detail.

Regarding the compliance with LOTW scaling and the recovery of the representative shear stress profile and friction velocity, Brasseur and Wei (2010) identified three criteria that should be satisfied when running ABL simulations:

$$\begin{array}{}\text{(A1)}& {\displaystyle }& {\displaystyle }\mathcal{R}/{\mathcal{R}}^{*}>\mathrm{1},\text{(A2)}& {\displaystyle }& {\displaystyle }{\mathit{Re}}_{\mathrm{LES}}/{\mathit{Re}}_{\mathrm{LES}}^{*}>\mathrm{1},\text{(A3)}& {\displaystyle }& {\displaystyle }{N}_{\mathit{\delta }}/N_{\mathit{\delta }}^{*}>\mathrm{1},\end{array}$$

where N_δ is the number of cells used to resolve the boundary layer, Re_LES is a Reynolds number calculated with a spurious length scale δ_LES arising due to spurious frictional forces in the sub-grid-scale (SGS) model, and ℛ is the ratio between the fluctuating resolved stresses and the modeled SGS stresses at the first cell center. The critical values, identified with ^*, roughly correspond to ${\mathcal{R}}^{*}=\mathrm{1}$, ${\mathit{Re}}_{\mathrm{LES}}^{*}\approx \mathrm{300}$, and $N_{\mathit{\delta }}^{*}\approx \mathrm{50}$ for the neutral ABL. While the above criteria depend on the specific SGS closure employed within the LES, for an eddy-viscosity closure that employs the Smagorinsky model,

$$\begin{array}{}\text{(A4)}& {\displaystyle }& {\displaystyle }\mathcal{R}\approx {\displaystyle \frac{(N_{\mathit{\delta }}-\mathrm{1}){\mathit{\kappa }}^{\mathrm{2}}}{\mathrm{1.05}{N}_{\mathit{\delta }}{C}_{\mathrm{s}}^{\mathrm{2}}{A}_{\mathrm{R}}^{\mathrm{4}/\mathrm{3}}}},\text{(A5)}& {\displaystyle }& {\displaystyle }{\mathit{Re}}_{\mathrm{LES}}\approx {\displaystyle \frac{N_{\mathit{\delta }}}{\mathit{\kappa }}}\left(\mathcal{R}+\mathrm{1}\right),\end{array}$$

where C_s is the Smagorinsky constant, κ is the von Kármán constant, and $A_{\mathrm{R}}=\mathrm{\Delta }h/\mathrm{\Delta }z$ is the cell aspect ratio at the wall, with $\mathrm{\Delta }h=\text{max}(\mathrm{\Delta }x,\mathrm{\Delta }y)$. The criteria expressed by Eqs. (A1) to (A3), calculated assuming an eddy-viscosity closure using the Smagorinsky SGS model with C_s=0.1 and a boundary layer height of H=750 m, are summarized in Table A1 for different values of the grid spacing. Notably, the quantity $\mathcal{R}/{\mathcal{R}}^{*}$ strongly depends on the employed model coefficient and the cell aspect ratio at the wall. As a rule of thumb, to improve LOTW scaling, one should increase the horizontal resolution while keeping the vertical resolution constant or reduce the model coefficient. The criteria expressed by Eqs. (A1) to (A3) can be used to estimate if the LES setup is suitable to capture LOTW scaling, and adhering to these criteria is advisable, particularly when considering the ABL flow. However, numerous studies exist, mainly focused on the wind farm flow, wherein strict adherence to the Brasseur and Wei (2010) criteria is not observed (as seen in Stieren and Stevens, 2022; Stipa et al., 2024c; Lanzilao and Meyers, 2023) but which still offer valuable insight regarding the underlying flow physics. Moreover, it should be noted that while failing to respect Brasseur and Wei (2010) criteria in the precursor simulation may lead to a consistent LOTW mismatch, as the flow is recycled multiple times over the domain, doing so in the successor wind farm simulation only leads to a potential mismatch that develops from the inlet, where the provided inflow is prescribed, to the outlet, where the wind has evolved according to the specific simulation setup. In general, the impact of satisfying Eqs. (A1) to (A3) only in the precursor and not in the wind farm simulation depends on additional factors, such as the size of the computational domain over which the mismatch accumulates as well as the employed SGS closure and numerical schemes. Moreover, it may also depend on the adopted simulation code.

Table A1 Brasseur and Wei (2010) criteria for different precursor mesh sizes, calculated assuming Smagorinsky model coefficient C_s=0.1 and ABL height H=750 m (these are representative values for wind farm LESs). Critical values are set to ${\mathcal{R}}^{*}=\mathrm{1}$, ${\mathit{Re}}_{\mathrm{LES}}^{*}=\mathrm{300}$, and $N_{\mathit{\delta }}^{*}=\mathrm{50}$.

Regarding the mapping from the precursor inflow data to the successor domain inlet, this requires interpolation in both space and time, as the precursor and successor may not have the same mesh at the inlet or may not have advanced with an identical time step size. When this is the case, the two-dimensional inflow data mapped at the successor inlet lose the property of being divergence-free. Consequently, when dealing with an incompressible code, non-solenoidal fluctuations in velocity arising from the mapping will be advected into the internal cells and rendered solenoidal by the Poisson iteration, ultimately modifying the incoming profile of the resolved Reynolds stresses. The result is an imbalance in the momentum equation, wherein the driving pressure gradient is no longer balanced by the resolved shear stress, causing the mean flow in the successor to accelerate or decelerate from the mapped inlet plane depending on whether the error in the resolved stress is positive or negative. We observe this behaviour to be more prominent when interpolating from a finer to coarser mesh – i.e. when mapping yields a loss of information – and when the difference in mesh size is more than a factor of 2. To avoid this problem, continuity-preserving B-spline interpolation proposed by Schroeder et al. (2022) can be used for mapping inlet flow data instead of the classic bi-linear interpolation. However, we have found that this only yields marginal improvement, as spatial interpolation along x is required to render the interpolated flow truly divergence-free and B-spline interpolation along x requires saving three flow sections per time step, tripling the I/O overhead associated with precursor–successor mapping.

Based on the above discussion, it might be concluded that precursor and successor solutions should be conducted with identical spatial and temporal resolution to avoid altering the successor result by the inlet mapping procedure. However, a resolution that sufficiently captures the ABL turbulence in the precursor might be very restrictive for a successor simulation of a very large wind farm in terms of computational cost. This is the case of Maas (2023) and Cheung et al. (2023), who both used a grid size of $\mathrm{20}\times \mathrm{20}\times \mathrm{20}$ m for their precursors and successor analyses. Such a resolution is expected to capture the LOTW according to Eqs. (A1) to (A3) but led to 6.8 and 21.14×10⁹ mesh elements for the above-mentioned studies, respectively. Cheung et al. (2023) additionally used mesh refinement around the wind turbines, making these wind farm simulations perhaps the largest conducted to date. The Brasseur and Wei (2010) criteria were also satisfied for the LES performed by Wu and Porté-Agel (2017), who used a grid spacing of $\mathrm{40}\times \mathrm{16}\times \mathrm{10}$ m for their precursor and successor simulations. An alternative approach is to conduct the precursor using a coarser mesh, coincident with the lowest resolution required by the ADM to resolve the wind turbines. This is the case of Stieren and Stevens (2022) and Lanzilao and Meyers (2023), who used $\mathrm{30}\times \mathrm{15}\times \mathrm{5}$ and $\mathrm{31.25}\times \mathrm{21.74}\times \mathrm{5}$ m, respectively. Stipa et al. (2024c) conducted wind farm simulations using a $\mathrm{30}\times \mathrm{15}\times \mathrm{10}$ m grid spacing, but precursors were carried out on a domain characterized by a $\mathrm{15}\times \mathrm{15}\times \mathrm{10}$ m cells. While these values may not strictly adhere to the criteria defined by Brasseur and Wei (2010), they are reasonable approximations. Any potential LOTW mismatch is expected to be minimal and unlikely to significantly alter the simulation results with respect to a fully LOTW-compliant LES.

If the AFM is employed, the successor mesh can be further coarsened up to a grid spacing of the order of 40–60 m in the horizontal directions. In this case, interpolating the inflow data from a precursor that satisfies the Brasseur and Wei (2010) criteria may lead to the alteration of the shear stress profile described above. On the other hand, when a similar resolution to the successor mesh is employed in the precursor domain, the LOTW mismatch becomes large, especially at the first cell center, potentially affecting the wall shear stress if a wall model is used. In fact, widely adopted wall models based on the classic Monin and Obukhov (1954) similarity theory compute the shear stress at the wall by applying the LOTW locally at the first cell center. The velocity at the first cell is used to compute the friction velocity u^* as

$$\begin{array}{}\text{(A6)}& u^{*}={\displaystyle \frac{\mathit{\kappa }\sqrt{u_{\mathrm{1}}^{\mathrm{2}}+v_{\mathrm{1}}^{\mathrm{2}}}}{\ln \left(\frac{z_{\mathrm{1}}}{z_{\mathrm{0}}}\right)}}\end{array}$$

for the neutral ABL, where subscript 1 indicates quantities evaluated at the first cell center and z₀ is the equivalent roughness height. The wall shear stress is then calculated as

$$\begin{array}{}\text{(A7)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{xz}^{\text{w}}=-u^{*\mathrm{2}}{\displaystyle \frac{u_{\mathrm{1}}}{\sqrt{u_{\mathrm{1}}^{\mathrm{2}}+v_{\mathrm{1}}^{\mathrm{2}}}}},\text{(A8)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{yz}^{\text{w}}=-u^{*\mathrm{2}}{\displaystyle \frac{v_{\mathrm{1}}}{\sqrt{u_{\mathrm{1}}^{\mathrm{2}}+v_{\mathrm{1}}^{\mathrm{2}}}}}.\end{array}$$

From Eqs. (A6) to (A8), it is clear that a large LOTW mismatch, especially at the first cell, leads to an error in the wall shear stress and, in turn, in the vertical profile of resolved shear stress. While the LOTW mismatch only causes a departure of the wind profile from the logarithmic law of the wall, a mismatch in the shear stress profile is much more serious as it affects the turbulence intensity level experienced by the wind farm and possibly turbine and wind farm wake recovery. Hence, the wall model has to be modified such that the correct wall shear stress is applied regardless of the employed grid resolution.

Appendix B: Precursor simulations

Figure B1 Statistics gathered from neutral ABL simulations having grid spacings of $\mathrm{15}\times \mathrm{15}\times \mathrm{10}$, $\mathrm{50}\times \mathrm{50}\times \mathrm{10}$, and $\mathrm{50}\times \mathrm{50}\times \mathrm{30}$ m. The modifier u^* in the legend entries identifies those cases where the friction velocity in the wall model has been set equal to the $\mathrm{15}\times \mathrm{15}\times \mathrm{10}$ m case instead of being evaluated using Eq. (A6). All cases are normalized with u_ref=9 m s⁻¹ and $u^{*}=\mathrm{0.297}$ m s⁻¹ (the latter corresponds to the case with mesh resolution of $\mathrm{15}\times \mathrm{15}\times \mathrm{10}$ m). Panels (a) and (b) show the velocity magnitude and shear stress profile; (c) and (d) depict the LOTW scaling and non-dimensional shear (theoretical laws are identified by the dashed black lines); (e), (f), and (g) report the mean velocity variances; and (h) reports the wind veer resulting from the Coriolis force.

A precursor simulation that satisfies the criteria expressed by Eqs. (A1) to (A3) is first conducted on a fine grid, and the resulting friction velocity $u_{\text{fine}}^{*}$ is calculated using Eq. (A6). Then, a second precursor characterized by a coarser mesh where Δx_coarse>2Δx_fine and Δy_coarse>2Δy_fine is conducted and the wall model is modified such that Eqs. (A7) and (A8) are used with $u^{*}=u_{\text{fine}}^{*}$. This essentially renders the wall model independent of the employed grid size and ensures matching of the shear stress profile.

To verify this approach and to show the effects of coarsening the precursor simulation grid, we conducted simulations of a fully neutral atmospheric boundary layer using three different values of the grid spacing. The finer grid uses a resolution of $\mathrm{15}\times \mathrm{15}\times \mathrm{10}$ m, which is expected to satisfy all Brasseur and Wei (2010) criteria. The coarser cases use a resolution of $\mathrm{50}\times \mathrm{50}\times \mathrm{10}$ and $\mathrm{50}\times \mathrm{50}\times \mathrm{30}$ m. The former does not satisfy $\mathcal{R}/{\mathcal{R}}^{*}>\mathrm{1}$, while ${\mathit{Re}}_{\mathrm{LES}}/{\mathit{Re}}_{\mathrm{LES}}^{*}\approx \mathrm{1}$. However, a sufficient number of vertical grid nodes is used to resolve the ABL vertically. When the vertical grid spacing is increased, the first two criteria are satisfied, but $N_{\mathit{\delta }}/N_{\mathit{\delta }}^{*}<\mathrm{1}$. Both coarser cases are simulated using the conventional wall model and the modified wall model, where $u_{\text{fine}}^{*}$ is calculated from the precursor characterized by the finer grid, leading to five precursor cases in total.

As no stratification is present, the domain size is set to $\mathrm{4.2}\times \mathrm{4.2}\times \mathrm{0.7}$ km, essentially fixing the ABL height H to 700 m. The Coriolis parameter f_c is set to $\mathrm{1.184}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻¹, corresponding to a latitude of 54.5°. The equivalent roughness height z₀ is set to 0.001 m. Horizontal boundaries are periodic, while a stress-free condition is applied at the upper boundary. At the bottom, wall shear stress is directly applied in the momentum equation using Eqs. (A7) and (A8), while velocity at the ghost nodes is set such that the wall normal gradient at the boundary face coincides with the one evaluated at the first cell center. This boundary condition allows the velocity to be ultimately determined by the amount of shear stress applied by the wall model. Moreover, since the entire wall shear stress is modeled, the effective viscosity is set to zero at the wall to avoid double counting. A uniform driving pressure gradient is applied to the momentum equation such that the horizontally averaged velocity at h_ref=90 m is equal to 9 m s⁻¹ (Stipa et al., 2024b). All simulations use the dynamic Smagorinsky model (Lilly, 1992) with the Lagrangian averaging of the model coefficient proposed by Meneveau et al. (1996). Each simulation is carried out for 100 000 s, and statistics are horizontally and time averaged from 80 000 to 100 000 s.

Figure B1 compares zero- and first-order statistics from the five precursor simulations. Firstly, it is evident how the simulation in general depends on the grid size when Brasseur and Wei (2010) criteria are not satisfied. The case depicting the largest deviation from the finest precursor corresponds to the $\mathrm{50}\times \mathrm{50}\times \mathrm{10}$ m grid resolution without wall model correction. This case shows some difference in the mean velocity profile, produced by the LOTW mismatch, and, more importantly, it completely fails to capture the vertical shear stress profile and also strongly underestimates velocity variances. When the correct friction velocity is applied to the wall model, things improve substantially. The shear stress profile closely matches the finer case, and velocity variances also improve. The error in the mean velocity profile is reduced, but a LOTW mismatch is still observed. When the vertical mesh resolution is decreased from 10 to 30 m, things improve when employing the conventional wall model, and the shear stress almost matches with the fine precursor case. In fact, reducing the aspect ratio also reduces the LOTW mismatch according to Brasseur and Wei (2010), and so applying the correct friction velocity does not have a large impact on the flow statistics. Surprisingly, the wind veer seems to show little sensitivity to the wall model and the grid size except when very close to the wall.

These results suggest that a precursor simulation characterized by a very large grid spacing – thus not designed to fulfil the LOTW matching criteria – can still capture the shear stress profile when the correct friction velocity is imposed at the wall. Moreover, a small deviation from the logarithmic profile, located at the boundary layer top, is observed when the velocity profile is compared to that of an LES satisfying LOTW scaling. Therefore, results from both the $\mathrm{50}\times \mathrm{50}\times \mathrm{30}$ m case (with and without correction on u^*) and the $\mathrm{50}\times \mathrm{50}\times \mathrm{10}$ m case (only with u^* correction) are deemed suitable for providing an inflow condition to those wind farm simulations which, because they employ the AFM, are characterized by a comparably coarse grid.

Appendix C: Rotor center position relative to the grid

In this section, we analyze the AFM's sensitivity to the relative position between the rotor center and the surrounding grid points. Previously in Sect. 3, the rotor center position coincided with a cell vertex in the horizontal plane, meaning that the mesh arrangement was perfectly symmetrical with respect to the rotor (see Fig. C1). In this analysis, we horizontally shift the turbine location, leaving the grid unchanged and addressing the effect that this has on the velocity in the wake and on the turbine thrust and power. Specifically, we use the coarsest grid employed in Sect. 3, i.e. a horizontal resolution of 60×60 m, and we systematically displace the wind turbine by 20 and 40 m both in the streamwise and spanwise directions, leading to the five cases shown in Fig. C1. These correspond to a somewhat random placement of the wind turbine relative to the surrounding grid, similar to what is observed for large-scale wind farm computations. Notably, the chosen grid resolution, where $D/\mathrm{\Delta }\approx \mathrm{2}$, is the coarsest limit for employing the AFM. Consequently, results from the following analysis are to be considered as a worst-case scenario.

Figure C1 Top view of rotor position relative to the surrounding grid for the five cases conducted in this section. Panel (a) corresponds to the reference case already conducted in Sect. 3. Shift values are indicated at the top left of each panel.

Table C1 Turbine thrust and power, averaged over the last 250 s, for the cases conducted in Appendix C. The first two columns on the left indicate the streamwise and spanwise shift, respectively. The first line of the table corresponds to the reference case with no shift. Differences in thrust and power, in percent of the reference case, are also reported, together with the error in projecting the thrust force, calculated as $|T_p-T|/T$ %, where $T_p=\int bdV$ is the integrated body force b after projection and T is the turbine thrust before projection.

The numerical setup of the simulations is the same as Sect. 3, where $r_{\mathrm{1}/\mathrm{2}}=R$, s=6, and $\mathit{\sigma }/\mathrm{\Delta }x=\mathrm{1.5}$ have been set for the AFM. Figure C2 shows the spanwise velocity profiles at the turbine location and at three diameters downstream. In particular, while the wake shape exhibits some degree of asymmetry in the spanwise direction at both locations, the wake magnitude already loses the dependency on the streamwise shift at three diameters downstream. In general, it can be stated that the horizontal shift does not affect the magnitude of the velocity deficit in the far wake, while it slightly shifts the wake centerline in the spanwise direction. However, this effect is minimal, and the model seems pretty robust considering that only two cells are used to resolve the wake. In fact, when simulating an entire wind farm, as for example in Sect. 4, where the turbine position is inevitably random with respect to the grid, the model is able to correctly predict the wind speed deficit in the wake of the entire cluster, even when a very coarse resolution is employed.

Figure C2 Spanwise velocity profiles for each of the cases conducted in Appendix C, sampled at the rotor position (a) and three diameters downstream (b). Data points are interpolated in the streamwise direction between the two closest streamwise cell center coordinates and plotted at the spanwise cell center coordinates. Data are shifted along the spanwise direction by the opposite of the turbine shift so that the ideal wake centerline coincides in all cases.

Finally, Table C1 reports the thrust and power, their percentage difference with respect to the un-shifted case, and the error in projecting the thrust force for each individual case. As can be noticed, the difference in thrust and power due to the turbine shift is always below 5 % of the reference case, which is far less than the difference obtained when using erroneous model parameters, as reported in Fig. 4.

Overall, this analysis shows that the AFM has little sensitivity to the turbine position relative to the grid. Moreover, such sensitivity is expected to further decrease when the ratio $D/\mathrm{\Delta }$ is increased. For instance, as stated in Sect. 7, two to three cells are suggested when modeling an upstream wind farm, while four cells are required for thrust and power data comparable to the ADM's.