the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
A new RANSbased wind farm parameterization and inflow model for wind farm cluster modeling
Oscar GarcíaSantiago
Mark Kelly
Alexander Meyer Forsting
Camille DubreuilBoisclair
Knut Sponheim Seim
Marc Imberger
Alfredo Peña
Niels Nørmark Sørensen
PierreElouan Réthoré
Offshore wind farms are more commonly installed in wind farm clusters, where wind farm interaction can lead to energy losses; hence, there is a need for numerical models that can properly simulate wind farm interaction. This work proposes a Reynoldsaveraged Navier–Stokes (RANS) method to efficiently simulate the effect of neighboring wind farms on wind farm power and annual energy production. First, a novel steadystate atmospheric inflow is proposed and tested for the application of RANS simulations of large wind farms. Second, a RANSbased wind farm parameterization is introduced, the actuator wind farm (AWF) model, which represents the wind farm as a forest canopy and allows to use of coarser grids compared to modeling all wind turbines as actuator disks (ADs). When the horizontal resolution of the RANSAWF model is increased, the model results approach the results of the RANSAD model. A double wind farm case is simulated with RANS to show that replacing an upstream wind farm with an AWF model only causes a deviation of less than 1 % in terms of the wind farm power of the downstream wind farm. Most importantly, a reduction in CPU hours of 75.1 % is achieved, provided that the AWF inputs are known, namely, wind farm thrust and power coefficients. The reduction in CPU hours is further reduced when all wind farms are represented by AWF models, namely, 92.3 % and 99.9 % for the double wind farm case and for a wind farm cluster case consisting of three wind farms, respectively. If the wind farm thrust and power coefficient inputs are derived from RANSAD simulations, then the CPU time reduction is still 82.7 % for the wind farm cluster case. For the double wind farm case, the RANS models predict different wind speed flow fields compared to output from simulations performed with the mesoscale Weather Research and Forecasting model, but the models are in agreement with the inflow wind speed of the downstream wind farm. The RANSADAWF model is also validated with measurements in terms of wind farm wake shape; the model captures the trend of the measurements for a wide range of wind directions, although the measurements indicate more pronounced wind farm wake shapes for certain wind directions.
The growth of offshore wind energy has led to wind farm clustering, where wind farm interaction is unavoidable. Recently, the Danish government released a report with plans for a new 10 GW offshore wind farm cluster situated around an artificial energy island hub in the North Sea (COWI, 2020). This wind farm cluster consists of 10 wind farms of 1 GW, with a wind farm interdistance of only 8 km using a nonoptimized wind farm cluster layout. More examples of wind farm clusters can be found in other parts of the North Sea, the Baltic Sea, and the East China Sea (4coffshore.com, 2022). To our best knowledge, there are currently no International Electrotechnical Commission (IEC) standards regarding recommended minimal distances between wind farms to avoid losses due to wind farm wakes, and wind farm owners do not have any control over potential newly built neighboring wind farms. Hence, there is a need for improved numerical models that can estimate wind farm wake losses (and potential gains due to speed up effects) in order to establish standards for wind farm placement in wind farm clusters.
Lowfidelity engineering wind turbine wake models (Göçmen et al., 2016) can be used to model wind farm interaction; however, large model uncertainties exist due assumptions on wake shape and wake superposition methods. The performance of these models is case dependent and can be poor (Fischereit et al., 2022) or reasonable (Nygaard et al., 2020), which often depends on how the models are calibrated. The current stateoftheart numerical models employed to simulate wind farm clusters are based on mesoscale models, e.g., the Weather Research and Forecasting (WRF) model (Skamarock et al., 2019). The WRF model is normally run with nested domains, where the finest domain has a typical horizontal cell length of 1 km. The effect of a wind farm is modeled by a simple wind farm parameterization, where a drag force is added to the momentum equation (Fitch et al., 2012; Volker et al., 2015; Abkar and PortéAgel, 2015), and its magnitude is based on the wind turbine thrust curve. Fitch et al. (2012) also included a source of turbulent kinetic energy based on the windturbineextracted kinetic energy that is not converted to electricity. A turbulent kinetic energy source term was also motivated by Abkar and PortéAgel (2015) based on dispersive Reynolds stresses due to under resolving the windfarminduced turbulent kinetic energy. Volker et al. (2015) did not include such a source term but take into account subgridscale vertical wake expansion within one grid cell.
It is not trivial to verify the wind farm parameterizations in mesoscale models because, e.g., in the WRF model the wind farm parameterizations are implemented within 1D planetary boundary layer (PBL) schemes (Peña et al., 2022). These schemes are not scale aware and become horizontal resolution dependent when decreasing the grid size below ≈1 km. Only until recently, a new fully 3D PBL scheme was implemented in the WRF model (Kosović et al., 2020), and the Fitch parameterization was also already configured to run with it (Rybchuk et al., 2022). In this work, we propose a RANSbased (Reynoldsaveraged Navier–Stokes) wind farm model, the actuator wind farm model (AWF), which can be used to verify mesoscale wind farm parameterizations in a microscale environment because the minimal horizontal spacing is not limited. The AWF represents a wind farm as a forest canopy using distributed drag forces. The use of a forest canopy model to represent a single wind turbine wake has been employed by Réthoré (2009), although the model performed poorly compared to a highfidelity turbulenceresolving simulation. In the present work, we will show that the use of a forest canopy model for the entire wind farm can work quite well, as long as the correct wind farm drag force magnitude and distribution is applied. This is achieved by applying a precalculated wind farm drag force magnitude that is both wind speed and wind direction dependent and by employing a horizontal drag force distribution as a superposition of twodimensional Gaussian functions centered at the wind turbine locations. The latter represents a smoothed number of wind turbines per cell as opposed to counting the number of wind turbines per cell, as commonly used in WRF's wind farm parameterization, in order to overcome aliasing effects that can lead to artificial wind farm wake shapes. It is common to use Gaussian functions to distribute forces of simplified wind turbine models in computational fluid dynamics (CFDs), as performed for actuator line (Sørensen and Shen, 2002), sector (Storey et al., 2015), and disk (AD) (Mikkelsen, 2003; Wu and PortéAgel, 2011) models. The AWF model can also be used to calculate the annual energy production of a large wind farm cluster, which would become computationally very expensive if all wind turbines are modeled by ADs that require a horizontal cell size in the order of $D/\mathrm{8}$, where D is the rotor diameter. If the effect of neighboring wind farms on a single wind farm needs to be simulated, the wind farm of interest can be modeled as ADs using a fine horizontal spacing, while all the other wind farms can be modeled as AWFs using a much coarser horizontal spacing; we refer to this model as RANSADAWF. Such a numerical setup is only slightly more expensive in terms of computational costs compared to simulating the wind farm of interest without the neighboring wind farms. A further reduction in computational effort can be achieved by modeling all wind farms as AWFs, hereafter known as the RANSAWF model.
Wind farm (cluster) modeling with RANS is challenging when an inflow model with an atmospheric boundary layer (ABL) height is employed together with a large horizontal domain in the order of 200 km and more (van der Laan et al., 2017). This is because the eddy viscosity above the ABL height is very small, and numerical instabilities can appear if the domain is large enough. These numerical instabilities could be interpreted as numerical gravity waves because a solution with standing waves can be obtained when flow is solved with an unsteady RANS method. One could employ a more robust inflow model representing a neutral surface layer inflow; however, such a model is not very realistic for either large wind turbines that are expected to operate outside the surface or a wind farm cluster where the effect of the Coriolis force can become important (van der Laan and Sørensen, 2017b). An ABL inflow model can be obtained when using the global turbulence lengthscale limiter of Apsley and Castro (1997), since it does not require a temperature equation, and the ABL height is determined by setting a maximum turbulence length scale. In addition, the model is only dependent on two nondimensional numbers (van der Laan et al., 2020), which eases getting the desired inflow profile from a library of precalculated ABL profiles. However, the global turbulence lengthscale limiter also limits all turbulence length scales in the threedimensional domain as wind farm wake turbulence length scales. The latter can lead to a wind farm wake that stops recovering with downstream distance, which is a nonphysical result. Furthermore, the effective inversion strength is implicit in this model, and is relatively strong, which can escalate the problem of numerical instabilities in large horizontal domains. In the present work, we propose an alternative inflow model for wind farm cluster simulations in RANS for neutral conditions near the ground but with a stable inversion layer; such a model reflects a conventional neutral ABL as commonly applied in largeeddy simulations (Allaerts and Meyers, 2018; Kelly et al., 2019; Liu et al., 2021). The model does not use the global turbulence lengthscale limiter of Apsley and Castro (1997). Instead, the ABL height is set by a prescribed analytical temperature profile that includes an inversion height and strength. A temperature equation is not solved in order to maintain a steadystate inflow. Note that the use of an active temperature equation leads to an unsteady RANS model, and we are interested in a steadystate RANS model.
The wind farm cluster validation case and the corresponding supervisory control and data acquisition (SCADA) measurements are discussed in Sect. 2. The numerical setup of the RANS simulations using ADs, the proposed AWF model, and the new RANS inflow model are discussed in Sect. 3.1.1–3.1.3, respectively. The numerical setup of the WRF simulations and engineering wake model calculations are presented in Sect. 3.2 and 3.3, respectively. The AWF model is verified and compared with the AD model in Sect. 4.1. In Sect. 4.2, the AWF model is compared with outputs from realtime simulations using the WRF model (and a wind farm parameterization) and the TurbOPark engineering wake model (Nygaard et al., 2020; Pedersen et al., 2022), and validated with SCADA measurements of a wind farm cluster consisting of three wind farms: two are operated by Equinor, the other one by Ørsted.
Figure 1 depicts the investigated offshore wind farm cluster, which is located in the North Sea near the east coast of the United Kingdom. There are four wind farms located in an area of about 70×30 km^{2}. SCADA measurements are made available for the most eastern wind farm, Dudgeon, and wind turbine power data from the front row wind turbines are employed to probe the wind farm wakes of Sheringham Shoal and Race Bank, for wind directions around 235 and 270^{∘}, at wind farm interspacings of 16 and 26 km, respectively. The most western wind farm, located 54 km away from Dudgeon, consists of several smaller wind farms (Inner Dowsing, Lynn and Lincs), which are not included in the present study. More details on the investigated wind farms are listed in Table 1, where A_{wf} is wind farm area, computed by a concave hull method; N_{wt} is the number of wind turbines; and D and z_{H} are wind turbine rotor diameter and hub height, respectively. The mean wind turbine interspacing is calculated following Sørensen and Larsen (2021): $\sqrt{{A}_{\mathrm{wt}}}/\left(D\right[\sqrt{{N}_{\mathrm{wt}}}\mathrm{1}\left]\right)$. The wind turbines are propriety to Siemens Gamesa Renewable Energy (SGRE), and details on the thrust coefficient and power curves cannot be published. In the present work, we investigate a belowrated inflow wind speed, and the thrust curves of both wind turbines have typical belowrated thrust coefficient values.
2.1 Measurements
This study uses 3 years of SCADA data from Dudgeon, from 1 January 2018 to 1 January 2021. For each turbine, the data are first averaged to 10 min periods, and then each period is kept only if it does not contain any missing data, nonproduction or curtailment periods, low power values (<100 kW), or low wind speed values (<3 m s^{−1}). Following the filtering for each turbine, the 10 min period is kept if at least 64 out of 67 turbines remain, leading to a final availability of 56 % of the data for the whole period.
In parallel, ERA5 data interpolated from the nearest ERA5 grid points to the wind farm location, for the same 3year period, are used to estimate the Obukhov length L from the wind speed at a height of 10 m, the temperature at a height of 2 m, and the sea surface temperature, applying an iterative method following Ott and Nielsen (2014). The data are divided into three atmospheric stability classes, namely unstable (−200 m < L < 0 m), neutral ($\leftL\right\ge \mathrm{200}$ m), and stable (0 m ≤ L < 200 m). This classification is then applied to the production data to get three data sets. The data of each stability class are then binned by 1 m s^{−1} and 5^{∘} based on a reference wind speed and direction calculated by averaging the front row wind turbines for a specific sector. The latter is performed due to a lack of concurrent inflow measurements. The fact that we use the same wind turbine row to define the freestream and to probe the wake of an upstream wind farm means that is not possible to use the data set to validate the simulations results for the magnitude of wind farm wake losses. Therefore, the validation in Sect. 4.2 is only performed for the wind farm wake shape, where the wind speed of the front row wind turbines is normalized with the wind speed averaged over the same turbines.
3.1 RANS
The wind farm simulations are performed with PyWakeEllipSys v3.0 (DTU Wind Energy, 2022b), which is a Python wrapper of the inhouse flow solver EllipSys3D. EllipSys3D is a finitevolume CFD code that was initially developed by Michelsen (1992) and Sørensen (1994). The RANS model of EllipSys3D is used to simulate the steadystate wind farm flow under neutral atmospheric conditions including an ABL height, Coriolis forces, and a rough wall boundary with uniform roughness, representing homogeneous terrain. The wind turbine forces are represented by two different methods. The baseline approach is an AD model (Réthoré et al., 2014) and is used to verify the proposed AWF method. Each method is discussed separately in Sect. 3.1.1 and 3.1.2, respectively.
In this work, different RANS flow domains are used to model wind farm clusters containing single, double, and triple wind farms. A sketch of the flow domain types of the different applied actuator methods is depicted in Fig. 2. All flow domains are Cartesian grids, where the inflow direction at the reference height is aligned with the x axis, and different wind directions are modeled by rotating the wind farm cluster while maintaining the grid and the global inflow direction. The RANSAD method (Fig. 2a) represents all wind turbines in a wind farm by ADs. Each wind turbine is treated as a single AD model, which has its own polar grid that is connected to the flow domain and is also controlled independently. The cells around the ADs, as marked by the cyan box in Fig. 2a, are uniformly spaced in the horizontal plane, using a fine resolution in order to resolve the wind turbine wakes. The RANSAWF method (Fig. 2b) follows a similar flow domain as the RANSAD method; however, each wind farm is considered as a single AWF model that uses its own force controller. The wind turbine forces are distributed in a Cartesian grid that encapsulates the entire wind farm, and this Cartesian grid is then connected to the flow domain grid. The refined area in the RANSAWF method can be an order of magnitude coarser compared to the refined area for the RANSAD model, which is investigated in detail in Sect. 4.1. The third domain type is depicted in Fig. 2c and represents the RANSADAWF method, where both AD and AWF models are present. Since the AWF model does not require the fine spacing of the AD models, a second region in the flow domain is defined, as marked by the magenta box in Fig. 2c, where the cells are expanded in the horizontal direction while moving away from the cyan box, up to a maximum set spacing. It should be noted that the three methods as depicted in Fig. 2 can also be applied to simulate multiple wind farms. In this case, the size of the refined inner domain(s) is adjusted to resolve all wind turbine and farm wakes. An overview of the all the applied domains for the validation cases, as well as the computational effort, is provided in Table 4.3.
Figure 3 depicts a detailed description of the wind farm cluster validation case, where Dudgeon is modeled as ADs, and the other two wind farms, Sheringham Shoal and Race Bank, are modeled as AWFs (the RANSADAWF flow domain type of Fig. 2c). The domain includes two areas with a refined horizontal spacing as marked by the cyan and magenta rectangles, as shown in Fig. 3a. The cyan rectangle represents a uniform horizontal spacing of ${D}_{\mathrm{ref}}/\mathrm{8}$ around the Dudgeon wind farm (van der Laan et al., 2015c) and includes a sufficiently large area to resolve the wind farm flow for the wind directions of interests, namely, 185–200^{∘}. Here, D_{ref} is the reference rotor diameter used to scale the grid, and it is set to the smallest rotor diameter of all the wind turbines (107 m). The magenta rectangle includes cells that are stretched towards the boundary, and the maximum horizontal cell size is limited to 2 D_{ref}, which is the same as the horizontal resolution of the AWF model. The size of this area is set to cover both Sheringham Shoal and Race Bank for the wind directions of interest. Outside the magenta rectangle, the cells are further stretched, with a maximum expansion ratio of 1.2 covering a distance of 500 D_{ref} (53.5 km). The total number of cells in the x and y directions are 1920 and 2048, respectively. In the vertical, as depicted in Fig. 3b, the cells are stretched while moving away from the ground, and the first cell height is set to ${D}_{\mathrm{ref}}/\mathrm{200}$ (0.54 m). The maximum cell size in the first layer (up to z=3D_{ref}) is set to ${D}_{\mathrm{ref}}/\mathrm{6}$; however, the cells in the rotor area and below are smaller than ${D}_{\mathrm{ref}}/\mathrm{8}$. The height of the domain is set to 10 D_{ref}, and a total number of 64 cells are used in the vertical direction. This results in 252 million cells, which we divide into 960 blocks using a block edge size of 64 cells. We would have preferred to set a taller domain to avoid numerical blockage; however, there is a convergence issue when a tall domain height is applied together with an ABL inflow model that includes a low mixing region above the ABL height, which is further discussed in Sect. 3.1.3.
Figure 3a and b also depict the boundary conditions, where inlet conditions are enforced over the inflow and upper domain faces; consequently, the flow is lid driven. Periodic boundary conditions are used at the lateral boundaries to be able to include wind veer. An outlet boundary condition is used at the outflow boundary, at which all gradients in the normal direction are assumed to be zero. The ground is modeled by a rough wall boundary condition following Sørensen et al. (2007).
3.1.1 AD model
The RANSAD wind farm simulation setup is similar to the one used in a previous work (van der Laan et al., 2015b). The wind turbines are modeled as ADs, and each AD is represented by a polar grid using 10×32 cells in the radial and azimuthal directions, respectively. Each cell of the polar grid is connected to a set of flow domain cells, as discussed in detail by Réthoré et al. (2014). For each iteration, the thrust and tangential forces are calculated on the polar grid and are then injected in the flow domain grid as momentum sinks. The local velocities from the flow domain grid are returned to the polar grid to recalculate the AD forces and evaluate the wind turbine power. The magnitude of the thrust and tangential forces are controlled by a lookup table of alternative thrust and power coefficients that are based on a diskaveraged velocity to avoid the need for a freestream wind speed that is typically not known for interacting wind turbines; more details can be found in previous works (van der Laan et al., 2015a, 2019). For the model verification (Sect. 4.1), a general wind turbine model (van der Laan et al., 2022) is used with the same rotor diameter and hub height as the wind turbine from Dudgeon (SGRE 6 MW) but with a typical belowrated thrust coefficient of 0.8, a power coefficient of 0.45, and a tipspeed ratio of 8, which are all kept constant up to a rated wind speed of 10 m s^{−1}. The AD force distribution is modeled by the analytical AD force model from Sørensen et al. (2020). For the model validation (Sect. 4.2), we use an AD with a prescribed normalized thrust force distribution that is rescaled to obtain the desired thrust force magnitude (van der Laan et al., 2015a). For the SGRE 3.6 and 6 MW turbines, we use the thrust force distributions of the NREL5 MW (Jonkman et al., 2009) and DTU10 MW (Bak et al., 2013) reference wind turbines, respectively, calculated with previously performed detachededdy simulations (Réthoré et al., 2014; Bak et al., 2013) for a belowrated wind speed of 8 m s^{−1}. This is because we lack information on the rotational speed of the SGRE wind turbines, which is a required input for the analytical rotor model of Sørensen et al. (2020), and we are lacking the corresponding airfoil data that are needed for a higherfidelity AD model based on a blade element momentum method. However, we do not expect a large influence on the wind turbine and farm wakes when using the AD model based on normalized thrust force distribution, as opposed to using AD models based on the analytical rotor model or airfoil data, because the prescribed normalized thrust force distribution takes into account a realistic radial force distribution, and it employs the provided thrust coefficient curves of the SGRE wind turbines.
3.1.2 AWF model
The AWF model represents each wind farm as a single entity, and its effect on the flow is modeled by a distributed thrust force:
where ρ is the air density, C_{T,wf} is the wind farm thrust coefficient, $A(x,y,z)$ is the wind farm force distribution representative of the wind turbine density, and U_{i} is the local velocity vector. Equation (1) is the same as that used to model a forest canopy drag force using ${C}_{\mathrm{d}}={C}_{\mathrm{T},\mathrm{wf}}/\mathrm{2}$. In fact, we employ the forest canopy model from EllipSys3D as used in previous works (Boudreault, 2015; Dellwik et al., 2019), but here we employ the same turbulence model as in the RANSAD setup as opposed to the turbulence model developed for forest canopies (Sogachev et al., 2012; Boudreault, 2015; Dellwik et al., 2019).
One could employ additional terms in the turbulence model equations to account for the effect of under resolving the wind farm layout in the AWF model when using large horizontal cells (Abkar and PortéAgel, 2015). However, the literature is divided about whether this extra term should be zero (Volker et al., 2015), act as a source of turbulent kinetic energy (Fitch et al., 2012; Abkar and PortéAgel, 2015), or act as a sink of turbulent kinetic energy (Sogachev et al., 2012). Investigating this is not in the scope of this study, and so we do not use a source term of turbulence, partly because we already get reasonable results with only a momentum sink (Eq. 1).
AWF model force distribution $A(x,y,z)$
The drag force of the AWF model $A(x,y,z)$ is distributed over a Cartesian grid encapsulating the wind farm. The vertical drag force distribution represents the rotor area slices following Abkar and PortéAgel (2015). The horizontal drag force distribution represents the wind turbine density. Currently, when running the wind farm parameterization in the WRF model, the user sets the number of turbines per grid cell (for idealized cases) or inputs the turbine locations (realtime forcing cases), and so the number of turbines within the same grid cell are counted. Such counting can lead to griddependent results, particularly for large horizontal grid spacing and artificial velocity deficit shapes due to aliasing effects of the wind farm layout, as the WRF model is typically run at a horizontal grid spacing larger than 1 km. An example of this issue is illustrated in Sect. 4.1. To overcome this issue, we propose an alternative method to determine the horizontal drag force distribution of the wind farm by representing each wind turbine position as twodimensional Gaussian functions f(x,y) instead of points:
with x_{wt} and y_{wt} as the wind turbine positions and σ_{x} and σ_{y} as the standard deviations in the streamwise x and lateral y directions, respectively. We use σ_{x}=2Δ and ${\mathit{\sigma}}_{y}=\mathrm{max}(\mathrm{2}\mathrm{\Delta},D/\mathrm{4})$, with Δ as the horizontal grid spacing in the AWF model. The AWF horizontal grid spacing Δ does not have to be the same as the horizontal spacing in the CFD grid δ. The use of 2Δ is motivated in Sect. 4.1, and the max limiter is used in order to approach an AD model and avoid a tooconcentrated force; the latter could otherwise cause numerical problems. The AWF horizontal force distribution for each cell in the AWF model is obtained by superposing Eq. (2) for all wind turbines. This means that for a particular cell in the AWF grid (i,j), the wind turbines outside this cell can also contribute to the drag force, and cells that do not include a turbine can have nonzero force. In general, a smoothed wind farm layout is obtained as a smoothed wind farm force distribution, and the individual wind turbine positions are resolved when a fine enough horizontal spacing is set. Another property of the Gaussian method is that a regular wind farm layout will always have a uniform thrust force distribution in the middle of the farm, which is not guaranteed when binning the number of wind turbines per cell. The magnitude of the drag distribution $A(x,y,z)$ is not relevant because we calibrate the wind farm thrust coefficient to get the desired total wind farm thrust force, as discussed in the following section. The vertically integrated AWF horizontal force distribution of Race Bank and Sheringham Shoal are depicted in Fig. 4 for two different horizontal resolutions: Δ=D and Δ=2 D. Figure 4 shows how the regular wind farm layout of Sheringham Shoal results in a uniformly distributed drag force when Δ is coarsened from D to 2 D by comparing Fig. 4b and d. In addition, values below 0.01 of the maximum force distribution are set to zero for numerical efficiency. The distributed drag forces are stored in a separate Cartesian grid for each AWF model. The drag force of a CFD cell that lies within an AWF model is obtained by trilinear interpolation.
AWF model force magnitude C_{T,wf}
Equation (1) depends on the local velocity and a local wind farm thrust coefficient, while a single wind turbine modeled as an AD with fixed forces can be set by a known thrust coefficient based on the freestream wind speed. Dellwik et al. (2019) calibrated the drag coefficient of a single tree numerically, using a measured drag force as input. A related method was introduced by Abkar and PortéAgel (2015), where the thrust force was multiplied by a correction parameter obtained by largeeddy simulations to correct the cell wind speed in a wind farm parameterization model, to account for different wind farm layouts. For the AWF model, we follow a similar approach as Dellwik et al. (2019), where C_{T,wf} is calibrated using simulations of a single AWF model in order to obtain the wind farm thrust force determined by precursor RANSAD simulations of the same wind farm. The wind farm thrust force depends on both the inflow wind speed through the wind turbine thrust coefficient curves and the inflow wind direction because of wind turbine interaction. Hence, C_{T,wf} is a function of both inflow wind speed and wind direction. For wind farm cluster simulations, the freestream wind direction and wind speed are undefined. Therefore, we employ a force controller for each AWF model, where C_{T,wf} is updated every solver iteration based on an interpolation of a lookup table, where C_{T,wf} is a function of the AWF volumeaveraged wind speed and wind direction obtained from the AWF volumeaveraged velocity vector U_{AWF,i}:
U_{AWF,i} is weighted by the drag force distribution A, so it does not include the buffer zone around the forest canopy and accounts only for the volume around the wind turbine position for fine spacing AWF models. Without the drag force distribution weighting, U_{AWF,i} is typically over predicted. A similar approach was followed by Churchfield et al. (2017) for the application of the effective wind speed calculation for an actuator line model. The force control lookup table is created from the single AWF model calibration simulations used to determine C_{T,wf} from a desired total wind farm thrust force. The simulation steps necessary to perform the wind farm cluster validation case simulations are summarized below.

For each wind farm modeled as an AWF,
 a.
RANSAD simulations are performed to calculate the total wind farm thrust force for a range of wind speed and directions ($U=\mathrm{7},\mathrm{8},\mathrm{9}$ m s^{−1}, $\mathit{\varphi}=\mathrm{175}\mathrm{310}$^{∘} with a 5^{∘} interval);
 b.
RANSAWF simulations are performed using the same flow cases as step 1a in order to calculate C_{T,wf} and C_{P,wf} as a function of the AWF volumeaveraged wind speed and wind direction.
 a.

RANSADAWF or RANSAWF cluster simulations are performed.
In steps 1 and 2, we neglect the influence of inhomogeneous inflow conditions on the wind farm thrust and power, as for example partial wind farm wake effects of the neighboring wind farms. However, the AWF model (as applied in step 2 with a force controller) can partly respond to inhomogeneous conditions because the local thrust forces are dependent on the local velocity, although variations in wind turbine thrust coefficients are not captured due to the use of a global wind farm thrust coefficient. The impact of this assumption is investigated in Sect. 4.2.1. Step 1a is computationally expensive, especially when many different wind farms in a wind farm cluster are modeled as AWFs, and each of them requires input from precursor RANSAD wind farm simulations. The computational costs of step 1a could be alleviated by running RANSAD wind farm simulations more efficiently (van der Laan et al., 2022). For example, one could run consecutive wind direction cases (applicable to homogeneous terrain) and consecutive wind speed cases (applicable to inflow models that include Reynolds number similarity) both to save solver iterations, and employ wind farm layout symmetry to reduce the number of wind direction cases. In addition, one could also employ a fast engineering wake model to calculate the total wind farm thrust force, although this would most likely introduce a higher model uncertainty. The AWF model input can also be calculated with a RANSbased surrogate wind farm model, which is currently investigated in a followup work (van der Laan et al., 2023).
3.1.3 Atmospheric inflow and turbulence models
We employ a two equation turbulence model in the form of a k–ε model (Launder and Spalding, 1974):
with x_{j} as the Cartesian coordinates, $\mathit{\nu}=\mathrm{1.78406}\times {\mathrm{10}}^{\mathrm{5}}$ m^{2} s^{−1} as the molecular viscosity of air, ν_{T} as the turbulent eddy viscosity, k as the turbulent kinetic energy, and ε as the dissipation of k. In addition, 𝒫 is the mechanical production of turbulence, ℬ is the buoyant turbulence production or destruction, and S_{k,amb} and S_{ε,amb} are ambient turbulence sources. The remaining unknowns in Eqs. (4)–(6) (except f_{P}) are turbulence model coefficients, here set as constants $({C}_{\mathit{\mu}},{\mathit{\sigma}}_{k},{\mathit{\sigma}}_{\mathit{\epsilon}},{C}_{\mathit{\epsilon},\mathrm{2}})=(\mathrm{0.03},\mathrm{1.0},\mathrm{1.3},\mathrm{1.92})$. ${C}_{\mathit{\epsilon},\mathrm{1}}={C}_{\mathit{\epsilon},\mathrm{2}}{\mathit{\kappa}}^{\mathrm{2}}/\left({\mathit{\sigma}}_{\mathit{\epsilon}}\sqrt{{C}_{\mathit{\mu}}}\right)$ is used to enforce a balance with a logarithmic wind profile with κ=0.4 as the von Kármán constant, and C_{ε,3} is discussed later. In addition, the well known Boussinesq (1897) hypothesis is applied to define the relationship between the Reynolds stresses and the strainrate tensor. In absence of the buoyancy and without ambient sources, the model represents the k–ε–f_{P} turbulence model, which has been developed for a wind turbine wake simulation subjected to a neutral atmospheric surface layer (van der Laan et al., 2015c). The model includes a scalar function f_{P}, which acts as a local turbulence lengthscale limiter in regions with highvelocity gradients, i.e., the near wake. In previous works (van der Laan and Sørensen, 2017b; van der Laan et al., 2017), the k–ε–f_{P} model was combined with an atmospheric boundary layer inflow model employing a constant pressure gradient, Coriolis forces, and the global turbulence lengthscale limiter of Apsley and Castro (1997), in order to simulate the effect of an idealized ABL on a wind farm. In addition, ambient turbulence sources were used to avoid zero turbulence quantities above the ABL. This ABL model does not require a temperature equation nor a buoyancy source term and can represent stable conditions by setting a low value of the maximum turbulence length scale ℓ_{max} by replacing C_{ε,1} with ${C}_{\mathit{\epsilon},\mathrm{1}}^{*}={C}_{\mathit{\epsilon},\mathrm{1}}+({C}_{\mathit{\epsilon},\mathrm{2}}{C}_{\mathit{\epsilon},\mathrm{1}})\mathrm{\ell}/{\mathrm{\ell}}_{\mathrm{max}}$, as introduced by Apsley and Castro (1997). Note that the turbulence length scale is a model parameter defined as $\mathrm{\ell}\equiv {C}_{\mathit{\mu}}^{\mathrm{3}/\mathrm{4}}{k}^{\mathrm{3}/\mathrm{2}}/\mathit{\epsilon}$. When the global turbulence lengthscale limiter of Apsley and Castro (1997) is applied to threedimensional flows containing turbulence length scales larger than the maximum values set by ℓ_{max}, one can observe nonphysical behavior in the solution and also numerical problems (van der Laan et al., 2017). This can occur for flow over complex terrain with large hills, for large wind farms, and for clusters of wind farms, as investigated in the present work. The nonphysical behavior typically leads to an artificially slow wake recovery. Therefore, we propose an alternative atmospheric inflow that sets the boundary layer height by an inversion layer using a prescribed potential temperature profile in combination with a nonzero ℬ instead of using the global turbulence lengthscale limiter of Apsley and Castro (1997). The new model does not employ an active temperature equation in order to maintain a steadystate model because the use of an active temperature equation in combination with a nonlinear temperature profile results in an unsteady inflow model, where the boundary layer height keeps growing with time.
We employ an analytically prescribed potential temperature profile θ by defining a zero and a constant temperature gradient at the ground and above the inversion height z_{i}, respectively, including a smooth transition in between the two regions:
where $\partial \mathit{\theta}/\partial z{}_{c}$ is the constant temperature gradient above z_{i}, and z_{T} is half the width of the transition zone around z_{i}. The latter is clear when taking a firstorder Taylor expansion of Eq. (7), $\partial \mathit{\theta}/\partial z\approx \mathrm{1}/\mathrm{2}[\mathrm{1}+(z{z}_{i})/{z}_{\mathrm{T}}]\partial \mathit{\theta}/\partial z{}_{c}$, and setting $z={z}_{i}\pm {z}_{\mathrm{T}}$. A temperature profile can be obtained by integrating over the height z and using $\mathit{\theta}(z=\mathrm{0})={\mathit{\theta}}_{\mathrm{0}}$ as the surface temperature:
There is a numerical issue with cosh (z) going to infinity for large z, which can be addressed by rewriting cosh (z) as $(\mathrm{1}+{e}^{\mathrm{2}x})/\left(\mathrm{2}{e}^{x}\right)$:
We set ${z}_{\mathrm{T}}/{z}_{i}$ as equal to 0.2; larger values would result in more smoothing (wider transition), and other values could be investigated in future work.
The effect of the prescribed temperature is represented by a buoyancy source term in the k and ε transport Eqs. (5) and (6):
with σ_{θ}=0.74 as the Prandtl number for temperature and g=9.81 m s^{−2} as the gravitational acceleration constant. We do not use a buoyancy source term in the momentum equation, and a constant air density of 1.225 kg m^{−3} is employed in the present work, as we only simulate wind farms in flat terrain. The buoyancy source term in the ε equation (Eq. 6) is multiplied by a constant C_{ε,3}:
which is based on the transient ABL model of Sogachev et al. (2012) without the global turbulence lengthscale limiter ℓ_{max} (or using ℓ_{max}=∞). In addition to the buoyancy sources, we also add ambient sources terms (Spalart and Rumsey, 2007) S_{k,amb} and S_{ε,amb} to the k and ε transport Eqs. (5) and (6), respectively, in order to prevent zero values of k and ε above the ABL, similar to van der Laan et al. (2020) but replacing ℓ_{max} with z_{i}:
with k_{amb}, ε_{amb}, ℓ_{amb}, and I_{amb} as the ambient values for k, ε, the turbulence length scale, and the turbulence intensity above the ABL, respectively. We set ${C}_{\mathrm{amb}}={\mathrm{10}}^{\mathrm{7}}$ and ${I}_{\mathrm{amb}}={\mathrm{10}}^{\mathrm{5}}$, which are lower values compared to previous work (van der Laan et al., 2020) but are necessary in order to get numerically stable results and are still low enough to not have an influence on the inflow profiles.
Since we lack measurements of the freestream, we need an alternative source to determine the input parameters for the inflow model. The New European Wind Atlas (NEWA) database (Dörenkämper et al., 2020) comprises 30 years (1989–2018) of simulated output of the WRF model for Europe. We extended the simulated period to cover 2019 and 2020 using the same model configuration of NEWA but only within NEWA's Great Britain subdomain. We extract data for a 3year time period, 1 January 2018–1 January 2021, at the grid point closest to the wind farm cluster center at a latitude and longitude of 53.21647 and 1.10947^{∘}. The grid point is only 1.12 km from the cluster center. The data are filtered for the same criteria as the measurements: a wind speed bin of 7–9 m s^{−1} (at z=100 m) and a Obukhov length range of $\leftL\right\ge \mathrm{200}$ m. Based on the WRF model output for these 3 years, the boundary layer height is 512 m, the wall temperature (skin temperature) is 285 K, and the turbulence intensity based on k (I_{ref}) at z=100 is 4.4 %. For the RANS ABL inflow model, we use the same wall temperature and I_{ref} at z=102 m, an inversion height of z_{i}=1000 m, and an inversion strength of 0.005 K m^{−1}. This results in a wind speed ABL height of about 700 m, as depicted in Fig. 5, where results of a precursor simulation are shown employing EllipSys1D (van der Laan and Sørensen, 2017a). Figure 5 also includes the selected NEWA data simulated by the WRF model (only data up to 500 m were extracted) and a neutral atmospheric surface layer (ASL) using the same I_{ref} and U_{ref}. The results from the WRF model and the prescribed temperature predict a similar wind speed, wind veer, and turbulence intensity around the rotor area, but the inflow profiles are quite different above the wind turbine. The potential temperature profile from the WRF model only matches near the ground and indicates stable conditions above 75 m or a very shallow inversion height. Setting a lower z_{i} in the prescribed temperature model in order to obtain an ABL height closer to that output by the WRF simulations is possible. However, numerical convergence problems can occur when such a shallow ABL inflow is applied for wind farms with large wind turbines, as discussed in more detail in Appendix B. The Coriolis parameter is set to $\mathrm{1.168}\times {\mathrm{10}}^{\mathrm{4}}$ s^{−1} based on the latitude of the wind farm cluster center. The geostrophic wind speed and roughness length are found by an optimization procedure employing the previously mentioned parameters and using a reference wind speed, U_{ref} of 8 m s^{−1} at z=102 m. The input and derived parameters are summarized in Table 2.
In the present work, we only use one inflow profile for all simulations, for simplicity and also when the inflow wind speed cases are not equal to 8 m s^{−1}. For example, the ADs are controlled using a lookup table based on force calibration simulations (van der Laan et al., 2015a), for which the entire wind speed range of the operational regime is computed for a single wind turbine by using different C_{T} values while keeping the inflow constant. This means that we assume Reynolds number similarity, which is not valid for the prescribed temperature inflow model because it follows a Rossby and Zilitinkevich number similarity, as shown in Appendix A. However, since the AD controller is only applied for wind farm flows around a wind speed of 8 m s^{−1}, we neglect the effect of small variations in these dimensionless numbers.
The new inflow model still includes a low eddy viscosity region above the ABL height, similar to the ABL model based on ℓ_{max}, as shown by the small turbulence length scale in Fig. 5d for z>7z_{ref}. If this region is included in the inflow for a threedimensional simulation of a large wind farm cluster (e.g., in order of 200×200 km^{2} or larger), then numerical instabilities can occur (van der Laan et al., 2017). Here, the largest horizontal domain is 169×208 km^{2}, which is just small enough to avoid numerical instabilities if the domain height is set to 10 D_{ref} (about 1 km). Furthermore, if the prescribed temperature model is employed with an inversion height that results in a shallow ABL, where the wind turbines are (partially) operating in the low eddy viscosity region, then numerical problems can occur as well. Appendix B illustrates the issue of employing shallow ABL inflows to wind farm cluster simulations, and the new prescribed temperature inflow model is also compared with the inflow model based on ℓ_{max}. We find that the use of Rayleigh damping in a steadystate method as RANS does not remove the numerical instabilities, possibly because Rayleigh damping methods have been developed for transient simulations (Durran and Klemp, 1983). The numerical issues associated with the low eddy viscosity region above the ABL need to be further investigated in future work.
3.2 WRF
In the previous section, output from a mesoscale long run using the WRF model was used to derive input parameters for the RANS inflow model. Another set of mesoscale simulations using the WRF model was performed on the Dudgeon wind farm area (including the Dudgeon, Sheringham Shoal, and Race Bank wind farms) but using an implementation of the WRF model version 3.7.1, in which the explicit wake parameterization (EWP) is included (Volker et al., 2015). The EWP represents mesoscale wind farm effects on the atmospheric flow. Turbineinduced drag forces are formulated as gridcellaveraged forces, while windturbineinduced turbulence is treated via an explicit subgrid scale turbulence diffusion formulation (Volker et al., 2015). In contrast to WRF's native wind farm parameterization (Fitch et al., 2012), no explicit source of turbulent kinetic energy (TKE) is added to the TKE equation by the EWP. The simulation uses a oneway nested domain setup with three domains (Fig. 6). The innermost domain has a grid spacing of 1 km with 216×216 grid points. More details about model configurations as well as initial and boundary conditions are provided in Appendix D. This simulation is aimed to be compared with PyWakeEllipSysbased modeling approaches of wind farm cluster wakes, which use idealized representations of the ABL and do not include mesoscale effects. To provide a fair comparison, the simulation period (calendar year of 2018) has been screened for weather episodes that fulfill the following conditions at the WRF gridcell center closest to the L02 turbine of Dudgeon (Fig. 1): (1) inflow wind speed at hub height of 7–9 m s^{−1}, (2) inflow direction of 232.5–237.5^{∘}, and (3) nearneutral stability conditions (Obukhov length L>200 m). To avoid the inclusion of isolated cases, all conditions are required to persist for at least 1 h. This resulted in a sample size of 31 temporal snapshots of 10 min (instantaneous wind speeds) from three distinct weather events in 2018. Furthermore, a reference simulation without wind farms covering the same calendar period is run and the same 31 temporal snapshots of 10 min are extracted for the normalization of the wind speed results from the simulation with wind farms. This is performed to correct for the impact of the background mesoscale flow, e.g., coastal effects or largescale wind speed gradients in the horizontal wake profile.
3.3 TurbOPark engineering wake model
The TurbOPark model from Nygaard et al. (2020) and Pedersen et al. (2022) is an analytical wake model optimized for simulating large offshore wind farms and farmtofarm interaction. It employs a Gaussian wake deficit profile, but differs in its definition of the wake expansion from other analytical models by using nonlinear streamwise wake expansion in contrast to the more commonly used linear expansion. The nonlinearity derives from assuming that the wake expansion rate is proportional to the local turbulence intensity, which is assumed to be a combination of atmospheric and turbineadded turbulence. In combination with Frandsen's model for the streamwise attenuation of turbulence intensity (Frandsen, 2007), the wake expands quickly near the rotor but expands progressively slower moving further downstream. This conserves wakes deficits over large distances, making the model suitable for assessing wake losses in large arrays or between farms.
The TurbOPark wake model is implemented in DTU Wind's opensource wind farm simulation tool PyWake (DTU Wind Energy, 2022a). PyWake is fully modular, so it gives complete freedom how wake, blockage, and other submodels are combined to define the engineering wind farm model for the annual energy production (AEP) simulation. The local averaged disk velocity is used to determine thrust, power, and deficits, whereas the freestream streamwise turbulence intensity of 5.5 %, estimated from the turbulence intensity based on the TKE as $\mathrm{4.4}/\mathrm{0.8}$ (van der Laan et al., 2015c), is used at all turbines to determine wake expansion, and the superposition of wakeadded turbulence is not taken into account. A windfarmoptimized version of the selfsimilar single turbine blockage model (Troldborg and Meyer Forsting, 2017), derived from RANS simulations of multiple fullscale turbines, accounts for wind farm blockage effects, including speedups. A mirror plane models the ground effect for both wakes and blockage. Wake deficits are superposed by the root sum square, whereas blockage effects are linearly summed, as they constitute small perturbations. The total wind farm flow field is obtained by summing blockage and wake contributions. The flow solution is obtained iteratively, as deficit information needs to be passed up and downstream. The computation is accelerated by initializing the solution with a wakesonly simulation. Around cutin and cutoff wind speed, the discontinuous nature of the thrust curve may spoil convergence, as certain turbines switch on and off between iterations. This is avoided by identifying unstable turbines and switching them off permanently.
TurbOPark, as implemented in PyWake, is used with two different setups: one reflects the original model setup from Nygaard et al. (2020) and Pedersen et al. (2022), and the other is a revised setup, where the ground model for the wake deficit is switched off, and a wake expansion coefficient of 0.06 (instead of 0.04) is employed, as we find that the original model setup overpredicts wind farm wake effects. The revised setup is one way of reducing these effects such that the results compare better with the simulated results from RANS and the WRF model for the investigated double wind farm case, as shown in Sect. 4.2.1. However, other setups, with linear wake summation for instance, could give as convincing predictions if tuned. An overview of the PyWake setups are given in Appendix C.
4.1 Model verification
4.1.1 Horizontal drag distribution: binning vs. Gaussian methods
The horizontal drag distribution in the AWF model represents the wind turbine density per cell. If large horizontal cells are used in the AWF model, e.g., in the order of the wind turbine interspacing in a wind farm, then the AWF grid orientation with respect to the wind farm layout can have a large influence on the wind farm wake shape when counting or binning the number wind turbines per cell. We refer to this method as the binning method. One can solve the issue by using a superposition of twodimensional Gaussian functions for the wind turbines to represent the number of wind turbines per cell. To illustrate the difference between the binning and the Gaussian methods, a square wind farm with 8×8 wind turbines using a spacing of 8 D is simulated for a wind direction aligned with the diagonal. In practice, we simply rotate a square wind farm layout by 45^{∘}. The wind farm is simulated using RANSAD (as a reference) and five RANSAWF models, with a large horizontal grid spacing in the drag distribution grid (Δ=8 D), one using the binning method and four using the Gaussian method with smoothing by employing different standard deviations: $\mathit{\sigma}={\mathit{\sigma}}_{x}={\mathit{\sigma}}_{y}=\mathit{\{}\mathrm{\Delta}/\mathrm{4},\mathrm{\Delta}/\mathrm{2},\mathrm{\Delta},\mathrm{2}\mathrm{\Delta}\mathit{\}}$. The flow domain for the AWF models is the same as the one used for the AD model, namely a horizontal spacing of $D/\mathrm{8}$ at the wind farm location and a maximum of 1 D in the wind farm wake up to four wind farm lengths L_{WF} downstream. Figure 7 shows results of these simulations in terms of streamwise velocity normalized by the freestream (Fig. 7f–k) and added wake turbulence intensity (Fig. 7l–q); both quantities are extracted at hub height. The vertically integrated drag distribution of the AWF models (after interpolation to the CFD grid) is shown in Fig. 7a–e. When the binning method is applied, the drag distribution becomes inhomogeneous due to aliasing effects, see Fig. 7a. This results in artificial wind farm wake shapes, as depicted in Fig. 7g and m, compared to the AD simulation shown in Fig. 7f and l. The same problem occurs with the Gaussian method when not enough smoothing is applied, e.g., the case for $\mathit{\sigma}=\mathrm{\Delta}/\mathrm{4}$ (Fig. 7b, h, and n). However, when more smoothing is used, i.e., $\mathit{\sigma}=\mathrm{\Delta}/\mathrm{2}$ and σ=Δ, the horizontal drag force distribution becomes uniform, as shown in Fig. 7c–e (as expected for a regular layout), and the artificial wake shape disappears (Fig. 7i–k and o–q).
The results shown in Fig. 7 are also depicted in Fig. 8 but in terms of lateral profiles of streamwise velocity at hub height and at three different downstream locations. It is clear that by using the binning and the Gaussian method without enough smoothing $\mathit{\sigma}=\mathrm{\Delta}/\mathrm{4}$, artificial wake shapes appear, which are not present in the simulations using ADs. The results using the largest smoothing σ=2Δ reduce the accuracy with respect to the AD simulation because the effective resolution also reduces with an increased level of smoothing. However, using σ=Δ can lead to a checkerboard horizontal drag force distribution for a particular setup, e.g., for a nonrotated square wind farm layout with 8×8 wind turbines, 8 D spacing, and Δ=2 D. This problem disappears when using σ=2Δ. From a numerical point of view, smoothing body forces with a Gaussian in a numerical simulation should use a minimal value of σ=2Δ in order to prevent wiggles in the solution. This is based on the Gaussian smoothing used to represent the blade forces in actuator line models (Troldborg, 2008; Forsting and Troldborg, 2020). Therefore, we apply a smoothing of σ=2Δ for all other simulations using the AWF model. If one would like to increase the accuracy of the AWF model, then it is recommended to refine Δ with σ=2Δ and not solely reduce σ.
4.1.2 Horizontal drag distribution: AWF grid and CFD grid spacing
The effect of the horizontal resolution in the AWF grid Δ is depicted in Fig. 9 using a wind farm with 8×8 wind turbines, a turbine spacing of 8 D, and a rowaligned wind direction. Results of six AWF models with different Δ values are shown, and the results are compared with results from the AD model (using the same CFD grid spacing as the AD model $\mathit{\delta}=D/\mathrm{8}$) in terms of streamwise velocity normalized by the freestream (Fig. 9g–m) and added wake turbulence intensity (Fig. 9n–t), both extracted at hub height. The vertically integrated drag distribution of the AWF models are shown in Fig. 9a–f. When refining the horizontal spacing in the AWF model, the solution approaches the one of the AD model. However, there are still differences between the AWF model with $\mathrm{\Delta}/\mathrm{8}$ and the AD model mainly at the wind farm location, as shown in Fig. 9g, h, n, and o. If a closer match is desired, one would need to further refine Δ. The need for a finer AWF grid spacing in order to approach an AD simulation result is not a surprise because the drag distribution in the AWF model is represented by a Cartesian grid, and thus more cells than eight are required to represent a rotor disk area by square cells. This is also why our AD model uses a polar grid for each wind turbine in the RANSAD setup. When Δ≥D or $\mathrm{\Delta}\ge s/\mathrm{8}$, with s as the wind turbine interspacing, the individual wind turbines are no longer visible in the horizontal drag distribution of the AWF model, as shown in Fig. 9.
The difference between the AD and the AWF simulation results are depicted in Fig. 10 in terms of horizontal velocity, Fig. 10a and b, and turbulence intensity, Fig. 10c and d. A range of wind directions are simulated between 270–315^{∘} with an interval of 5^{∘}. For each wind direction, the mean and maximum absolute difference between the AD and AWF models is calculated from cross planes centered at $(y,z)=(\mathrm{0},{z}_{\mathrm{H}})$, with a width L_{WF} and a height D, for a range of downstream locations to estimate the impact of Δ on a potential downstream wind farm. Subsequently, the results are either averaged over the wind directions, Fig. 10a and c, or filtered for the maximum absolute difference, Fig. 10b and d. In addition, results are shown using the same CFD grid as the AD model $\mathit{\delta}=D/\mathrm{8}$, and when the AWF model spacing is the same as the CFD grid spacing Δ=δ. The first provides a more fair comparison, while the latter is used in practice when applying the AWF model in a wind farm cluster. As expected, the difference between the AWF and AD models increases with coarsening the horizontal AWF model resolution Δ. However, the difference is not much influenced by coarsening the CFD grid by using Δ=δ; the main effects are observed for the large AWF model spacing of Δ=4 D and 8 D. The errors are the largest at the wind farm area, which is not our area of interest when upstream wind farms are modeled as AWFs. One could select a value for Δ depending on the desired accuracy. For the present work, we select Δ=2 D, which results in a mean error in streamwise velocity and turbulence intensity of ±0.4 % and 0.3 %, respectively, for the wind farm wake at x>2L_{WF}. The maximum absolute errors at x>2L_{WF} in streamwise velocity and turbulence intensity are in the order of 6 % and 1 %, respectively, but often only occur locally. A similar grid study has been performed for the same wind farm layout but with a wind turbine spacing of 4 D; the results are depicted in Fig. 11. In general, the differences between the AWF and AD models are larger compared to the lowerdensity wind farm. However, the mean error in streamwise velocity for the chosen setup using Δ=2 D and Δ=δ is still only ±0.7 % for x>2 L.
4.2 Model validation
4.2.1 Double wind farm case: Sheringham Shoal and Dudgeon
The effect of Sheringham Shoal on Dudgeon for a wind direction of 235±2.5^{∘} and a wind speed of 8±0.5 m s^{−1} is simulated with RANS, TurbOPark, and the WRF model. Three RANS setups are used: two employ either all ADs or all AWF models for both wind farms, while the third uses ADs for the downstream wind farm (Dudgeon) and an AWF model for the upstream wind farm (Sheringham Shoal). The horizontal canopy spacing of the AWF models is set to 2 D. The RANS models are simulated for a wind speed of 8 m s^{−1} and every 5^{∘} between 220–250^{∘} and post processed by a Gaussian filter of 5^{∘}. PyWake's TurbOPark implementation follows the original formulation by Nygaard et al. (2020) (as close as possible) and a revised setup where the ground model is switched off, and a wake expansion coefficient of 0.06 is employed instead of 0.04. TurbOPark is used for the same wind speeds as the RANS setups and for every 1^{∘}and is also post processed using the same Gaussian filter and a subsequent linear average between 235±2.5^{∘}. The Gaussian filter is applied to represent the uncertainty of the measured wind direction (Gaumond et al., 2014; van der Laan et al., 2015b); a mathematical description can be found in Antonini et al. (2019). The results of the double wind farm case are depicted in Fig. 12; in terms of horizontal wind speed contours Fig. 12a–e; wake magnitude, Fig. 12f; and shape, Fig. 12g; the latter two are obtained from the wind turbine power of the front row of Dudgeon using the wind turbine power curve. The wake shape represents the wind farm wake velocity deficit normalized by its mean value, while the wake magnitude refers to the wind farm wake velocity deficit normalized by the freestream. A wake shape is used for validation because the SCADA measurements of the front row turbines are both used to determine the freestream wind speed U_{ref} and to measure the upstream wind farm wake due to the lack of freestream measurements, as discussed in Sect. 2. Figure 12h also shows the results of the wake of Sheringham Shoal extracted along a 145–325^{∘} transect (perpendicular to the main wind direction of 235^{∘}), at 1.4 km upstream of Dudgeon. These results are used in order to compare results of all RANS models and TurbOPark with those using the WRF model because it is not possible to use the wind turbine power output of the front row of wind turbines from the WRF simulation to extract the wake effect of Sheringham Shoal due to the large horizontal resolution. The upstream distance of 1.4 km is chosen to avoid cells that include a wind turbine in the WRF model setup. In addition, U_{∞} for the WRF model represents the wind speed from a simulation without turbines, while U_{∞} is the actual inflow wind speed in the RANS and TurbOPark models. Finally, results from simulating the Dudgeon wind farm alone using RANSAD are also depicted in Fig. 12e, which show that wind farm induction is in the order of 1 % of the freestream and much smaller than the effect of Sheringham Shoal on Dudgeon predicted by the RANS and TurbOPark simulations that also include Sheringham Shoal.
The contour plots in Fig. 12a–e are rotated to better visualize the incoming flow for the front row wind turbines and transect, for which the results are depicted in Figs.12f–h. Figure 12b–d shows comparable wind farm wakes between the RANS models. The horizontal wind speed flow fields from the WRF model results (Fig. 12e) are quite different compared to those obtained from the RANS models. For example, the WRF model simulation lacks regions of increased wind speed at the sides of Sheringham Shoal, although the flow in between the wind farms is more similar. The horizontal wind speed deficits of TurbOPark (using the revised setup) shown in Fig. 12a have distinct single wind turbine wakes due to the applied superposition method because the depicted contours in Fig. 12a reflect a single wind direction of 235^{∘}.
There is a maximum of 0.9 % difference between the RANSAD and RANSADAWF models in terms of the equivalent wind speed extracted from the first row of wind turbines (with respect to the freestream wind speed). Note that it is not possible to perform this exercise with the RANSAWF simulations, since both wind farms are represented by an AWF model. This difference in wind speed is smaller along the transect, and all three RANS models predict a wake magnitude in the range of the WRF model results, as shown in Fig. 12e. The results from the RANSADAWF and RANSAWF models are nearly identical (Fig. 12h) because the upstream effect on Dudgeon is small at the transect. TurbOPark in its original formulation predicts stronger wind farm wake effects compared to RANS for the front row of Dudgeon (Fig. 12f) and compared to the WRF model results for the transect (Fig. 12h). However, in terms of wake shape (Fig. 12g), all models compare relatively well with the SCADA measurements, which shows the difficulty of validating the models with SCADA without freestream measurements. TurbOPark in its revised setup predicts much better results compared to those of RANS and the WRF model, partially due to removing the mirror plane for the wake deficits. Whilst not relevant in the near wind farm wake, the mirror deficits become active at large distances from the turbine of origin when hitting Dudgeon.
The wind farm power loss of Dudgeon due the presence of Sheringham Shoal is depicted in Fig. 13, where results of the three RANS methods are shown as a function of wind direction in Fig. 13a. In addition, the RANSADAWF and RANSAWF models are also simulated with a finer horizontal resolution of Δ=1 D. Furthermore, the difference with respect to the RANSAD simulation is plotted in Fig. 13b. When the downstream wind farm is also modeled as an AWF (RANSAWF), then the difference with the RANSAD model results is generally larger (up to 1.5 %) except for the wind directions of 245 and 250^{∘}. The error in the RANSAWF model is not strongly related to the grid size, since refining the horizontal spacing by a factor of two only changes the results of the RANSAWF by a maximum of 0.1 %. We suspect that the main source of error is the fact that the AWF model is controlled as an entire object instead of individual wind turbine control, as performed for the RANSAD model. When the downstream wind farm is operating in a half wind farm wake situation, mostly applicable for the wind directions 220–225 and 245–250^{∘}, then the entire AWF is still affected if the change in the AWF volumeaveraged wind speed or wind direction changes the wind farm thrust coefficient, which would not be the case when the downstream wind farm is represented by ADs. If the latter is important, then one could consider splitting each wind farm into several AWF models or simply model the wind farm of interest with ADs following the RANSADAWF approach.
4.2.2 Wind farm cluster
The wind farm cluster consisting of the three wind farms Dudgeon, Sheringham Shoal, and Race Bank is simulated with RANSADAWF, where Dudgeon is represented by ADs, and the other two wind farms are modeled with the AWF model. A range of wind directions between 190 and 300^{∘} is simulated for every 5^{∘} using a wind speed of 8 m s^{−1}. The wind turbine power of the first row of Dudgeon is used to calculate the effective wind speed, and the results are Gaussian averaged using the same standard deviation as employed in Sect. 4.2.1. The RANS simulations are compared with the SCADA measurements for the southern and western front rows of Dudgeon in Figs. 14 and 15, respectively. The results are normalized by U_{ref}, which represents the mean row wind speed that is used to determine the freestream in the SCADA measurements due to the lack of freestream measurements. Hence, the wake shape is only validated, while the wake magnitude is not, as also discussed in Sect. 4.2.1.
Figures 14 and 15 indicate that the RANS simulations capture the overall trend of the wind farm wake shape as a function of wind directions. For the southern row of Dudgeon, the measurements indicated a more pronounced wind farm wake shape, best visible by the difference in simulation and measurements at the corner wind turbines J05 and A05 for the wind directions between 205–245^{∘} and 235–250^{∘}, respectively. This could indicate stronger wake effects in the measurements, possibly due to the effect of near stable conditions, which is not accounted for in the simulations.
The RANS results of the western row in Fig. 15 are more difficult to interpret because the western row width is not large enough to capture the entire wind farm wake of Race Bank, which leads to flat wake shapes when the western front row of Dudgeon is operating in the full wind farm wake of Race Bank, as seen for the wind direction of 265^{∘}.
4.3 Computational effort
The main motivation for the AWF model is the reduction in computational effort compared to using ADs. Table 3 provides an overview of the CFD grid sizes and computational effort per flow case (i.e., a single inflow wind speed and direction) for running one, two, and three wind farms with different combinations of wind farms modeled as AD and/or AWFs. The underlying CFD solver in PyWakeEllipSys, EllipSys3D, is a flow solver based on a block structured grids, which uses a message passing interface (MPI) for communication. In the present work, we use block (edge) sizes of 32 and 64, corresponding to 32^{3} and 64^{3} cells per block, respectively. The largest part of the computational effort is used to solve for the pressure correction equation. A multigrid algorithm is applied to solve the pressure correction equation, and the coarsest level is made parallel on a single node using a sharedmemory MPI. The most recent version of EllipSys3D is expected to scale efficiently and be able to run with one block per processor for grids with 100 000 blocks of 4^{3}. Four shared memory CPUs are used for the wind farm cluster simulations using three wind farms and only ADs due to the large number of blocks (4550), while the other simulations are not using shared memory CPUs. The AD model in EllipSys3D was recently made more scalable by distributing the same number of ADs to each CPU as much as possible; initial tests have shown good scalability up to 1000 ADs (van der Laan et al., 2023). Further optimizing the computational effort and scalability of large wind cluster simulations using ADs is planned in the future.
The RANSAD grids for multiple wind farms could be reduced by using more complex grid topologies, where each wind farm is situated in a refined region instead of using a single large refined area that includes all wind farms. However, it is more challenging to generate such a grid topology compared to boxtype domains. Furthermore, it should be noted that the simulations are performed on an inhouse shared computer cluster. Hence, the listed CPU hours and wall clock time in Table 3 should be interpreted as indicative. The computer cluster consists of 516 computational nodes, each one equipped with two 16core AMD EPYC 7351 2.9 GHz processors and 128 GB RAM of memory (Technical University of Denmark, 2019).
When the investigated wind farm cluster is simulated with Dudgeon as ADs and the other wind farms as AWF models, 76 % of cells can be saved compared to modeling the wind farm cluster with only ADs due to a reduction in horizontal spacing outside the Dudgeon wind farm area, as listed in Table 3. For the double wind case, 63 % of the cells can be saved, and the number of CPU hours per flow case is reduced by 75 %. When the wind farm cluster consisting of three wind farms is solely modeled by AWF models, then the grid size is reduced by 99.5 %, and a flow case can be solved in about 1.5 min. using 200 CPUs, which is a reduction of 99.9 % in terms of CPU hours. However, each AWF model requires information of a wind farm thrust coefficient C_{T,wf} and the wind farm power coefficient (if the AWF model is used to predict wind farm power) C_{P,wf} as a function of the AWF volumeaveraged wind speed and wind direction. If this input is derived from a set of RANS simulations using ADs, as performed in the present work, then the computational effort will be dominated by this step. For example, Table 3 shows that the Sheringham Shoal wind farm modeled by ADs takes about 0.20 h per flow case or 88 CPU hours using a grid of 115 million cells. The wind farm cluster modeled by three AWFs requires about 1206 CPU hours per flow case when including the wind farm precursor steps, which is still a reduction of 83 % compared to simulating all wind farms as ADs (6958 CPU hours per flow case), as shown by the last column of Table 3. One could attempt to use a simplified wake model to determine C_{T,wf} and C_{P,wf}, but one would risk a loss of accuracy. An initial study shows that is possible to predict C_{T,wf} and C_{P,wf} within a mean absolute error of around 1 % when using a RANSbased surrogate wind farm model (van der Laan et al., 2023).
A RANSbased wind farm parameterization, the AWF model, is proposed. It uses a wind farm thrust force as a momentum sink similar to a forest canopy model. The AWF model can be used as an obstacle model to a downstream wind farm of interest represented by ADs, or it can be used to estimate wind farm power when the downstream wind farm is also modeled as an AWF. The AWF model simulates wind farm interaction by a calibrated controller of wind farm thrust and power as a function of the wind farm volumeaveraged wind direction and wind speed.
When the horizontal spacing in the AWF model is refined, the wind farm flow approaches the results from a RANSAD wind farm simulation. This is achieved by calibrating the thrust force magnitude with precursor RANSAD wind farm simulations and employing a horizontal thrust force distribution in the form of wind turbine density using a superposition of the wind turbine coordinates represented by twodimensional Gaussian functions. The verification study showed that the Gaussian superposition method solves the problem of artificial wind farm wake effects that can occur when the number of wind turbines are binned for large horizontal cells, as current wind farm parameterizations implemented in numerical weather models, such as the WRF model, do.
A new atmospheric inflow model is introduced that is potentially more suited for wind farm cluster simulations because it does not rely on an ABL height set by a global turbulence lengthscale limiter that can result in nonphysical wind farm wakes. The proposed inflow model relies on a prescribed analytical temperature profile including an inversion height and inversion strength, while a temperature equation is not solved for in order to maintain a steadystate inflow model. The model is shown to be dependent on three nondimensional numbers. The new inflow model does not (and is not expected to) solve the problem of numerical instabilities related to the low eddy viscosity region above the ABL in combination with large horizontal domains associated with wind farm clusters. The problem is mitigated in the present work by using a relatively low domain height, which can introduce additional numerical blockage, thereby negatively influencing the predictions. An alternative solution needs to be found in the future to allow the more preferred, taller domain heights. In addition, the new ABL model requires more validation with measurements.
The proposed RANSAWF and inflow models are employed to simulate two neighboring wind farms and a cluster consisting of three wind farms, where either one of the wind farms is modeled by ADs, and the remaining wind farms are represented by AWF models (RANSADAWF) or all wind farm are AWF models (RANSAWF). The results for the double wind farm case are compared with TurbOPark, WRF, and RANSAD simulations (for the latter, all wind turbines are modeled by ADs) using the wind speed derived from the front row turbines of the downstream wind farm and the horizontal wind speed extracted from a transect 1.4 km upstream of the wind farm farthest downstream. The latter is performed because the chosen resolution for the WRF model simulations was not sufficient to resolve the front row wind turbine wind speed. While the overall horizontal wind speed at the wind farms in the WRF model simulations is quite different with respect to the RANS results, the horizontal wind speed at the transect from the WRF results compares well with those from the RANSAD, RANSADAWF, and RANSAWF models. In addition, the front row wind speed in the RANSADAWF only deviated by 0.9 % compared to the RANSAD results with respect to the freestream. The original formulation of TurbOPark shows stronger wind farm wake effects compared to the other models, but its result in terms of wind farm wake shape compared well with all models. This indicates that a comparison in terms of wind farm wake shape should not be the only type of validation. A revised TurbOPark setup, where the ground model is switched off and a larger wake expansion coefficient of 0.06 is used, predicts much better results compared to the higherfidelity models. Unfortunately, the RANSADAWF simulations of the wind farm cluster could only be validated with the shape of the front row wind speed because the SCADA measurements of this row are used to both measure the wake of the upstream wind farm and determine the freestream conditions due to a lack of concurrent freestream measurements. The trends of the RANSADAWF simulation results of the upstream wind farm wake shapes compare reasonably well with the results from the SCADA, although the measured shapes indicate stronger wind farm wake effects, possibly due to nearstable conditions that were not filtered out from the SCADA in order to maintain a large enough data set. More validation of both the AWF model and prescribed temperature inflow model is required. In addition, we need SCADA with concurrent inflow measurements in order to validate the magnitude of wind farm wakes and its impact on neighboring wind farms.
With the AWF model, one can simulate large wind farm clusters with RANS; the wind farm cluster validation case showed a reduction of 92.3 % and 99.9 % in CPU hours when two of three wind farms or all wind farms are represented by AWF models instead of using ADs, respectively, when the input wind farm thrust and power coefficients are known. If the wind farm power and thrust coefficients are calculated from a RANSAD simulation of each wind farm, as performed in the present work, then the reduction in computational effort is 82.7 % when the wind farm cluster validation case is solely modeled by AWF models. Simpler and faster models that can generate the AWF input are currently investigated in a followup work (van der Laan et al., 2023). This also enables a larger wind farm cluster to be efficiently simulated with the RANSAWF methodology in future work. Such a study would also benefit from a comparison with WRF model results to further investigate the effect of neglecting the mesoscales in the RANSAWF approach.
The AWF model can be used to further develop wind farm parameterizations in the WRF; this can be achieved by implementing the twodimensional Gaussian superposition method of the wind turbine density and by using a wind farm drag coefficient that is dependent on the wind direction instead of using the wind turbine thrust curves. Idealized WRF simulations could be used to compare the resulting wind farm wakes with RANSAWF simulation results.
The current implementation of the AWF model does not include any sources in the turbulence transport equations because the employed horizontal spacing of 2 D turned out to be sufficiently fine. If larger horizontal cell sizes are desired, one could investigate the use of additional turbulencerelated source terms.
In a previous work (van der Laan et al., 2020), it was shown that the ABL model of Apsley and Castro (1997), including Coriolis forces, follows a Rossby similarity; here, four dimensional input parameters (namely the geostrophic wind speed G, roughness length z_{0}, Coriolis parameter f_{c}, and ℓ_{max}) can be reduced to two Rossby numbers, $R{o}_{\mathrm{0}}\equiv G/\left(\right{f}_{\mathrm{c}}\left{z}_{\mathrm{0}}\right)$ and $R{o}_{\mathrm{\ell}}\equiv G/\left(\right{f}_{\mathrm{c}}\left{\mathrm{\ell}}_{\mathrm{max}}\right)$. We typically generate a library of normalized ABL profiles for all possible solutions, which can then be used to find values for G and ℓ_{max} in order to get a desired wind speed and turbulence intensity at a reference height for a fixed roughness length and Coriolis parameter. The new prescribed temperature ABL model from Sect. 3.1.3 depends on five dimensional parameters – G, z_{0}, f_{c}, z_{i}, and $\partial \mathit{\theta}/\partial z{}_{c}$ – however, the normalized profiles only depend on three nondimensional numbers:
Here, Ro_{0} is the surface Rossby number, $R{o}_{{z}_{i}}$ is a Rossby number based on the inversion height, and N_{f} is a dimensionless number based on the ratio of the Brunt–Väisälä frequency and the Coriolis parameter, which is also referred to as the Zilitinkevich number (Esau, 2004; Kelly et al., 2019; Liu et al., 2021). The similarity is numerically proven by the collapse of simulation results for different values of G and f_{c}, as depicted in Fig. A1. The numerical setup is similar to van der Laan et al. (2020), but we have chosen a smaller time step set to $\mathrm{1}/\mathrm{N}$ instead of $\mathrm{1}/\left{f}_{\mathrm{c}}\right$ in order to maintain numerical convergence. From the figure, one can note that N_{f} mostly changes the effective sharpness of the inversion, as seen in ϕ(z), as well as the ABL turning angle ϕ(z=z_{0}); the latter effect is more pronounced for larger ${z}_{i}/{z}_{\mathrm{0}}$. It is also possible to use an alternative dimensionless number for quantifying the inversion strength, using an ABLwide bulk Richardson number (instead of N_{f}):
A new atmospheric inflow model is proposed in Sect. 3.1.3 and used to perform the RANS wind farm cluster simulations in Sect. 4. The new inflow model employs a prescribed temperature profile with an inversion that sets the ABL height. Previous work employed an ABL inflow model based on a global turbulence lengthscale limiter, ℓ_{max} (Apsley and Castro, 1997; van der Laan and Sørensen, 2017b; van der Laan et al., 2017). For a low value of ℓ_{max}, a lower ABL height is obtained. A major issue with the global turbulence lengthscale limiter is that it can also limit the wind farm wake turbulence length scales that could lead to a nonphysical slow wake recovery. This issue was mitigated by switching off the global lengthscale limiter in regions where highvelocity gradients are present (van der Laan and Sørensen, 2017b), although it is unclear if this ad hoc solution is sufficient for wind farm cluster simulations. Another challenge is that for low ABL heights, one can obtain numerical oscillations, which both the new model and the model based on ℓ_{max} can suffer from. This section is meant to illustrate these issues and to show the difference between the inflow models when applied to the wind farm cluster validation case.
Table B1 lists the input and derived parameters of a shallow and a tall ABL inflow using the new prescribed temperature ABL model and the ABL model based on ℓ_{max}. The tall ABL case is the same as used through out the article (Sects. 3.1.3 and 4), while the shallow ABL is an additional case meant to illustrate the challenges with the two inflow models. For the ℓ_{max} ABL inflow model, the derived values for G and ℓ_{max} are found by interpolating a precalculated library of ABL profiles (van der Laan et al., 2020) for given z_{0}, f_{c}, I_{ref}, and U_{ref} values. The prescribed temperature model uses an optimizer for G and z_{0} instead, as discussed in Sect. 3.1.3. The tall ABL profile from the ℓ_{max} ABL inflow model uses a similar roughness length (${z}_{\mathrm{0}}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{5}}$ m) as was calculated for the inflow model based on the prescribed temperature (Table 2), which results in a relative large value of ℓ_{max}, namely 53.3 m. The shallow ABL profile of the prescribed temperature ABL model is generated with a lower z_{i} (300 m), resulting in a larger derived z_{0} (Table B1). In addition, the value of ${z}_{T}/{z}_{i}$ as used in the prescribed temperature ABL model is increased for the shallow ABL to avoid an inflection of the wind speed profile below the super geostrophic jet. To lower the ABL height for a model based on ℓ_{max}, one can lower z_{0} or I_{ref}; we choose to set ${z}_{\mathrm{0}}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ m. This results in a derived value of ℓ_{max} of only 7.62 m. Finally, the ABL model based on ℓ_{max} is employed with a lower ambient turbulence intensity value I_{amb} compared to the prescribed temperature model because the ℓ_{max} ABL model is more sensitive to this parameter. The chosen values for each ABL model are sufficient to not influence the numerical solution (van der Laan et al., 2020). The results of all profiles are depicted in Fig. B1. Here it is clear that over the rotor swept area the profiles of wind speed, direction, and normalized turbulence intensity resulting from the prescribed temperature model are similar to those resulting from the ℓ_{max} model for the tall ABL case. The main difference can be found in the turbulence model length scale (Fig. B1d), where the prescribed temperature ABL model predicts larger length scales over the rotor swept area compared to the ℓ_{max} based ABL model. For taller altitudes, the models produce quite different ABL solutions, probably because the prescribed temperature model sets an explicit inversion strength, while the ABL model based on ℓ_{max} calculates this implicitly. For the shallow ABL case, both models predict very different profiles. This is because the model based on ℓ_{max} represents stable conditions for low values of ℓ_{max}, while the prescribed temperature ABL model reflects a conventionally neutral ABL.
The four inflow profiles listed in Table B1 are applied to the wind farm cluster validation case using the RANSADAWF model for a wind direction and wind speed of 235^{∘} and 10 m s^{−1}, respectively. The results of these simulations in terms of horizontal wind speed contours at hub height (z=102 m) are depicted in Fig. B2. Here, the prescribed temperature model for the tall ABL case (Fig. B2a) is same as shown earlier in Fig. 14g, and its results are similar to the results of the inflow model based on ℓ_{max} (Fig. B2b) for the tall ABL case. This is because a large ABL of about 1 km is set by using either a large z_{i} or a large value of ℓ_{max}. Figure B3 also shows a similar wake magnitude calculated by the two inflow models, where the difference between the models is in the order of 1 %. The wake magnitude is based on the front row wind turbine power of Dudgeon, as performed in Sect. 4.2. Hence, the results presented in Sect. 4 would not have been significantly different if the inflow model based on ℓ_{max} was used instead of the new inflow model. The shallow ABL inflow profiles applied to the wind farm cluster do not lead to a converged result using the domain height L_{z} set to 10 D_{ref}, as employed for the tall ABL case. This is because the low eddy viscosity region above the ABL causes numerical instabilities when it is included in the wind farm cluster domain. The issue can be mitigated by lowering the domain height to exclude the low eddy viscosity region, as performed in previous work (van der Laan et al., 2017). Here, we apply two different values, ${L}_{z}/{D}_{\mathrm{ref}}=\mathrm{2.5}$ D and ${L}_{z}/{D}_{\mathrm{ref}}=\mathrm{3}$ D, and the results are shown in Figs. B2c–B2e. It should be noted that these low domain heights are not desired because the increased numerical blockage causes artificial flow accelerations and reduced wake losses. When using ${L}_{z}/{D}_{\mathrm{ref}}=\mathrm{2.5}$ D the shallow ABL simulations converge for both models (Fig. B2c and d). However, for ${L}_{z}/{D}_{\mathrm{ref}}=\mathrm{3}$ D (Fig. B2e, prescribed temperature ABL model) and ${L}_{z}/{D}_{\mathrm{ref}}=\mathrm{4}$ D (Fig. B2f, ℓ_{max} ABL model), the wind farm cluster simulation starts to produce numerical instabilities towards the outlet. It is clear that both ABL models cannot be used to simulate a wind farm cluster subjected to a shallow ABL inflow without lowering the domain height to an undesired small value. In addition, it is not fully determined if the new inflow model actually performs better than the model based on ℓ_{max}, and more work is needed to improve the numerical behavior when considering low ABL heights. However, the prescribed temperature model offers a more physical method of setting an ABL height compared to using a global turbulence lengthscale limiter, and the new model also provides the possibility to explicitly set the inversion strength.
An overview of the numerical setup of the original and revised TurbOPark model, as implemented in PyWake (DTU Wind Energy, 2022a), is provided in Table C1.
An overview of the numerical setup of the WRF simulations, including wind farms, is provided in Table D1.
(Hersbach et al., 2020)(30 arcsec, Danielson and Gesch, 2011)(Poulter et al., 2015)(Donlon et al., 2012)(Hong et al., 2004)(Iacono et al., 2008)(Grell and Freitas, 2014)(Tewari et al., 2004)(Nakanishi and Niino, 2006)The numerical results are generated with proprietary software, although the data presented can be made available by contacting the corresponding author.
MPVDL drafted the article and produced the figures. MPVDL and OGS proposed the idea of a wind farm model with a drag coefficient as a function of wind direction. MPVDL proposed representing the wind turbine positions in the AWF model as twodimensional Gaussian functions. AMF suggested using a standard deviation that is twice the horizontal AWF grid size and proposed the use of Eq. (3) for the volumeaveraged velocity vector of an AWF model; both ideas are based on lessons learned from actuator line modeling. MPVDL proposed the inflow model using a prescribed temperature, and MK proposed the analytical form of the prescribed temperature profile. CDB and KSS prepared the SCADA measurements of Dudgeon. KSS proposed revising the setup of the TurbOPark model simulations. MI performed and postprocessed the WRF simulations. AMF performed preliminary RANSAD simulations of the double wind farm case including Dudgeon and Sheringham Shoal, and AMF also implemented the TurbOPark model in PyWake. AP post processed the NEWA data. NNS implemented the necessary improvements to make the large wind farm cluster RANS simulations possible with EllipSys3D in an efficient and scalable way. PER facilitated and supervised the research performed in this article. All authors contributed to the methodology and finalization of the paper.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The first author was inspired by Andrey Sogachev to model a wind farm as a forest canopy. Furthermore, the authors would like to thank Mads Mølgaard Pedersen for his support with the implementation of TurbOPark in PyWake. We also gratefully acknowledge the computational and data resources provided on the Sophia HPC Cluster at the Technical University of Denmark, https://doi.org/10.57940/FAFC6M81.
This is work has been cofinanced by Equinor ASA. Oscar GarcíaSantiago was supported by the European Union's Horizon 2020 research and innovation program under grant agreement no. 861291 as part of the Train2Wind (https://www.train2wind.eu/, last access: 23 May 2023) Marie SkłodowskaCurie Innovative Training Networks. Alfredo Peña's contribution is partly funded by Independent Research Fund Denmark through the “Multiscale Atmospheric Modeling Above the Seas” (MAMAS) project (no. 021700055B).
This paper was edited by Cristina Archer and reviewed by Gonzalo Pablo Navarro Diaz and two anonymous referees.
Abkar, M. and PortéAgel, F.: A new windfarm parameterization for largescale atmospheric models, J. Renew. Sustain. Energ., 7, 013121, https://doi.org/10.1063/1.4907600, 2015. a, b, c, d, e, f
Allaerts, D. and Meyers, J.: Gravity Waves and WindFarm Efficiency in Neutral and Stable Conditions, Bound.Lay. Meteorol., 166, 269–299, https://doi.org/10.1007/s1054601703075, 2018. a
Antonini, E. G., Romero, D. A., and Amon, C. H.: Improving CFD wind farm simulations incorporating wind direction uncertainty, Renew. Energy, 133, 1011–1023, https://doi.org/10.1016/j.renene.2018.10.084, 2019. a
Apsley, D. D. and Castro, I. P.: A limitedlengthscale k–ε model for the neutral and stablystratified atmospheric boundary layer, Bound.Lay. Meteorol., 83, 75–98, https://doi.org/10.1023/A:1000252210512, 1997. a, b, c, d, e, f, g, h
Bak, C., Zahle, F., Bitsche, R., Kim, T., Yde, A., Henriksen, L. C., Natarajan, A., and Hansen, M. H.: Description of the DTU 10 MW Reference Wind Turbine, Tech. Rep. I0092, Technical University of Denmark, https://orbit.dtu.dk/files/55645274/The_DTU_10MW_Reference_Turbine_Christian_Bak.pdf (last access: 23 May 2023), 2013. a, b
Boudreault, L.E.: Reynoldsaveraged NavierStokes and LargeEddy Simulation Over and Inside Inhomogeneous Forests, PhD thesis, Wind Energy Department, Technical University of Denmark, https://orbit.dtu.dk/files/117912168/PHD_0042_leboudreault.pdf (last access: 23 May 2023), 2015. a, b
Boussinesq, M. J.: Théorie de l'écoulement tourbillonnant et tumultueux des liquides, GauthierVillars et fils, Paris, France, 1897. a
Churchfield, M. J., Schreck, S. J., Martinez, L. A., Meneveau, C., and Spalart, P. R.: An Advanced Actuator Line Method for Wind Energy Applications and Beyond, in: 35th Wind Energy Symposium, 9–13 January 2017, Grapevine, Texas, https://doi.org/10.2514/6.20171998, 2017. a
COWI: Vindressource, layouts og energiproduktion for Bornholm I + II, Nordsøen II + III og området vest for nordsøen II + III, Tech. Rep. A13299423, COWI A/S, Lyngby, Denmark, https://ens.dk/sites/ens.dk/files/Vindenergi/23_vindressource_layouts_og_energiproduktion.pdf (last access: 23 May 2023), 2020. a
Danielson, J. J. and Gesch, D. B.: Global multiresolution terrain elevation data 2010 (GMTED2010), USGS, https://doi.org/10.3133/ofr20111073, 2011. a
Dellwik, E., van der Laan, M., Angelou, N., Mann, J., and Sogachev, A.: Observed and modeled nearwake flow behind a solitary tree, Agr. Forest Meteorol., 265, 78–87, https://doi.org/10.1016/j.agrformet.2018.10.015, 2019. a, b, c, d
Donlon, C. J., Martin, M., Stark, J., RobertsJones, J., Fiedler, E., and Wimmer, W.: The Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA) system, Remote Sens. Environ., 116, 140–158, https://doi.org/10.1016/j.rse.2010.10.017, 2012. a
Dörenkämper, M., Olsen, B. T., Witha, B., Hahmann, A. N., Davis, N. N., Barcons, J., Ezber, Y., GarcíaBustamante, E., GonzálezRouco, J. F., Navarro, J., SastreMarugán, M., Sīle, T., Trei, W., Žagar, M., Badger, J., Gottschall, J., Sanz Rodrigo, J., and Mann, J.: The Making of the New European Wind Atlas – Part 2: Production and evaluation, Geosci. Model Dev., 13, 5079–5102, https://doi.org/10.5194/gmd1350792020, 2020. a
DTU Wind Energy: PyWake v2.5.0, https://topfarm.pages.windenergy.dtu.dk/PyWake/ (last access: 23 May 2023), 2022a. a, b
DTU Wind Energy: PyWakeEllipSys v3.2, https://topfarm.pages.windenergy.dtu.dk/cuttingedge/pywake/pywake_ellipsys/ (last access: 23 May 2023), 2022b. a
Durran, D. and Klemp, J.: A compressible model for the simulation of moist mountain waves, Am. Meteorol. Soc., 111, 2341–2361, 1983. a
Esau, I. N.: Parameterization of a surface drag coefficient in conventionally neutral planetary boundary layer, Ann. Geophys., 22, 3353–3362, https://doi.org/10.5194/angeo2233532004, 2004. a
Fischereit, J., Schaldemose Hansen, K., Larsén, X. G., van der Laan, M. P., Réthoré, P.E., and Murcia Leon, J. P.: Comparing and validating intrafarm and farmtofarm wakes across different mesoscale and highresolution wake models, Wind Energ. Sci., 7, 1069–1091, https://doi.org/10.5194/wes710692022, 2022. a
Fitch, A. C., Olson, J. B., Lundquist, J. K., Dudhia, J., Gupta, A. K., Michalakes, J., and Barstad, I.: Local and Mesoscale Impacts of Wind Farms as Parameterized in a Mesoscale NWP Model, Mon. Weather Rev., 140, 3017–3038, https://doi.org/10.1175/MWRD1100352.1, 2012. a, b, c, d
Forsting, A. M. and Troldborg, N.: Generalised grid requirements minimizing the actuator line angleofattack error, J. Phys.: Conf. Ser., 1618, 052001, https://doi.org/10.1088/17426596/1618/5/052001, 2020. a
4coffshore.com: Global offshore wind farm database, https://map.4coffshore.com/offshorewind/ (last access: 23 August 2022), 2022. a
Frandsen, S.: Turbulence and turbulencegenerated structural loading in wind turbine clusters, PhD thesis, risøR1188(EN), https://orbit.dtu.dk/files/12674798/ris_r_1188.pdf (last access: 23 May 2023), 2007. a
Gaumond, M., Réthoré, P.E., Ott, S., Peña, A., Bechmann, A., and Hansen, K. S.: Evaluation of the wind direction uncertainty and its impact on wake modeling at the Horns Rev offshore wind farm, Wind Energy, 17, 1169–1178, https://doi.org/10.1002/we.1625, 2014. a
Göçmen, T., van der Laan, M. P., Réthoré, P. E., Peña Diaz, A., Larsen, G. C., and Ott, S.: Wind turbine wake models developed at the technical university of Denmark: A review, Renew. Sustain. Energ. Rev., 60, 752–769, 2016. a
Grell, G. A. and Freitas, S. R.: A scale and aerosol aware stochastic convective parameterization for weather and air quality modeling, Atmos. Chem. Phys., 14, 5233–5250, https://doi.org/10.5194/acp1452332014, 2014. a
Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz‐Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J.: The ERA5 global reanalysis, Q. J. Roy. Meteorol. Soc., 146, 1999–2049, https://doi.org/10.1002/qj.3803, 2020. a
Hong, S.Y., Dudhia, J., and Chen, S.H.: A Revised Approach to Ice Microphysical Processes for the Bulk Parameterization of Clouds and Precipitation, Mon. Weather Rev., 132, 103–120, https://doi.org/10.1175/15200493(2004)132<0103:aratim>2.0.co;2, 2004. a
Iacono, M. J., Delamere, J. S., Mlawer, E. J., Shephard, M. W., Clough, S. A., and Collins, W. D.: Radiative forcing by longlived greenhouse gases: Calculations with the AER radiative transfer models, J. Geophys. Res., 113, D13103, https://doi.org/10.1029/2008jd009944, 2008. a
Jonkman, J., Butterfield, S., Musial, W., and Scott, G.: Definition of a 5MW Reference Wind Turbine for Offshore System Development, Tech. rep., National Renewable Energy Laboratory, https://www.nrel.gov/docs/fy09osti/38060.pdf (last access: 23 May 2023), 2009. a
Kelly, M., Cersosimo, R. A., and Berg, J.: A universal wind profile for the inversioncapped neutral atmospheric boundary layer, Q. J. Roy. Meteorol. Soc., 145, 982–992, https://doi.org/10.1002/qj.3472, 2019. a, b
Kosović, B., Munoz, P. J., Juliano, T. W., Martilli, A., Eghdami, M., Barros, A. P., and Haupt, S. E.: ThreeDimensional Planetary Boundary Layer Parameterization for HighResolution Mesoscale Simulations, J. Phys.: Conf. Ser., 1452, 012080, https://doi.org/10.1088/17426596/1452/1/012080, 2020. a
Launder, B. and Spalding, D.: The numerical computation of turbulent flows, Comput. Meth. Appl. Mech. Eng., 3, 269–289, https://doi.org/10.1016/00457825(74)900292, 1974. a
Liu, L., Gadde, S. N., and Stevens, R. J. A. M.: Universal Wind Profile for Conventionally Neutral Atmospheric Boundary Layers, Phys. Rev. Lett., 126, 104502, https://doi.org/10.1103/PhysRevLett.126.104502, 2021. a, b
Michelsen, J. A.: Basis3D  a platform for development of multiblock PDE solvers, Tech. Rep. AFM 9205, Technical University of Denmark, Lyngby, Denmark, https://orbit.dtu.dk/files/272917945/Michelsen_J_Basis3D.pdf (last access: 23 May 2023), 1992. a
Mikkelsen, R.: Actuator Disc Methods Applied to Wind Turbines, PhD thesis, Technical University of Denmark, https://orbit.dtu.dk/files/5452244/Robert.PDF (last access: 23 May 2023), 2003. a
Nakanishi, M. and Niino, H.: An improved MellorYamada Level3 model: Its numerical stability and application to a regional prediction of advection fog, Bound.Lay. Meteorol., 119, 397–407, https://doi.org/10.1007/s1054600590308, 2006. a
Nygaard, N. G., Steen, S. T., Poulsen, L., and Pedersen, J. G.: Modelling cluster wakes and wind farm blockage, J. Phys.: Conf. Ser. 1618, 062072, https://doi.org/10.1088/17426596/1618/6/062072, 2020. a, b, c, d, e
Ott, S. and Nielsen, M.: Developments of the offshore wind turbine wake model Fuga, Tech. Rep. I0046, Technical University of Denmark, https://orbit.dtu.dk/files/118472784/DTU_Wind_Energy_E_0046.pdf (last access: 23 May 2023), 2014. a
Pedersen, J. G., Svensson, E., Poulsen, L., and Nygaard, N. G.: Turbulence Optimized Park model with Gaussian wake profile, J. Phys. Conf. Ser., 2265, 022063, https://doi.org/10.1088/17426596/2265/2/022063, 2022. a, b, c
Peña, A., Mirocha, J. D., and van der Laan, M. P.: Evaluation of the Fitch windfarm wake parametrization with largeeddy simulations of wakes using the Weather Research and Forecasing model, Mon. Weather Rev., 150, 3051–3064, 2022. a
Poulter, B., MacBean, N., Hartley, A., Khlystova, I., Arino, O., Betts, R., Bontemps, S., Boettcher, M., Brockmann, C., Defourny, P., Hagemann, S., Herold, M., Kirches, G., Lamarche, C., Lederer, D., Ottlé, C., Peters, M., and Peylin, P.: Plant functional type classification for earth system models: results from the European Space Agency's Land Cover Climate Change Initiative, Geosci. Model Dev., 8, 2315–2328, https://doi.org/10.5194/gmd823152015, 2015. a
Réthoré, P.E.: Wind Turbine Wake in Atmospheric Turbulence, PhD thesis, Aalborg University, Risø DTU, Roskilde, Denmark, https://orbit.dtu.dk/files/4548747/risphd53.pdf (last access: 23 May 2023), 2009. a
Réthoré, P.E., van der Laan, M. P., Troldborg, N., Zahle, F., and Sørensen, N. N.: Verification and validation of an actuator disc model, Wind Energy, 17, 919–937, https://doi.org/10.1002/we.1607, 2014. a, b, c
Rybchuk, A., Juliano, T. W., Lundquist, J. K., Rosencrans, D., Bodini, N., and Optis, M.: The sensitivity of the Fitch wind farm parameterization to a threedimensional planetary boundary layer scheme, Wind Energ. Sci., 7, 2085–2098, https://doi.org/10.5194/wes720852022, 2022. a
Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. O., Liu, Z., Berner, J., Wang, W., Powers, J. G., Duda, M. G., Barker, D. M., and Huang, X.Y.: A Description of the Advanced Research WRF Version 4, Tech. rep., NCAR Tech. Note NCAR/TN556+STR, NCAR, https://doi.org/10.5065/1dfh6p97, 2019. a
Sogachev, A., Kelly, M., and Leclerc, M. Y.: Consistent TwoEquation Closure Modelling for Atmospheric Research: Buoyancy and Vegetation Implementations, Bound.Lay. Meteorol., 145, 307–327, https://doi.org/10.1007/s1054601297265, 2012. a, b, c
Sørensen, J. N. and Larsen, G. C.: A Minimalistic Prediction Model to Determine Energy Production and Costs of Offshore Wind Farms, Energies, 14, 2, https://doi.org/10.3390/en14020448, 2021. a
Sørensen, J. N. and Shen, W. Z.: Numerical Modeling of Wind Turbine Wakes, J. Fluids Eng., 124, 393–399, https://doi.org/10.1115/1.1471361, 2002. a
Sørensen, J. N., Nilsson, K., Ivanell, S., Asmuth, H., and Mikkelsen, R. F.: Analytical body forces in numerical actuator disc model of wind turbines, Renew. Energy, 147, 2259, https://doi.org/10.1016/j.renene.2019.09.134, 2020. a, b
Sørensen, N. N.: General purpose flow solver applied to flow over hills, PhD thesis, Risø National Laboratory, Roskilde, Denmark, https://orbit.dtu.dk/files/12280331/Ris_R_827.pdf (last access: 23 May 2023), 1994. a
Sørensen, N. N., Bechmann, A., Johansen, J., Myllerup, L., Botha, P., Vinther, S., and Nielsen, B. S.: Identification of severe wind conditions using a Reynolds Averaged Navier–Stokes solver, J. Phys.: Conf. Ser., 75, 1–13, https://doi.org/10.1088/17426596/75/1/012053, 2007. a
Spalart, P. and Rumsey, C.: Effective inflow conditions for turbulence models in aerodynamic calculations, AIAA J., 45, 2544–2553, 2007. a
Storey, R. C., Norris, S. E., and Cater, J. E.: An actuator sector method for efficient transient wind turbine simulation, Wind Energy, 18, 699–711, https://doi.org/10.1002/we.1722, 2015. a
Technical University of Denmark: Sophia HPC Cluster, https://doi.org/10.57940/FAFC6M81, 2019. a
Tewari, M., Chen, F., Wang, W., Dudhi, J., LeMone, M., Mitchell, K., Ek, M., Gayno, G., Wegiel, J., and Cuenca, R. H.: Implementation and verification of the unified Noah land surface model in the WRF model, in: 20th Conference on Weather Analysis and Forecasting/16th Conference on Numerical Weather Prediction, Seattle, WA, USA, http://n2t.net/ark:/85065/d7fb523p (last access: 23 May 2023), 2004. a
Troldborg, N.: Actuator Line Modeling of Wind Turbine Wakes, PhD thesis, DTU, https://orbit.dtu.dk/files/5289075/niels_troldborg.pdf (last access: 23 May 2023), 2008. a
Troldborg, N. and Meyer Forsting, A. R.: A simple model of the wind turbine induction zone derived from numerical simulations, Wind Energy, 20, 2011–2020, https://doi.org/10.1002/we.2137, 2017. a
van der Laan, M. P. and Sørensen, N. N.: A 1D version of EllipSys, Tech. Rep. DTU Wind Energy E0141, Technical University of Denmark, https://orbit.dtu.dk/files/130854208/DTU_Wind_Energy_E_0141.pdf (last access: 23 May 2023), 2017a. a
van der Laan, M. P. and Sørensen, N. N.: Why the Coriolis force turns a wind farm wake clockwise in the Northern Hemisphere, Wind Energy Science, 2, 285–294, https://doi.org/10.5194/wes22852017, 2017b. a, b, c, d
van der Laan, M. P., Sørensen, N. N., Réthoré, P.E., Mann, J., Kelly, M., and Troldborg, N.: The k–ε–f_{P} model applied to double wind turbine wakes using different actuator disk force methods, Wind Energy, 18, 2223–2240, https://doi.org/10.1002/we.1816, 2015a. a, b, c
van der Laan, M. P., Sørensen, N. N., Réthoré, P.E., Mann, J., Kelly, M., Troldborg, N., Hansen, K. S., and Murcia, J. P.: The k–ε–f_{P} model applied to wind farms, Wind Energy, 18, 2065–2084, https://doi.org/10.1002/we.1804, 2015b. a, b
van der Laan, M. P., Sørensen, N. N., Réthoré, P.E., Mann, J., Kelly, M., Troldborg, N., Schepers, J. G., and Machefaux, E.: An improved k–ε model applied to a wind turbine wake in atmospheric turbulence, Wind Energy, 18, 889–907, https://doi.org/10.1002/we.1736, 2015c. a, b, c
van der Laan, M. P., Peña, A., Volker, P., Hansen, K. S., Sørensen, N. N., Ott, S., and Hasager, C. B.: Challenges in simulating coastal effects on an offshore wind farm, J. Phys.: Conf. Ser., 854, 012046, https://doi.org/10.1088/17426596/854/1/012046, 2017. a, b, c, d, e, f
van der Laan, M. P., Andersen, S. J., and Réthoré, P.E.: Brief communication: Windspeedindependent actuator disk control for faster annual energy production calculations of wind farms using computational fluid dynamics, Wind Energ. Sci., 4, 645–651, https://doi.org/10.5194/wes46452019, 2019. a
van der Laan, M. P., Kelly, M., Floors, R., and Peña, A.: Rossby number similarity of an atmospheric RANS model using limitedlengthscale turbulence closures extended to unstable stratification, Wind Energ. Sci., 5, 355–374, https://doi.org/10.5194/wes53552020, 2020. a, b, c, d, e, f, g
van der Laan, M. P., Andersen, S. J., Réthoré, P.E., Baungaard, M., Sørensen, J. N., and Troldborg, N.: Faster wind farm AEP calculations with CFD using a generalized wind turbine model, J. Phys.: Conf. Ser., 2265, 022030, https://doi.org/10.1088/17426596/2265/2/022030, 2022. a, b
van der Laan, M. P., GarcíaSantiago, O., Sørensen, N. N., Troldborg, N., Criado Risco, J., and J., B.: Simulating wake losses of the Danish Energy Island wind farm cluster, J. Phys.: Conf. Ser., accepted, 2023. a, b, c, d
Volker, P. J. H., Badger, J., Hahmann, A. N., and Ott, S.: The Explicit Wake Parametrisation V1.0: a wind farm parametrisation in the mesoscale model WRF, Geosci. Model Dev., 8, 3715–3731, https://doi.org/10.5194/gmd837152015, 2015. a, b, c, d, e
Wu, Y.T. and PortéAgel, F.: LargeEddy Simulation of WindTurbine Wakes: Evaluation of Turbine Parametrisations, Bound.Lay. Meteorol., 138, 345–366, https://doi.org/10.1007/s105460109569x, 2011. a
 Abstract
 Introduction
 Wind farm cluster validation case
 Simulation methodology
 Results and discussion
 Conclusions
 Appendix A: Similarity of the prescribed temperature inflow model
 Appendix B: Comparison and challenges of RANS inflow models applied to the wind farm cluster validation case
 Appendix C: Numerical setup of TurbOPark (as implemented in PyWake)
 Appendix D: Numerical setup of WRF simulations including wind farms
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Wind farm cluster validation case
 Simulation methodology
 Results and discussion
 Conclusions
 Appendix A: Similarity of the prescribed temperature inflow model
 Appendix B: Comparison and challenges of RANS inflow models applied to the wind farm cluster validation case
 Appendix C: Numerical setup of TurbOPark (as implemented in PyWake)
 Appendix D: Numerical setup of WRF simulations including wind farms
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References