the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Setpoint optimization in wind farms to mitigate effects of flow blockage induced by atmospheric gravity waves
Johan Meyers
Recently, it has been shown that flow blockage in large wind farms may lift up the top of the boundary layer, thereby triggering atmospheric gravity waves in the inversion layer and in the free atmosphere. These waves impose significant pressure gradients in the boundary layer, causing detrimental consequences in terms of a farm's efficiency. In the current study, we investigate the idea of controlling the wind farm in order to mitigate the efficiency drop due to windfarminduced gravity waves and blockage. The analysis is performed using a fast boundary layer model which divides the vertical structure of the atmosphere into three layers. The windfarm drag force is applied over the whole windfarm area in the lowest layer and is directly proportional to the windfarm thrust setpoint distribution. We implement an optimization model in order to derive the thrustcoefficient distribution, which maximizes the windfarm energy extraction. We use a continuous adjoint method to efficiently compute gradients for the optimization algorithm, which is based on a quasiNewton method. Power gains are evaluated with respect to a reference thrustcoefficient distribution based on the Betz–Joukowsky set point. We consider thrust coefficients that can change in space, as well as in time, i.e. considering timeperiodic signals. However, in all our optimization results, we find that optimal thrustcoefficient distributions are steady; any timeperiodic distribution is less optimal. The (steady) optimal thrustcoefficient distribution is inversely related to the vertical displacement of the boundary layer. Hence, it assumes a sinusoidal behaviour in the streamwise direction in subcritical flow conditions, whereas it becomes a Ushaped curve when the flow is supercritical. The sensitivity of the power gain to the atmospheric state is studied using the developed optimization tool for almost 2000 different atmospheric states. Overall, power gains above 4 % were observed for 77 % of the cases with peaks up to 14 % for weakly stratified atmospheres in critical flow regimes.
Today, it is well known that turbines strongly interact when clustered together in large arrays, increasing the momentum deficit in the lowest region of the atmospheric boundary layer (ABL). These turbine–turbine interactions, such as reduced wind speed and increased turbulence intensity, occur within the windfarm area and can lead to detrimental consequences in terms of a farm’s efficiency (Barthelmie et al., 2010). However, it has been recently discovered that nonlocal effects such as gravity waves may also have strong implications on the windfarm energy extraction (Allaerts and Meyers, 2018, 2019).
In a stable atmosphere, an air parcel which is vertically perturbed will have the tendency to fall back to its original position. In such a case, an oscillation is initiated that is driven by gravity and inertia; this is called a gravity wave. Mountains are examples of orographic obstacles that trigger vertical flow displacement and consequently gravity waves (Smith, 1980). The drag force exerted by the mountain is usually transported upward by these waves. At the point of breakdown, the drag force is released in the upper levels of the atmosphere, causing a slowdown of the largescale flow (Eliassen and Palm, 1960; Durran, 1990). Moreover, when air is lifted in a stable atmosphere, a cold anomaly is created, which induces horizontal pressure gradients (Smith, 2010).
In a wind farm, the upward displacement of the boundary layer, caused by diverging fluid streamlines due to flow deceleration by the turbines, can trigger gravity waves in the stable free atmosphere above the boundary layer as well (Smith, 2010; Allaerts and Meyers, 2017). As a result, an adverse pressure gradient develops in the induction region of the wind farm, which slows down the windfarm inflow velocity (Allaerts and Meyers, 2018). The size of this region scales with the length of the farm. This phenomenon is one possible cause of flow blockage. Note that it differs from classical hydrodynamic blockage caused by the turbine induction, which typically scales with the turbine rotor diameter and which has also been studied recently in much detail (Bleeg et al., 2018; Segalini and Dahlberg, 2019). The goal of the current study is to determine a windfarm thrustcoefficient distribution that minimizes the gravitywaveinduced blockage effects, maximizing the flow wind speed and therefore the power production. Moreover, we investigate the impact of different atmospheric conditions on the optimal thrustcoefficient distribution and corresponding power gains.
Gravity waves have been related to flows over mountains for a long time (Queney, 1948). However, the cumulated blockage effect induced by the wind farm in the induction region was associated with windfarminduced gravity waves only in recent years. In the pioneering work of Smith (2010), a quasianalytical model of atmospheric response to windfarm drag was used for modelling gravitywave excitation due to diverging streamlines above the windfarm area. Results have shown that gravitywave excitation is strongly dependent upon the height of the boundary layer and the stability of the atmosphere aloft. Later, a fast boundarylayer model was proposed by Allaerts and Meyers (2019), who highlighted the crucial role of the inversion layer in determining gravitywave patterns. The authors also used this model for an annual energy production study of the Belgian–Dutch offshore windfarm cluster, showing that the annual energy loss due to the effect of selfinduced gravity waves might be on the order of 4 % to 6 % (Allaerts et al., 2018).
Gravity waves were also observed in mesoscale and largeeddy simulation (LES) models. Fitch et al. (2012) and Volker (2014) proposed two different windfarm parameterizations for the Weather Research and Forecasting model (WRF). Windfarminduced gravity waves were observed in both cases, causing flow deceleration several kilometres upstream of the farm. Allaerts and Meyers (2017, 2018) have investigated the interaction between an “infinitely” wide wind farm and both a conventionally neutral and stable boundary layer in typical offshore conditions in a LES framework. They found that for low ABL heights, gravity waves induce strong pressure gradients and play an important role in the distribution of the kinetic energy within the farm. Wu and PortéAgel (2017) considered a large finitesize wind farm operating in a conventionally neutral boundary layer (CNBL) with different freeatmosphere stratification, and they conclude that strongly stratified atmospheres decrease the turbine power output up to 35 % with respect to the weakly stratified cases. Windfarm flow blockage was also detected in field measurements. Wind speed data taken before and after the placement of three wind farms showed that there was a reduction in wind speed of about 3 % in the induction region of each wind farm after turbines were installed (Bleeg et al., 2018).
In the last decades, a considerable amount of research has focused on windfarm control strategies that allow the maximization of the farm power output. We refer to Kheirabadi and Nagamune (2019) for a recent comprehensive overview. However, earlier studies all focus on influencing wake dynamics and wake mixing, which occur at a much smaller scale than windfarminduced gravity waves, to improve power extraction in waked turbines. Important control mechanisms include wake redirection (by yawing and tilting of the turbine), and turbine derating strategies. Control actions that influence windfarm physics on a much larger scale, such as selfinduced gravity waves, have not been explored to date.
In the current work, we concentrate on using windfarm control to alter/improve the interaction between the wind farm and its selfinduced gravitywave system. To this end, we use the fast boundarylayer model proposed by Allaerts and Meyers (2019), which divides the vertical structure of the atmosphere into three layers (from here on named the threelayer model in the paper), and we reformulate it as an optimization problem. The objective function is defined as the windfarm energy extracted over a time period T, while the constraints are the model equations plus a box constraint for the windfarm thrust setpoint distribution C_{T}(x, y, t). Note that we do not use the tipspeed ratio and/or the pitch angle distribution as control parameters. Instead, we directly control the thrust setpoint distribution. In fact, the former approach would not add further insight into the study performed in the current manuscript. The model equations are derived following the theory for interacting gravity waves and boundary layers developed by Smith et al. (2006) and Smith (2007, 2010). Consequently, the optimal thrustcoefficient distribution computed using the optimization formulation of the threelayer model takes into account the effects of selfinduced gravity waves. Hence, we investigate whether it is possible to mitigate gravitywaveinduced blockage effects by varying the thrust setpoint distribution within the windfarm area.
The remainder of this paper is formulated as follows. The threelayer model and its optimization formulation are introduced in Sect. 2. Next, Sect. 3 describes the numerical setup, windfarm layout and atmospheric state. Thereafter, Sect. 4 presents optimization results. The optimal thrust setpoint distributions obtained in two different flow cases are discussed in Sect. 4.1. The sensitivity of the power gain to the atmospheric state is carried out in Sect. 4.2. Finally, conclusions and suggestions for further research are given in Sect. 5.
We now introduce the approach used for modelling windfarminduced gravity waves and the method applied for maximizing the windfarm energy output. The threelayer model is described in Sect. 2.1, and its optimization formulation is derived in Sect. 2.2.
2.1 Threelayer model
In the work of Smith (2010), the atmospheric response to windfarm drag is simulated by dividing the vertical structure of the atmosphere in two layers: the ABL and the free atmosphere aloft. This approach has strong limitations. In fact, the author is implicitly assuming that the turbine drag is mixed homogeneously between turbine level and the top of the ABL. In real wind farms, the turbine drag slows down the flow only within a few hundred metres of the ground level, triggering the formation of an internal boundary layer (Wu and PortéAgel, 2013; Allaerts and Meyers, 2017). To overcome the limitations of Smith's model, the threelayer model divides the ABL into two layers: the windfarm layer in which the turbine forces are felt directly (a layer's height of twice the turbine hub height has been used by Allaerts and Meyers (2019) based on insights from LES in Allaerts et al., 2018) and a second layer up to the top of the ABL. Finally, the third layer models the free atmosphere above the ABL following the approach of Smith (2010).
The threelayer model has been validated against LES results by Allaerts and Meyers (2019) (see Sect. 3 VAL2) on a twodimensional (x–z) domain (i.e. all spanwise derivatives are set to zero). The model shows a mean absolute error (MAE) of 1.3 % and 1.8 % in terms of maximum displacement of the inversion layer and maximum pressure disturbance, respectively. Moreover, the model underestimates the velocity over the windfarm area with a MAE of 5.6 %. Note that the threelayer model is a linearized model; hence the discrepancies with LES results increase with increasing perturbation values. In fact, the model agrees very well with LES data when perturbations are small (i.e. when nonlinear effects are negligible). From this perspective, it may be expected that errors decrease slightly in optimized settings in which perturbation magnitudes are typically lower. For further details about model validation, we refer to Allaerts and Meyers (2019).
The model equations are derived starting from the incompressible threedimensional Reynoldsaveraged Navier–Stokes (RANS) equations for the ABL (Stull, 1988). A depth integration over the wind farm and upper layer height is further computed, which removes the vertical velocity from the equations. Hence, the basic equation system is reduced to a set of only three equations: the continuity equation and the momentum equations in horizontal directions. Subsequently, the governing equations are linearized with respect to the background state variables, using some additional modelling assumptions for the turbulent stresses (see Allaerts and Meyers (2019) for more details). As a result, we use the following equations for the two layers in the ABL.
When the wind farm is not operating (i.e. the windfarm drag force is zero), a horizontally invariant reference state of (U_{1}, H_{1}) and (U_{2}, H_{2}) characterizes the wind farm and upper layer, where U_{1}=(U_{1}, V_{1}) and U_{2}=(U_{2}, V_{2}) are the heightaveraged horizontal components of the background velocity and H_{1}, H_{2} represent the height of the two layers. Whenever the farm extracts power from the flow, small velocity and height perturbations (u_{1}, η_{1}) and (u_{2}, η_{2}) are triggered. The equations derived by Allaerts and Meyers (2019) predict the spatial evolution of these perturbations. In this article, we also consider the temporal evolution, and thus, the relevant time derivatives are added to the equations. Furthermore, ρ_{0} denotes the air density, assumed constant within the ABL; ν_{t,1} and ν_{t,2} are the depthaveraged turbulent viscosity; f_{c}=2Ωsin ϕ is the Coriolis frequency, with Ω the angular velocity of the earth and ϕ the latitude; and $\mathbf{J}={\mathit{e}}_{x}\otimes {\mathit{e}}_{y}{\mathit{e}}_{y}\otimes {\mathit{e}}_{x}$ is the twodimensional rotation dyadic with e_{x} and e_{y} twodimensional unit vectors in the x and y directions, respectively. Finally, the perturbation of the friction at the ground and at the interface between both layers is described by the matrices C^{′} and D^{′}.
The righthand side of Eq. (1) is characterized by the windfarm drag force f. We use a boxfunction windfarm force model similar to that of Smith (2010) in our study. This allows us to avoid the complexity of wake models while gaining in computational time. In fact, this model uniformly spreads the force over the simulation cells in the windfarm area and does not represent the disturbances caused by each turbine in detail. The force magnitude depends on the windfarm layout, the wind speed and the thrust setpoint distribution (i.e. the C_{T} value in every grid cell within the farm). As for the flow equations, the windfarm drag force model is linearized around a constant background state. We retain the first two terms of the Taylor expansion; both scale linearly with the thrustcoefficient distribution. Hence, the drag force is given by $\mathit{f}={\mathit{f}}^{\left(\mathrm{0}\right)}+{\mathit{f}}^{\left(\mathrm{1}\right)}$ with
where
and with B(x, y) a box function equal to 1 within the windfarm area and zero outside. The x and y axes denote the streamwise and spanwise directions, respectively. The windfarm drag force magnitude in Eqs. (5) and (6) scales with
where s_{x} and s_{y} are the streamwise and spanwise turbine spacings relative to the rotor diameter, η_{w} is the wake efficiency and $\mathit{\gamma}={u}_{\mathrm{r}}^{\mathrm{2}}/\parallel {\mathit{U}}_{\mathrm{1}}{\parallel}^{\mathrm{2}}$ is a velocity shape factor with u_{r} the rotoraveraged wind speed (Allaerts and Meyers, 2018). Moreover, $\mathbf{I}={\mathit{e}}_{x}\otimes {\mathit{e}}_{x}+{\mathit{e}}_{y}\otimes {\mathit{e}}_{y}$ denotes the unit dyadic. Finally, C_{T}(x, y, t) represents the thrustcoefficient distribution. To compute the thrust coefficient ${\stackrel{\mathrm{\u0303}}{C}}_{\mathrm{T},k}\left(t\right)$ of a turbine at location (x_{k}, y_{k}), it is possible to evaluate the thrust setpoint distribution C_{T}(x_{k}, y_{k}, t). A more accurate connection between ${\stackrel{\mathrm{\u0303}}{C}}_{\mathrm{T},k}\left(t\right)$ and the drag force f would, for example require the use of an analytical windfarm wake model. This is however not considered in the current work, so that wake effects are not explicitly incorporated in the optimization. Rather, we consider the optimization of the gravitywave system, while presuming that the wake efficiency parameter η_{w} does not change as a result of the optimization. Relation (Eq. 6) is nonlinear since it contains a product between time and spacedependent variables (i.e. C_{T} and u_{1}). We decide to retain this term because it allows us to include gravitywave feedback on windfarm energy extraction. In fact, f^{(1)}≥0 so that it reduces the drag force that the farm exerts on the flow, thereby reducing effects of blockage in the model. We note that Allaerts and Meyers (2019) have shown that the flow perturbations computed with this simple drag force model have similar trends and orders of magnitude as the ones computed using a drag model that relies on the more detailed analytical wake model of Niayifar and PortéAgel (2016). Therefore, we believe that the model adopted is a reasonable representation of reality.
The total vertical displacement of the inversion layer ${\mathit{\eta}}_{\mathrm{t}}={\mathit{\eta}}_{\mathrm{1}}+{\mathit{\eta}}_{\mathrm{2}}$ triggers gravity waves which induce pressure perturbations p. The relation between these two quantities is given by Smith (2010):
where ℱ^{−1} and ∗ denote the inverse Fourier transform and the convolution product, respectively. The pressure p is evaluated at the capping inversion height, and it is assumed to be constant through the whole ABL (using the classical boundarylayer approximation $\partial p/\partial z=\mathrm{0}$). The complex stratification coefficient $\widehat{\mathrm{\Phi}}$ in Fourier components is expressed as
Relation (Eq. 10) is obtained from linear threedimensional, nonrotating, nonhydrostatic gravitywave theory (Nappo, 2002) under the assumption of constant wind speed U_{g}=(U_{g}, V_{g}) and Brunt–Väisälä frequency N. The reduced gravity ${g}^{\prime}=g\mathrm{\Delta}\mathit{\theta}/{\mathit{\theta}}_{\mathrm{0}}$ accounts for twodimensional trapped lee waves (from here on named inversion waves) which corrugate the capping inversion layer. The potentialtemperature difference Δθ denotes the strength of the capping inversion, and θ_{0} is a reference potential temperature. The effect of internal gravity waves is represented by the second term of relation (Eq. 10), where m denotes the vertical wavenumber which is given by
According to the sign of m^{2}, we can have propagating or evanescent waves. Moreover, $\mathrm{\Omega}=\mathit{\omega}\mathit{\kappa}\cdot {\mathit{U}}_{\mathrm{g}}$ denotes the intrinsic wave frequency with κ=(k, l) the horizontal wavenumber vector.
Finally, combining Eq. (9) with Eqs. (2) and (4), we can write the continuity equations for the wind farm and upper layer as
where $p={p}_{\mathrm{1}}+{p}_{\mathrm{2}}$ is intended as the sum of the pressure perturbations induced by the vertical displacements η_{1} and η_{2}, respectively. This form will be used in the remainder of the manuscript.
2.2 Optimization model
The goal of the optimization framework is to find a timeperiodic optimal thrustcoefficient distribution ${C}_{\mathrm{T}}^{\mathrm{O}}(x$, y, t) that minimizes the gravitywaveinduced blockage effects, maximizing the flow wind speed and consequently the windfarm energy extraction over a selected time period T. The background atmospheric state is presumed to be steady, which is the reason why we use a timeperiodic control (i.e. leading to a moving time average of the optimal control that is steady and does not lead to endoftime effects). The windfarm layout and the atmospheric state are inputs of the optimization model and are detailed in Sect. 3. Note that the relation between overall windfarm drag and windfarm blockage is nontrivial. On the one hand, increased windfarm drag leads to increased windfarm blockage induced by gravity waves. This results from mass conservation and the upward displacement of the free atmosphere. On the other hand, increased windfarm blockage reduces windfarm drag. Thus, the aim of the optimization is to find the optimal balance between these two opposing trends.
By using axial momentum theory (Burton et al., 2001), we find that the power coefficient C_{p}(x, y, t) depends upon the thrust coefficient according to the following nonlinear relationship:
The objective function of the optimization model consists in the energy extracted by the farm in the time interval [0, T]; hence it is defined as
where $\mathrm{\Omega}={D}_{x}\times {D}_{y}$ is the computational domain area. The nonlinear relationship between C_{p} and C_{T} and the product between control and state variables in Eq. (15) imply that the objective function 𝒥 is nonconvex.
The windfarm optimal configuration that maximizes the energy output (note that the objective function is defined with a minus sign) is then obtained by solving the following nonlinear timeperiodic optimization problem constrained by a system of partial differential equations (PDEs):
The constraints are the state (or forward) equations presented in the previous paragraph. Since Eq. (14) is defined only for C_{T}∈[0, 1), we added a box constraint to the optimization model. Moreover, the time periodicity is imposed by assuming C_{T}(x, y, 0)=C_{T}(x, y, T). The system state ψ=[u_{1}, v_{1}, u_{2}, v_{2}, p_{1}, p_{2}] includes the velocity and pressure perturbations in the wind farm and upper layer, which also define the unknowns of the threelayer model. The control parameters consist of the value of the thrust set point in each grid cell within the windfarm area. Hence, the size of the control space is proportional to ${N}_{x}^{\mathrm{wf}}{N}_{y}^{\mathrm{wf}}{N}_{t}$, where N_{t} represents the number of time steps within the time horizon T, while ${N}_{x}^{\mathrm{wf}}$ and ${N}_{y}^{\mathrm{wf}}$ denote the number of grid points within the windfarm area along the x and y directions, respectively.
It is common practice in a PDEconstrained optimization problem to not optimize the cost functional 𝒥(ψ, C_{T}) directly because such a problem would span both the state and control space. To avoid exploring the entire feasibility region, we require ψ(C_{T}) to be the solution of the state equations throughout the optimization process. In other words, defining an operator 𝒩(ψ, C_{T}) that denotes the state equations, we are enforcing 𝒩(ψ(C_{T}), C_{T})=0 during optimization iterations. This technique leads us to a reduced optimization problem with feasibility region given by the control space (De Los Reyes, 2015). The reduced optimization problem is written as
where the only remaining constraints are the ones on the control parameters.
The gradient of the reduced objective function $\mathrm{\nabla}\stackrel{\mathrm{\u0303}}{\mathcal{J}}$ is needed for the solution of the reduced optimization problem. To this end, we use the continuous adjoint method. The adjoint (or backward) equations are given by (see Appendix A2 for detailed derivation)
Note that the minus sign in the argument of ${\mathcal{F}}^{\mathrm{1}}\left(\widehat{\mathrm{\Phi}}\right)(\mathit{x}$, −t) is not a result of classical integration by parts, but rather arrives from applying Fubini's theorem to the convolution term in Eqs. (12) and (13) (see Appendix A2 for details). The adjoint variables are grouped in the vector ${\mathit{\psi}}^{*}=[{\mathit{\zeta}}_{\mathrm{1}}$, ζ_{2}, Π_{1}, Π_{2}], where ζ_{1}=(ζ_{1}, χ_{1}) and ζ_{2}=(ζ_{2}, χ_{2}) are the adjoint velocity perturbation fields in the wind farm and upper layer, respectively, while Π_{1} and Π_{2} are the adjoint pressure perturbations. Using the solution of the adjoint equations, the gradient of the cost function is expressed as (see Appendix A3 for details)
where $\mathrm{d}{C}_{\mathrm{p}}/\mathrm{d}{C}_{\mathrm{T}}$ is computed from Eq. (14). To compute the gradient $\mathrm{\nabla}\stackrel{\mathrm{\u0303}}{\mathcal{J}}$, we need to solve the forward and backward equations. Since the cost for solving the adjoint equations is roughly the same as for the forward equation, the computational cost for evaluating $\mathrm{\nabla}\stackrel{\mathrm{\u0303}}{\mathcal{J}}$ is proportional to the cost of solving twice the state equations. To verify the approach, we compare the adjoint gradient to a standard finitedifference approximation in Appendix A4.
We define the model setup used to assess the potential of setpoint optimization in mitigating selfinduced gravitywave effects in this section. We discuss the numerical setup in Sect. 3.1. Next, the selected windfarm layout is presented in Sect. 3.2. Finally, the atmospheric state is discussed in Sect. 3.3.
3.1 Numerical setup
Both the forward and adjoint equations are discretized using a Fourier–Galerkin spectral method in space and time. The use of Fourier modes in time automatically results in satisfying the periodicity conditions that we are aiming for in our optimization setup. All terms in the equations are linear, except for the firstorder term of the windfarm drag force (and its adjoint). These terms contain products between temporarily and spatially dependent variables. To avoid expensive convolution sums, these products are computed in physical space. Full dealiasing is obtained by padding and truncation according to the $\mathrm{3}/\mathrm{2}$ rule (Canuto et al., 1988). The use of Fourier modes in space forces periodic boundary conditions at the edges of the computational domain. Therefore, the domain has a sufficiently large dimension ${D}_{x}\times {D}_{y}=\mathrm{1000}\times \mathrm{400}$ km^{2}, so that the perturbations die out before being recycled. The grid has ${N}_{x}\times {N}_{y}=\mathrm{4000}\times \mathrm{1600}$ grid points, which corresponds to a space resolution of Δ=250 m or 6.4×10^{6} DOF per layer. Finally, different time horizon values are used spanning from T=10 min to T=10 h with the number of time steps ranging from N_{t}=12 to N_{t}=120. The discretized forward and backward equations form two systems of the dimension 6N_{x}N_{y}N_{t}×6N_{x}N_{y}N_{t}, which are solved using the LGMRES algorithm (Baker et al., 2005).
For the optimization, two different algorithms are compared in Fig. 1. The LBFGSB (limitedmemory Broyden–Fletcher–Goldfarb–Shanno with box constraint) algorithm (Byrd et al., 1995) is an iterative quasiNewton method. In the current application, the step length is evaluated with the inexact line search Wolfe condition (Wolfe, 1969). The truncated Newton method (TNC) computes the search direction by solving the Newton equation iteratively, applying the conjugate gradient method. This inner loop is stopped (truncated) as soon as a termination criterion is satisfied (Nocedal and Wright, 1999). In both cases, the system matrix of the Newton equation consists of an approximate Hessian matrix, while the righthand side needs gradient information to be computed, which is provided by the continuous adjoint method (see Appendix A for derivation and validation). Figure 1 shows that the cost function decreases rapidly in the first two to three algorithm iterations, reaching a plateau afterwards. The use of a quasiNewton method in combination with the limited complexity of our optimization model (for instance, the constraints are linearized equations) allow us to reach such a fast convergence. Moreover, the continuous adjoint method limits the number of function evaluations, since it is not necessary to evaluate $\stackrel{\mathrm{\u0303}}{\mathcal{J}}({C}_{\mathrm{T}}+\mathit{\alpha}\mathit{\delta}{C}_{\mathrm{T}})$ for all directions δC_{T} in the control space (at the expense of solving an auxiliary set of equations). In particular, Fig. 1a shows that the cost functional is converged after six LBFGSB iterations. Apart from the first iteration, the line search method needs three “function evaluations” before updating the cost functional. Hence, we need to solve 20 times the forward and backward equations for reaching convergence. On the other hand, Fig. 1b illustrates that the cost functional is mainly minimized within the first TNC iteration, and convergence is reached after only four iterations. However, 63 function evaluations are needed. Hence, we will use the LBFGSB algorithm for solving the PDEconstrained optimization problem in the remainder of the article. To limit computational effort, a maximum of four LBFGSB iterations will be performed.
The solver (which is not parallelized) takes a couple of hours to solve the equations for a grid with resolution of 250 m (6.4×10^{6} DOF per layer). Since convergence is reached after approximately 20 function evaluations (which means that we solve state and adjoint equations 20 times), the optimizer takes a couple of days to compute an optimal thrust setpoint distribution. However, after this work was performed, we have upgraded the forward solver which is now approximately 1000 times faster than our previous version. Optimization of the backward solver is planned for the future, and we expect that this will lead to an optimization algorithm that will only take several minutes for the same case.
3.2 Windfarm layout
Allaerts and Meyers (2019) conducted a sensitivity study on the effects of windfarm layout on gravitywaveinduced power losses. They show that these power losses monotonically increase with the size of the farm. Also, they mention that the losses are at a maximum when the windfarm ratio ${L}_{y}/{L}_{x}$ is close to $\mathrm{3}/\mathrm{2}$, while being negligible for a very wide but short farm, and vice versa (i.e. negligible for ${L}_{y}/{L}_{x}\ll \mathrm{1}$ and ${L}_{x}/{L}_{y}\gg \mathrm{1}$). Since we are interested in optimal thrustcoefficient distributions in the presence of gravity waves, we have selected the “worstcase” windfarm layout (i.e. a windfarm width and length of L_{y}=30 km and L_{x}=20 km, respectively). We note that this was also the farm layout chosen by Allaerts et al. (2018) and Allaerts and Meyers (2019), which in size resembles the Belgian–Dutch windfarm offshore cluster located in the North Sea, but simplified to a rectangularshaped wind farm. Also, Smith (2010), Fitch et al. (2012) and Wu and PortéAgel (2017) have used a farm with similar dimensions in their studies. The wind turbine relative spacings along the x and y directions are ${s}_{x}={s}_{y}=\mathrm{5.61}$ (both nondimensionalized with respect to the turbine rotor diameter D), so that the density of turbines in the farm is similar to the one of Allaerts and Meyers (2019) (i.e. leading to β=0.01 in Eq. (8), setting both the wake efficiency η_{w} and γ to 0.9 as in Allaerts and Meyers, 2018). Note that we do not define a specific layout or a number of turbines, but we only fix the density of turbines in the farm. The turbine dimensions are based on a DTU 10 MW IEA wind turbine (Bortolotti et al., 2019) with rotor diameter D=198 m and turbine hub height z_{h}=119 m.
3.3 Background state variables
The governing equations are linearized around a constant background state. To determine this state, we need vertical profiles of potential temperature, velocity, shear stress and eddy viscosity plus the surface roughness z_{0} and the friction velocity u_{*}. We describe the techniques used in determining these profiles in the remainder of this section. Similar to Allaerts and Meyers (2019), we select two atmospheric states for initial testing of the optimizer, corresponding to sub and supercritical flow conditions.
We choose a temperature profile that corresponds to a conventionally neutral ABL. The potential temperature in the neutral ABL is fixed to θ_{0}=288.15 K. A capping inversion strength Δθ of 5.54 and 3.7 K is used, which leads to a sub and supercritical flow, respectively (see below). Finally, a freeatmosphere lapse rate Γ=1 K km^{−1} is chosen. The lapse rate also defines the Brunt–Väisälä frequency N.
The velocity and stress profiles within the ABL are obtained following the approach of Nieuwstadt (1983). The nondimensional surface roughness length ${\stackrel{\mathrm{\u203e}}{z}}_{\mathrm{0}}={z}_{\mathrm{0}}/H$ and the nondimensional boundarylayer height ${h}_{*}=H{f}_{\mathrm{c}}/{u}_{*}$ are the input parameters of Nieuwstadt model, where f_{c} is the Coriolis frequency and $H={H}_{\mathrm{1}}+{H}_{\mathrm{2}}$ is the ABL height. The windfarm layer height is assumed to be twice the turbine hub height, so H_{1}=2z_{h}. The ABL height is fixed to H=1000 m and the friction velocity is set to ${u}_{*}=\mathrm{0.6}$ m s^{−1}. Finally, a surface roughness of ${z}_{\mathrm{0}}={\mathrm{10}}^{\mathrm{1}}$ m is adopted. Using ${h}_{*}=\mathrm{0.166}$ and ${\stackrel{\mathrm{\u203e}}{z}}_{\mathrm{0}}={\mathrm{10}}^{\mathrm{4}}$ as input values for the Nieuwstadt model, we derive the velocity U_{1}, U_{2}, the eddy viscosity ν_{t,1}, ν_{t,2} for the windfarm and upper layer, and the friction coefficients C and D (used for computing the matrices C^{′} and D^{′}; see Allaerts and Meyers, 2019). Besides the friction coefficients C and D, which are given at z=0 and z=H_{1}, all other physical quantities are depthaveraged over the height H_{1} and H_{2}. Finally, the wind profile is oriented such that the wind in the windfarm layer is always directed along the x axis (i.e. V_{1}=0 m s^{−1}).
The pressure gradient strengths induced by inversion and internal gravity waves are dependent upon the Froude number $Fr={U}_{\mathrm{B}}/\sqrt{{g}^{\prime}H}$ and a nondimensional group ${P}_{\mathrm{N}}={U}_{\mathrm{B}}^{\mathrm{2}}/NH\parallel {\mathit{U}}_{\mathrm{g}}\parallel $, respectively (Smith, 2010; Allaerts and Meyers, 2019), where the velocity scale U_{B} is defined as
The chosen background state defines a Froude number of 0.9 for Δθ=5.54 K, which implies subcritical flow conditions (Fr<1), and a Froude number of 1.1 for Δθ=3.7 K, which leads to supercritical flow conditions (Fr>1). Further, P_{N} expresses the impact of internal waves in the troposphere, which increases when P_{N} decreases. The background state defined in Table 1 leads to P_{N}=1.92. The numerical setup, windfarm layout and background state variables are summarized in Table 1.
We discuss the results of the optimization problem in the current section. Firstly, the optimal thrustcoefficient distributions and relative power gains are illustrated for two specific flow conditions in Sect. 4.1. Thereafter, the sensitivity of the power gain to the atmospheric state is studied in Sect. 4.2.
4.1 Optimal thrustcoefficient distributions
The optimization model described in Sect. 2.2 is time and space dependent. Hence, the model is capable of finding a timeperiodic optimal thrustcoefficient distribution over the windfarm area in a fixed time interval [0, T]. However, all optimal thrust setpoint distributions found for the different combinations of time horizons and time steps reported in Table 1 are constant in time. We have verified this using a range of steady and unsteady starting conditions for C_{T} in the algorithm, but we did not find any unsteady optimum. We believe that this is due to the use of steadystate inflow conditions, meaning that we neglect mesoscale temporal variations in the velocity field (these could lead to timedependent optimal control signals, but are not included in the current work). Since we do not observe any unsteady behaviour in our optimal solutions, we show only steadystate results in the remainder of the manuscript, and we conclude for the time being that unsteady timeperiodic excitation is less effective than a stationary spatially optimal distribution in this context.
We also note that our findings are in contrast with recent works of Goit and Meyers (2015), Munters and Meyers (2018), and Frederik et al. (2020), in which the authors illustrated the benefits of dynamic induction control over yaw and static induction control. However, the characteristic timescale of gravitywave effects is estimated to be approximately 1 h (Gill, 1982; Allaerts and Meyers, 2019), which is an order of magnitude above the typical timescale of wake convection between turbines, and turbulent mixing in turbine wakes (this also justifies the larger sampling time used). Hence, while unsteadiness of the thrust coefficient (with a typical timescale of 50 s for largescale turbines) can lead to improved wake mixing (Goit and Meyers, 2015; Munters and Meyers, 2018; Frederik et al., 2020), it has no impact on phenomena that occur at larger timescales, such as windfarminduced gravity waves.
The steadystate optimal thrustcoefficient distributions obtained in sub and supercritical conditions are analysed in the remainder of this section. To improve the understanding of such distributions, gravitywaveinduced flow patterns obtained with ${C}_{\mathrm{T}}^{\mathrm{O}}(x$, y) are compared with a reference case. The setup of the reference model is the one reported in Table 1, but instead a uniform thrust setpoint distribution over the windfarm area is used, with ${C}_{\mathrm{T}}^{\mathrm{R}}(x$, $y)={C}_{\mathrm{T}}^{\mathrm{Betz}}=\mathrm{8}/\mathrm{9}$.
Figure 2 illustrates a planform view of the perturbation flow patterns obtained with Fr=0.9 (top row) and Fr=1.1 (bottom row) using the reference model setup. The farm extracts energy from the flow, causing a momentum sink in the windfarm layer. Due to the continuity constraint, an upward flow displacement above the windfarm area takes place, which causes the boundary layer height to increase. Figure 2a shows that an inversionlayer vertical displacement of about 65 m takes place at the windfarm entrance region for the subcritical case. A second peak of lower magnitude is located in the downwind region. On the other hand, for the supercritical case, Fig. 2d displays a similar maximum value of η_{t} attained close to the windfarm centre. In both cases, the inversionlayer vertical displacement decreases in the windfarm exit region and assumes a wavy behaviour in the windfarm wake. The vertical displacement of air parcels triggers inversion waves on the 2D inversionlayer surface and internal waves in the free atmosphere (3D waves). These waves induce pressure gradients, as visible in Fig. 2b and e, where a region of high pressure builds up in correspondence with high η_{t} values, leading to flow blockage. However, Fig. 2b shows a stronger adverse pressure gradient in the windfarm induction region than the one in Fig. 2e. In fact, inversion waves travel upstream in subcritical conditions, which leads to more slowdown in the induction region. In both the sub and supercritical cases, favourable pressure gradients reduce the velocity deficits in the bulk of the farm. Finally, Fig. 2c and f illustrate relative velocity reductions in the windfarm layer. The stronger inversion strength found in the subcritical flow case transforms the inversion layer in a quasirigid lid, which limits vertical displacements. The lower streamlines' divergence over the windfarm area implies lower velocity reductions. Moreover, the favourable pressure gradient is stronger when Fr=0.9, allowing for lower velocity deficits within the windfarm area. This explains the higher velocity reduction (up to 20 %) seen in Fig. 2f. Such a strong response could be on the limit of our small amplitude assumption. The planform view of pressure and velocity perturbations in the windfarm and upper layer in subcritical flow conditions are also illustrated on a wider domain in Appendix A (see Fig. A1).
The goal of our study is to find an optimal setpoint distribution which reduces the velocity perturbations displayed in Fig. 2c and f. While maximizing the flow wind speed through the farm, we also maximize the windfarm energy extraction. To this end, we solve the optimization problem discussed in Sect. 2.2. The inputs of the optimization model are the windfarm layout and the atmospheric conditions, which are detailed in Table 1. Moreover, an initial thrustcoefficient distribution needs to be specified. We have verified that for many different initial conditions the algorithm always converges to the same optimal solution. Therefore, a random initial thrust setpoint distribution is chosen. The optimal configurations obtained for different Froude numbers are illustrated in Fig. 3. We find that the optimal thrustcoefficient distributions are nonuniform in space and assume different spatial distributions according to the atmospheric state. In particular, when the flow is subcritical the optimal thrust setpoint distribution assumes a sinusoidal behaviour in the streamwise direction while it becomes a Ushaped curve when the flow is supercritical. In both cases, ${C}_{\mathrm{T}}^{\mathrm{O}}$ is almost invariant along the spanwise direction.
We denote with ${\mathcal{P}}^{\mathrm{R}}=\stackrel{\mathrm{\u0303}}{{\mathcal{J}}^{\mathrm{R}}}/T$ and ${\mathcal{P}}^{\mathrm{O}}={\stackrel{\mathrm{\u0303}}{\mathcal{J}}}^{\mathrm{O}}/T$ the power extracted using ${C}_{\mathrm{T}}={C}_{\mathrm{T}}^{\mathrm{R}}=\mathrm{8}/\mathrm{9}$ and ${C}_{\mathrm{T}}={C}_{\mathrm{T}}^{\mathrm{O}}$, respectively. Further, we define
where 𝒢 denotes the relative power gain obtained using an optimal thrustcoefficient distribution instead of the reference one. Note that the optimal distributions are steady state; therefore the power gain definition is not dependent on the choice of the time horizon T. The power gains attained in sub and supercritical flow conditions are 5.3 % and 7 %, respectively. Clearly, power gains are also strongly dependent on the atmospheric conditions. Therefore, a sensitivity study is carried out in Sect. 4.2. To assess the benefit of an optimal nonuniform distribution over an optimal uniform one, we have applied the optimization framework developed in Sect. 2.2 assuming a spatially invariant C_{T}. Results are discussed in Appendix B.
The optimal setpoint distributions displayed in Fig. 3 are related to the vertical displacement of the inversion layer over the windfarm area. Figure 4 shows streamwise profiles of η_{t} and ${C}_{\mathrm{T}}^{\mathrm{O}}$ through the centre of the farm for Fr=0.9 and Fr=1.1. To reduce gravitywave excitation, ${C}_{\mathrm{T}}^{\mathrm{O}}$ is seen to be inversely related with η_{t}. In fact, Fig. 4a shows that the streamwise profile of η_{t} has a sinusoidal behaviour. Hence, the optimal setpoint distribution is sinusoidal as well, explaining the pattern displayed in Fig. 3a. On the other hand, η_{t} assumes a Ushaped profile through the wind farm in supercritical conditions (see Fig. 4b), a profile that is also found in ${C}_{\mathrm{T}}^{\mathrm{O}}$ (see Fig. 3b). Moreover, Fig. 2a and d show that the gradient of η_{t} along the spanwise direction is much smaller than the one along the streamwise direction, explaining the almost constant thrust setpoint distributions along the y direction. Figure 4 also shows that ${\mathit{\eta}}_{\mathrm{t},\mathrm{max}}^{\mathrm{O}}<{\mathit{\eta}}_{\mathrm{t},\mathrm{max}}^{\mathrm{R}}$ in both sub and supercritical conditions, meaning that the optimal thrust setpoint distribution decreases the upward flow displacement over the windfarm area. The maximum inversionlayer displacement is located at the entrance region of the farm. If we compare ${\mathit{\eta}}_{\mathrm{t}}^{\mathrm{R}}$ and ${\mathit{\eta}}_{\mathrm{t}}^{\mathrm{O}}$ in this region, a displacement reduction of 14.5 % and 16.8 % is attained with the optimal configuration for the sub and supercritical cases, respectively.
A lower vertical displacement of the inversion layer reduces gravitywave excitation; therefore we also expect a lower strength of the adverse pressure gradient at the entrance of the farm compared to the one obtained with ${C}_{\mathrm{T}}^{\mathrm{R}}$. Figure 5a and c confirm this hypothesis, showing streamwise profiles of pressure perturbations p^{R} and p^{O} through the centre of the farm for Fr=0.9 and Fr=1.1. The pressure peak is located at the entrance of the farm and a pressure peak reduction of 14.3 % and 16.2 % is attained with the optimal configuration for the sub and supercritical cases, respectively. Figure 5b and d show streamwise profiles of velocity perturbations ${u}_{\mathrm{1}}^{\mathrm{R}}$ and ${u}_{\mathrm{1}}^{\mathrm{O}}$ through the centre of the farm for Fr=0.9 and Fr=1.1. The lower adverse pressure gradient strength attained with the optimal configuration allows for a lower velocity perturbation u_{1} in the induction region with respect to the reference case. Moreover, the optimal configuration also reduces the streamline divergence, accounting for higher flow wind speeds through the farm. Consequently, a velocity perturbation reduction of 13.4 % and 15.5 % is attained for the sub and supercritical cases, which explains the higher power gain obtained for Fr=1.1. The relative change in percentage between optimal and reference maximum flow perturbation values is summarized in Table 2.
The optimal thrustcoefficient distributions and power gains discussed in this section are obtained with data listed in Table 1. However, the atmospheric state changes in real case scenarios, and we have seen that the optimal configuration strongly depends upon the atmospheric parameters. Therefore, the sensitivity of the power gain to the atmospheric state is performed in the next section.
4.2 Sensitivity study
Allaerts and Meyers (2019) pointed out that gravitywaveinduced power loss is significant only for certain atmospheric states. Since our aim is to recover its power loss, we also expect the power gain to be sensitive to the atmospheric conditions. We note that gravitywave patterns are also sensitive to the windfarm layout. However, a sensitivity study over the windfarm layout is beyond the scope of the article.
The nondimensionalization of the threelayer model equations with respect to the boundary layer height H and the friction velocity u_{*} highlights four nondimensional groups that govern the atmospheric state, which are as follows.

The nondimensional boundary layer height ${h}_{*}=H{f}_{\mathrm{c}}/{u}_{*}$. Values of ${h}_{*}\approx \mathrm{0.1}$ denote shallow boundary layers typically found over sea, while ${h}_{*}\approx \mathrm{0.35}$ rather relates to a deep landbased boundary layer. We vary h_{*} between 0.16 and 0.4.

The nondimensional surface roughness length ${\stackrel{\mathrm{\u203e}}{z}}_{\mathrm{0}}={z}_{\mathrm{0}}/H$. This number varies several orders of magnitude according to the sea state or land surface. We vary ${\mathrm{log}}_{\mathrm{10}}\left({\stackrel{\mathrm{\u203e}}{z}}_{\mathrm{0}}\right)$ between −4.2 and −2.8 in the current study.

The nondimensional Brunt–Väisälä frequency $N/{f}_{\mathrm{c}}$. The Brunt–Väisälä frequency is an important parameter in gravitywave theory which expresses the highest possible frequency for internal gravity waves (Gill, 1982). Typical values of freeatmosphere lapse rate Γ range between 1 and 10 K km^{−1}. Low and high Γ values are associated with weakly and strongly stratified atmospheres, respectively. We vary Γ between 0.03 and 12 K km^{−1} corresponding to $\mathrm{10}\le N/{f}_{\mathrm{c}}\le \mathrm{200}$.

The inversion parameter ${g}^{\prime}H/A{u}_{*}^{\mathrm{2}}$. According to Csanady (1974), the height of the inversion layer is determined by a balance of surface stress and buoyancy. Equilibrium conditions are reached when ${g}^{\prime}H/A{u}_{*}^{\mathrm{2}}\approx \mathrm{1}$, with A=500 being an empirical constant. We vary the inversion parameter between 0.5 and 1.5.
Allaerts and Meyers (2019) conducted a similar sensitivity study on the gravitywaveinduced power loss on a wider range of nondimensional numbers. However, since we are optimizing turbine thrust set points, we need to ensure that ${U}_{\mathrm{1}}/{U}_{\mathrm{r}}<\mathrm{1}$ (U_{r}=11 m s^{−1} is the rated wind speed of the DTU 10 MW IEA wind turbine), otherwise turbines would operate in an aboverated wind speed regime and it would not make any sense to optimize their power production. The choice of the four ranges for the nondimensional groups discussed above ensures that ${U}_{\mathrm{1}}/{U}_{\mathrm{r}}\le \mathrm{0.9}$ for all atmospheric states.
Using the atmospheric state reported in Table 1, the nondimensional numbers assume values ${h}_{*}=\mathrm{0.166}$, ${\stackrel{\mathrm{\u203e}}{z}}_{\mathrm{0}}={\mathrm{10}}^{\mathrm{4}}$ and $N/{f}_{\mathrm{c}}=\mathrm{58}$. The inversion parameter is equal to 1.046 and 0.691 in the sub and supercritical cases, respectively. The optimal thrustcoefficient distributions discussed in Sect. 4.1 were obtained using these dimensionless group values. The sensitivity of the power gain to atmospheric conditions is performed by varying h_{*}, ${\stackrel{\mathrm{\u203e}}{z}}_{\mathrm{0}}$ and $N/{f}_{\mathrm{c}}$ against the inversion parameter ${g}^{\prime}H/A{u}_{*}^{\mathrm{2}}$, similarly to Allaerts and Meyers (2019). The numerical setup is the one detailed in Table 1. However, we use a grid cell size which is 4 times bigger ($\mathrm{\Delta}x\times \mathrm{\Delta}y=\mathrm{1000}\times \mathrm{1000}$ m^{2}), meaning that we use 4×10^{5} cells instead of 6.4×10^{6}, so that the necessary computational resources remain reasonable. To assess the validity of this choice, we performed a grid sensitivity study in Appendix C showing that the power gain value changes about 1 % when the number of grid cells is increased by 1 order of magnitude (see Fig. C1). The high computational efficiency of the threelayer model allowed us to perform a sensitivity study of the optimization results over 1960 different atmospheric conditions (thus effectively running an optimization problem for every atmospheric state). Since the windfarm layout impact on power gains is beyond the scope of our study, we impose the wind direction to be along the x axis in the windfarm layer in all simulations (V_{1}=0 m s^{−1}).
To better understand the power gain sensitivity to atmospheric conditions, we examine how the nondimensional parameters Fr and P_{N} impact the flow fields. The pressure gradients induced by inversion waves scale with g^{′}; therefore high inversion strengths correspond to strong inversionwave feedback and low Froude number values. These twodimensional waves are nondispersive with phase speed $\sqrt{{g}^{\prime}H}$ (Sutherland, 2010). Therefore, Fr also represents the ratio of the bulk wind speed within the ABL to the velocity of the inversion waves. If Fr<1 (subcritical flow) the twodimensional waves can affect the upstream flow, while they can travel only downstream if Fr>1 (supercritical flow). The flow is said to be critical when Fr=1. On the other hand, internalwaveinduced pressure gradients are governed by the second nondimensional group P_{N}. Strong internalwave feedback corresponds to low P_{N} values. In fact, strongly stratified atmospheres imply high N values, meaning that they account for higher internalwave oscillation frequencies and phase speed (Sutherland, 2010).
Two different flow regimes can be identified.

Regime 1 represents lowP_{N} flows. The strongly stratified free atmosphere limits vertical displacement of air parcels; hence reduced streamline divergence over the windfarm area is observed. This results in low velocity reductions and η_{t} values. Moreover, the flow fields are Fr independent in these atmospheric states (Smith, 2010).

Regime 2 represents highP_{N} flows. The inversionlayer strength determines the flow fields' properties since the influence of internal waves is negligible. The weakly stratified atmosphere makes the ABL behave like an idealized shallowwater system for Fr≃1 (choking effect; Smith, 2010). Moreover, the perturbations' magnitude is strongly dependent upon the Froude number.
Smith (2010) and Allaerts and Meyers (2019) defined a third regime where N=0 and ${g}^{\prime}=\mathrm{0}$, which would correspond to Fr, P_{N}→∞ or a purely neutral atmosphere. Gravity waves are not excited in this particular flow condition, and only drag forces and frictional effects play a role in the flow behaviour. Since we are interested in finding optimal thrust setpoint distributions which allow the recovery of gravitywaveinduced power loss, we did not investigate this regime in the current study.
Figure 6a–c illustrate the sensitivity of Fr to changes in h_{*}, ${\stackrel{\mathrm{\u203e}}{z}}_{\mathrm{0}}$ and $N/{f}_{\mathrm{c}}$ against the inversion parameter. In all cases, the Froude number ranges from approximatively 0.5 to 1.4. The black line denotes critical flow conditions. Lines of constant Froude number run parallel to this line, meaning that Fr is invariant and quasiinvariant to $N/{f}_{\mathrm{c}}$ and h_{*}, respectively. On the other hand, changes in ${\stackrel{\mathrm{\u203e}}{z}}_{\mathrm{0}}$ have a strong impact on the wind profile convexity and therefore on Fr. The sensitivity of P_{N} to the atmospheric state is displayed in Fig. 6d–f. P_{N} is not dependent on the inversion parameter. Hence, lines of constant P_{N} values are vertical and parallel to the dashed black line, which denotes atmospheric conditions for which P_{N}=1.5. This line divides the domain in regions where the internalwave effects are important (regime 1, right side) or limited (regime 2, left side). However, internal waves still play a crucial role in softening the flow perturbations' magnitude when P_{N} values are only slightly greater than 1, as in Fig. 6d and e. On the other hand, very high P_{N} numbers (P_{N}>10) are attained in weakly stratified conditions (see Fig. 6f). We will use the abovementioned regime classification as a proxy for the interpretation of the power gain sensitivity patterns (note that the terms high and low in the regime characterizations refer to the maximum and minimum Fr and P_{N} values found over the sensitivity domain).
Figure 7a–c and d–f illustrate the sensitivity of the optimal inversionlayer verticaldisplacement reduction 𝒢_{η} and power gain 𝒢 to changes in h_{*}, ${\stackrel{\mathrm{\u203e}}{z}}_{\mathrm{0}}$ and $N/{f}_{\mathrm{c}}$, against the inversion parameter. The displacement reduction is defined as
where ${\mathit{\eta}}_{\mathrm{t},\mathrm{max}}^{\mathrm{O}}$ and ${\mathit{\eta}}_{\mathrm{t},\mathrm{max}}^{\mathrm{R}}$ denote the maximum inversionlayer displacement attained with the optimal and reference model configurations, respectively. As we discussed in Sect. 4.1, the lowering of the inversionlayer vertical displacement reduces the strength of the adverse pressure gradient, increasing the flow wind speed and consequently the windfarm power output. Figure 7 confirms this statement. In fact, regions of high vertical displacement reduction strictly correspond to regions of high power gain.
Allaerts and Meyers (2017, 2018, 2019) found that for low ABL heights, gravity waves induce strong pressure gradients and play an important role in the distribution of the kinetic energy within the farm. Indeed, the large geostrophic wind angle found in shallow boundary layers redirects the favourable pressure gradient seen over the windfarm area of 90^{∘} for ${h}_{*}\to \mathrm{0}$, decreasing the dispersive impact of internal gravity waves. Figure 7d shows that the maximum power gain is indeed attained for ${h}_{*}=\mathrm{0.17}$ (i.e. for shallow boundary layer) in supercritical flow conditions, with gains of about 7.5 % corresponding to a displacement reduction of 17.5 %. A similar pattern is seen in Fig. 7e, where a maximum power gain of 8.4 % is attained again in supercritical conditions for ${\mathrm{log}}_{\mathrm{10}}\left({\stackrel{\mathrm{\u203e}}{z}}_{\mathrm{0}}\right)=\mathrm{4.2}$, corresponding to a displacement reduction of 19%. Both 𝒢 and 𝒢_{η} show higher sensitivity to changes in ${\stackrel{\mathrm{\u203e}}{z}}_{\mathrm{0}}$ than in h_{*}, decreasing rapidly for increasing value of surface roughness. Interestingly, power gains are close to zero in the case of high ${\stackrel{\mathrm{\u203e}}{z}}_{\mathrm{0}}$ values. This is due to the additional frictional drag which dissipates perturbation energy, limiting gravitywave excitation and consequently the potential of our optimization.
The sensitivity of 𝒢 and 𝒢_{η} to changes in freeatmosphere stability is shown in Fig. 7c and f. The high P_{N} sensitivity to changes in N (from P_{N}≈11 to P_{N}≈0.5 for increasing values of $N/{f}_{\mathrm{c}}$) accounts for a clear distinction between regime 1 and regime 2. The former shows power gains of about 5 % while the latter attains gains of 14 % corresponding to inversion displacement reductions of 24 %. Figure 7f illustrates that the power gain peak is obtained in critical flow conditions (Fr=1), differently from the previous cases. The very high P_{N} values (hence, the limited presence of internal waves) attained when Fr=1 allow for the choking effect to take place (Smith, 2010; Allaerts and Meyers, 2019). Very large flow perturbations are triggered in these atmospheric conditions, leaving greater potential for power recovery. The choking effect is not visible in Fig. 7d and e since there P_{N}≈2 when Fr=1 (the flow perturbations are softened by internal waves).
Overall, higher inversionlayer displacement reductions and power gains are attained in critical and supercritical flow conditions for high P_{N} values (regime 2), that is for low h_{*}, ${\stackrel{\mathrm{\u203e}}{z}}_{\mathrm{0}}$ and $N/{f}_{\mathrm{c}}$. This is not surprising due to the strong impact that gravity waves have on a farm's performance in such conditions (see Sect. 4.1 or Smith, 2010, and Allaerts and Meyers, 2019). Moreover, we observe strong gradients of 𝒢 and 𝒢_{η} along contours of P_{N} in regime 2 and weak gradients in regime 1. This suggests that the flow properties are Fr independent for low P_{N} values, confirming the observations of Smith (2010).
In the current study, we for the first time investigated the potential of thrust setpoint optimization in large wind farms for mitigating gravitywaveinduced blockage effects, with the aim of increasing the windfarm energy extraction. Thus, a fast boundary layer model proposed by Allaerts and Meyers (2019) was adopted. The threelayer model simulates the atmospheric response to turbine drag in large wind farms by dividing the vertical structure of the atmosphere into three layers. This approach accurately captures the effects of regional pressure gradients induced by large wind farms at low computational expenses. We first added the time dependency to the model so that timeperiodic gravitywave patterns could be reproduced. Further, we reformulated the model as an optimization framework with the objective of maximizing the windfarm energy output at all costs. Gradient information was derived using the continuous adjoint method. To limit the computational cost, a boxfunction windfarm force model was used, which assumes that the force is distributed over the whole windfarm area. The windfarm layout was inspired by the works of Allaerts et al. (2018) and Allaerts and Meyers (2019), roughly representing the Belgian–Dutch offshore windfarm cluster.
The optimization model was applied to two different atmospheric states representative of subcritical (Fr=0.9) and supercritical (Fr=1.1) flow conditions. The optimal configurations were then compared with a reference model setup which uses a uniform thrustcoefficient distribution. We did not observe dynamic behaviour in the optimal thrust setpoint distributions for different choices of time horizon and time step, meaning that it is not necessary to excite nonstationary wave patterns to further increase the windfarm energy output. However, we observed interesting spatial patterns. The optimal thrust setpoint distributions turned out to be inversely related with the inversionlayer vertical displacement η_{t}. This has led to a sinusoidal and Ushaped ${C}_{\mathrm{T}}^{\mathrm{O}}$ distribution along the streamwise direction in sub and supercritical conditions, respectively. An inversionlayer displacement reduction of 14.5 % and 16.8 % was observed in sub and supercritical conditions, which lowered the adverse pressure gradient strength in the windfarm induction and entrance region. The reduced blockage effects allowed for higher flow wind speeds through the farm. The optimal configurations showed power gains of 5.3 % and 7 % in sub and supercritical conditions with respect to the reference model setup.
The atmospheric state is far from being constant in real case scenarios; therefore the power gain sensitivity to changes in atmospheric conditions was further studied. Thus, the developed thrust setpoint optimization tool was applied for several wind profiles, inversion strengths and atmosphere stratifications for a total of 1960 different atmospheric states. Regions of high inversionlayerdisplacement reduction in the sensitivity domain strictly corresponded to regions of high power gain. This has confirmed that it is essential to reduce the streamline divergence over the windfarm area for limiting gravitywaveinduced power loss. The strong gravitywave feedback in highP_{N} conditions made these atmospheric states the most suitable for power recovery purposes. Power gains of up to 14 % were found for weakly stratified atmospheres (P_{N}≈11) in correspondence with critical flow conditions (Fr=1). This is related to the large flow perturbations induced by the choking effect (Smith, 2010). Overall, power gains above 4 % were observed for 77 % of the cases. We note that the gravitywaveinduced power losses are also sensitive to the windfarm layout. Optimization of layout (including, e.g. relevant technoeconomical constraints) is however not considered here and can be an interesting topic for future research.
The results discussed in the current manuscript make windfarm setpoint optimization a promising tool for gravitywaveinduced power loss recovery. However, many challenges remain before this can be translated to real windfarm applications. In the current work, we did not include an explicit wake model in our model, and we have presumed that wake losses remain unchanged during optimization (i.e. η_{w} is assumed to be constant). In the future, an analytical windfarm wake model, such as, e.g. the one developed by Niayifar and PortéAgel (2016) and used by Allaerts and Meyers (2019), could be adopted for optimization. This would however also require better representation in the wake model of changing background variables and pressure gradients. For instance, gravitywaveinduced pressure gradient effects on turbine wake recovery could be included using the model proposed by Shamsoddin and PortéAgel (2018) that incorporates effects of pressure gradients. Furthermore, the use of a windfarm drag model which analytically computes the wake of each turbine would allow us to separately investigate the influence of wake effects and gravity waves on the optimal turbine set points. This is work for future research. Moreover, the threelayer model has been validated with LES results only (Allaerts and Meyers, 2019). We are planning to perform a more extensive validation of the model in the near future. Next, we also plan to apply the results obtained in this article to a higherfidelity model (i.e. our inhouse LES solver SPWind). However, this requires some work on the efficiency of nonreflecting boundary conditions in our LES solver (Allaerts and Meyers, 2017, 2018). Finally, we assumed that the free atmosphere is uniformly stratified and steady. The relaxation of these assumptions would extend the applicability of the model, e.g. to atmospheres with heightdependent Brunt–Väisälä frequency and geostrophic wind, among others.
The continuous adjoint method is briefly explained in Appendix A1. Next, the threelayer model adjoint equations and cost functional gradient are derived in Appendices A2 and A3, respectively. Finally, the comparison between a finitedifference approximation of the cost function gradient and the adjoint evaluation is performed in Appendix A4.
A1 Continuous adjoint method
We adopt the standard L^{2} inner product over the time interval [0, T] and simulation domain Ω:
where a and b are two generic vectors. Moreover, we denote with ψ=[u_{1}, v_{1}, u_{2}, v_{2}, p_{1}, p_{2}] the vector containing the state variables and with C_{T}=C_{T}(x, y, t) the control parameter.
The reduced cost functional is defined as
where
The gradient of the reduced cost functional $\mathrm{\nabla}\stackrel{\mathrm{\u0303}}{\mathcal{J}}$ is interpreted as the Riesz representation of the Gâteaux derivative operator at C_{T} in any arbitrary direction δC_{T}:
where ℋ denotes the control Hilbert space.
Next, we define the state constraints of the optimization problem (i.e. the threelayer model equations) with shorthand notation 𝒩(ψ, C_{T}). The reduced formulation of the optimization problem implies by definition that 𝒩(ψ(C_{T}), C_{T})=0; therefore we can write the reduced cost functional as
where ${\mathit{\psi}}^{*}=[{\mathit{\zeta}}_{\mathrm{1}}$, χ_{1}, ζ_{2}, χ_{2}, Π_{1}, Π_{2}] denotes the vector containing the adjoint variables which play the role of Lagrange multipliers. In fact, it is easy to notice that $\stackrel{\mathrm{\u0303}}{\mathcal{J}}\left({C}_{\mathrm{T}}\right)=\mathcal{L}({C}_{\mathrm{T}}$, ψ(C_{T}), ψ^{*}), where ℒ is the Lagrangian of the optimization problem in Eq. (16).
Using Eqs. (A4) and (A5), the gradient of the reduced cost functional can be expressed as
where $\mathit{\delta}\mathit{\psi}=\mathrm{d}\mathit{\psi}/\mathrm{d}{C}_{\mathrm{T}}\mathit{\delta}{C}_{\mathrm{T}}$. The adjoint of the operator $\partial \mathcal{N}/\partial \mathit{\psi}$ is given by
where the righthand side is found using integration by parts. Similarly, the adjoint of $\partial \mathcal{N}/\partial {C}_{\mathrm{T}}$ is expressed as
The boundary terms BT_{1} and BT_{2} arise as a result of the integration by parts. Due to spatial and time periodicity constraints, it is easy to show that ${\mathrm{BT}}_{\mathrm{1}}={\mathrm{BT}}_{\mathrm{2}}=\mathrm{0}$. Hence, substituting Eqs. (A7) and (A8) into Eq. (A6), we obtain
Further, we assume that the adjoint variables satisfy the following relation:
which defines the adjoint equations. Therefore, the adjoint gradient is given by
A2 Derivation of the adjoint equations
We apply relation Eq. (A7) for deriving the adjoint of the operator $\partial \mathcal{N}/\partial \mathit{\psi}$. Starting with the velocity perturbations in the windfarm layer, we have
and by computing an integration by parts we obtain
Note that the minus sign in the argument of ${\mathcal{F}}^{\mathrm{1}}\left(\widehat{\mathrm{\Phi}}\right)(\mathit{x}$, −t) does not come from classical integration by parts. In fact, given three functions f, g, h∈L^{1}(Ω), it can be shown that
where in the second passage we have changed the order of integration (Fubini's theorem). This property allows us to write
Similarly, for the velocity perturbations in the upper layer, we have that
Following the same procedure for the pressure perturbations p_{1} and p_{2}, we obtain
and
Using Eq. (A10), the resulting adjoint equations correspond to
The adjoint momentum equations of the upper layer are homogeneous since the adjoint windfarm drag force is felt only indirectly in this layer ($\partial \mathcal{K}/\partial {\mathit{u}}_{\mathrm{2}}=\mathrm{0}$). Moreover, $\partial \mathcal{K}/\partial {p}_{\mathrm{1}}=\partial \mathcal{K}/\partial {p}_{\mathrm{2}}=\mathrm{0}$. On the other hand, the adjoint momentum equations of the windfarm layer are driven by the cost function. Using Eq. (A3), we obtain
Figure A1 illustrates a planform view of the forward and adjoint solutions in subcritical flow conditions (Fr=0.9). Both solutions are derived assuming a steadystate formulation of the optimization problem. The numerical setup, windfarm layout and atmospheric state are the ones listed in Table 1. Due to integration by parts, the convective term is negative in the backward equations, causing the flow to propagate upstream (i.e. from right to left of our domain) as displayed in Fig. A1 (bottom row). Moreover, the wind farm acts as a source term, and it speeds up the adjoint solution instead of decelerating it, causing an acceleration within the windfarm area and in the wake region.
A3 Derivation of the gradient
The adjoint gradient of the cost function is derived using Eq. (A11). To compute the adjoint of the operator $\partial \mathcal{N}/\partial {C}_{\mathrm{T}}$, we need to evaluate the following inner product:
which is easily rewritten as
Moreover, we derive the first term on the righthand side of Eq. (A11) using Eq. (A3), which results in
Finally, we obtain the gradient expression by substituting Eqs. (A22) and (A23) in Eq. (A11), which gives
A4 Verification of the adjoint gradient
The aim of this paragraph is to assess the quality of the gradient through comparison with a finitedifference approximation. The comparison is done using a grid resolution of 500 m. All other parameters correspond to the ones listed in Table 1, with Fr=0.9.
We define with
the directional derivative of $\mathrm{\nabla}\stackrel{\mathrm{\u0303}}{\mathcal{J}}$ along δC_{T}, where $\mathrm{\nabla}\stackrel{\mathrm{\u0303}}{\mathcal{J}}$ is the gradient computed with Eq. (A24) and δC_{T} is a perturbation of the baseline control C_{T}. Using finite difference, the same directional derivative can be approximated as
The truncation error of Eq. (A26) is proportional to the order of magnitude of the step length α. Therefore, α should be as small as possible to limit the discretization error. However, small values of α induce roundoff errors due to finiteprecision floatingpoint arithmetic. In other words, relation Eq. (A26) provides accurate gradient information only for a lower and upperbounded range of step length values.
Next, we define
where R and ℰ represent the ratio and the relative error between gradient information computed with the adjoint and finitedifference method. If the continuous adjoint method provides correct gradient information, we expect R≃1 and ℰ to be sufficiently small.
The following generic baseline control is chosen:
where ${C}_{\mathrm{T}}^{\mathrm{Betz}}=\mathrm{8}/\mathrm{9}$, ${k}_{x}=\mathrm{2}\mathit{\pi}/{L}_{x}$ and ${k}_{y}=\mathrm{2}\mathit{\pi}/{L}_{y}$. Ideally, we should validate the adjointbased gradient against the finitedifference one for all possible perturbations δC_{T}. However, such validation would require us to solve the governing equations (forward and backward) 2.4×10^{3} times since the control space has such DOF using this numerical setup. This computation is too expensive; therefore we select a limited class of perturbations given by
for different values of a and b.
Results of the comparison are shown in Fig. A2. We can appreciate that for ${\mathrm{10}}^{\mathrm{11}}\le \mathit{\alpha}\le {\mathrm{10}}^{\mathrm{4}}$ the ratio R is very close to unity, and the relative error ℰ is on the order of 10^{−4}, showing the typical Ushaped curve (Nita et al., 2016). However, for smaller step length values the relative error increases due to the decreasing arithmetic accuracy of the finitedifferencebased gradient. The relative error also increases for $\mathit{\alpha}>{\mathrm{10}}^{\mathrm{4}}$ due to discretization errors. We can appreciate that Fig. A2b displays a firstorder truncation error in accordance with relation Eq. (A26).
In the current section, we use the optimization framework derived in Sect. 2.2 to find an optimal uniform and steady thrustcoefficient distribution that minimizes the gravitywaveinduced blockage effects. To avoid confusion, we will denote with ${C}_{\mathrm{T}}^{\mathrm{O}}$ and ${C}_{\mathrm{T}}^{\mathrm{O},\mathrm{u}}$ the optimal nonuniform and uniform distributions, respectively. The windfarm layout and the atmospheric state are the ones detailed in Sect. 3.
Figure B1a and b display the optimal spatially invariant ${C}_{\mathrm{T}}^{\mathrm{O},\mathrm{u}}$ together with the streamwise profile of ${C}_{\mathrm{T}}^{\mathrm{O}}$ through the centre of the farm and its averaged value over the windfarm area $\langle {C}_{\mathrm{T}}^{\mathrm{O}}\rangle $ for the sub and supercritical cases, respectively. Moreover, ${C}_{\mathrm{T}}^{\mathrm{R}}$ denotes the thrust distribution used in the reference model. Interestingly, ${C}_{\mathrm{T}}^{\mathrm{O},\mathrm{u}}$ corresponds to the average of the nonuniform distribution in both cases. Since ${C}_{\mathrm{T}}^{\mathrm{O}}$ is sensitive to the atmospheric conditions, we expect ${C}_{\mathrm{T}}^{\mathrm{O},\mathrm{u}}$ to also depend on the atmospheric state (in fact, we observe a different value of ${C}_{\mathrm{T}}^{\mathrm{O},\mathrm{u}}$ in sub and supercritical conditions).
In the current example, the power gain 𝒢 (see Eq. 21) over the reference model configuration obtained with the nonuniform distributions ${C}_{\mathrm{T}}^{\mathrm{O}}$ is 5.3 % and 7 % for the sub and supercritical cases, respectively. For the optimal uniform distributions, we obtain a power gain of 5 % and 6.6 %.
A grid sensitivity analysis is performed to determine the dependence of the optimization results on the grid cell size. To this end, we fix the size of the numerical domain to 1000H×400H, and we vary the grid resolution spanning from 5H to $H/\mathrm{3}$, or equivalently from 1.6×10^{4} to 3.6×10^{6} DOF per layer. The results obtained are compared with the ones derived on a finer grid with resolution equal to $H/\mathrm{4}$.
Figure C1a and b display the cost function and power gain relative error, respectively, which are computed as
where ${\mathcal{P}}^{\mathrm{F}}={\stackrel{\mathrm{\u0303}}{\mathcal{J}}}^{\mathrm{F}}/T$ and 𝒢^{F} are the cost function (scaled with the time horizon T) and power gain obtained with a $H/\mathrm{4}$ grid resolution while $\mathcal{P}=\stackrel{\mathrm{\u0303}}{\mathcal{J}}/T$ and 𝒢 are the ones obtained with coarser grids. Note that the optimal distributions are steady state; therefore ℰ_{𝒫} is not dependent on the choice of the time horizon T. The cost function is evaluated using the reference case setup. The power gain is obtained using the optimization model described in Sect. 2. The model setup is reported in Table 1.
Spectral methods are known to have exponential convergence when used for discretizing smooth functions (i.e. f∈C^{∞}). However, algebraic convergence is attained for functions f∈C^{p} with p≥0. Figure C1 illustrates that we obtain a firstorder convergence. This is due to the twodimensional Heaviside function B(x, y) used for representing the windfarm footprint, which is discontinuous with discontinuous derivatives. Figure C1 also confirms that the results of the optimization model are grid independent. In fact, the cost function and power gain values change about 1 % and 4 % when the number of grid cells is increased by 2 orders of magnitude (from 10^{4} to 10^{6}). This justifies the use of a coarser grid in the sensitivity study performed in Sect. 4.2.
The code used for the simulations is written in Python and can be provided by contacting the corresponding author.
The raw data of the simulation results can be provided by contacting the corresponding author.
LL and JM jointly set up the simulation studies in the current work. LL performed code implementations and carried out the simulations. LL and JM jointly wrote the manuscript.
The authors declare that they have no conflict of interest.
The authors acknowledge support from the Research Foundation Flanders (FWO, grant no. G0B1518N). The authors thank Nicole Van Lipzig for useful discussions. The computational resources and services in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation Flanders (FWO) and the Flemish Government department EWI.
This research has been supported by the Fonds Wetenschappelijk Onderzoek (grant no. G0B1518N).
This paper was edited by Katherine Dykes and reviewed by JanWillem van Wingerden, Alan Wai Hou Lio, Dries Allaerts, and one anonymous referee.
Allaerts, D. and Meyers, J.: Boundarylayer development and gravity waves in conventionally neutral wind farms, J. Fluid Mech., 814, 95–130, 2017. a, b, c, d, e
Allaerts, D. and Meyers, J.: Gravity Waves and WindFarm Efficiency in Neutral and Stable Conditions, Bound.Lay. Meteorol., 166, 269–299, 2018. a, b, c, d, e, f, g
Allaerts, D. and Meyers, J.: Sensitivity and feedback of windfarm induced gravity waves, J. Fluid Mech., 862, 990–1028, 2019. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z, aa
Allaerts, D., Vanden Broucke, S., Van Lipzig, N., and Meyers, J.: Annual impact of windfarm gravity waves on the BelgianDutch offshore windfarm cluster, J. Phys.: Conf. Ser., 1037, 072006, https://doi.org/10.1088/17426596/1037/7/072006, 2018. a, b, c, d
Baker, A. H., Jessup, E. R., and Manteuffel, T.: A technique for accelerating the convergence of restarted GMRES, SIAM J., 26, 962–984, 2005. a
Barthelmie, R., Pryor, S., Frandsen, S., Hansen, K., Schepers, J., Rados, K., Schlez, W., Neubert, A., Jensen, L., and Neckelmann, S.: Quantifying the impact of wind turbine wakes on power output at offshore wind farms, J. Atmos. Ocean. Tech., 27, 1302–1317, 2010. a
Bleeg, J., Purcell, M., Ruisi, R., and Traiger, E.: Wind Farm Blockage and the Consequences of Neglecting Its Impact on Energy Production, Energies, 11, 1609, https://doi.org/10.3390/en11061609, 2018. a, b
Bortolotti, P., Tarrés, H. C., Dykes, K., Merz, K., Sethuraman, L., Verelst, D., and Zahle, F.: IEA Wind Task 37 on Systems Engineering in Wind Energy – WP2.1 Reference Wind Turbines, Technical report, https://doi.org/10.2172/1529216, 2019. a
Burton, T., Sharpe, D., Jenkins, N., and Bossanyi, E.: Wind Energy Handbook, Wiley, New York, 2001. a
Byrd, R. H., Lu, P., Nocedal, J., and Zhu, C.: A limited memory algorithm for bound constrained optimization, SIAM J., 16, 1190–1208, 1995. a
Canuto, C., Hussaini, M. Y., Quarteroni, A., and Zang, T. A.: Spectral Methods in Fluid Dynamics, SpringerVerlag, Berlin, Germany, https://doi.org/10.1007/9783642841088, 1988. a
Csanady, G. T.: Equilibrium theory of the planetary boundary layer with an inversion lid, Bound.Lay. Meteorol., 6, 63–79, 1974. a
De Los Reyes, J. C.: Numerical PDEConstrained Optimization, in: Springer Briefs in Optimization, Springer Briefs in Optimization, Cham, Switzerland, https://doi.org/10.1007/9783319133959, 2015. a
Durran, D. R.: Mountain waves and downslope winds, Am. Meteorol. Soc., 23, 1990. a
Eliassen, A. and Palm, E.: On the transfer of energy in stationary mountain waves, Geophys. Norveg., 22, 3, 1960. 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, 2012. a, b
Frederik, J. A., Weber, R., Cacciola, S., Campagnolo, F., Croce, A., Bottasso, C., and Wingerden, J. W.: Periodic dynamic induction control of wind farms: proving the potential in simulations and wind tunnel experiments, Wind Energ. Sci., 5, 245–257, https://doi.org/10.5194/wes52452020, 2020. a, b
Gill, A. E.: AtmosphereOcean Dynamics, in: International Geophysics Series 30, Academic Press, San Diego, USA, 1982. a, b
Goit, J. P. and Meyers, J.: Optimal control of energy extraction in windfarm boundary layers, J. Fluid Mech., 768, 5–50, 2015. a, b
Kheirabadi, A. C. and Nagamune, R.: A quantitative review of wind farm control with the objective of wind farm power maximization, J. Wind Eng. Indust. Aerodynam., 192, 45–73, 2019. a
Munters, W. and Meyers, J.: Dynamic Strategies for Yaw and Induction Control of Wind Farms Based on LargeEddy Simulation and Optimization, Energies, 11, 177, https://doi.org/10.3390/en11010177, 2018. a, b
Nappo, C. J.: An Introduction to Atmospheric Gravity Waves, in: International Geophysics Series 85, Academic Press, Waltham, USA, 2002. a
Niayifar, A. and PortéAgel, F.: Analytical modeling of wind farms: A new approach for power prediction, Energies, 9, 741, https://doi.org/10.3390/en9090741, 2016. a, b
Nieuwstadt, F.: On the solution of the stationary, baroclinic Ekmanlayer equations with a finite boundarylayer height, Bound.Lay. Meteorol., 26, 377–390, 1983. a
Nita, C., Vandewalle, S., and Meyers, J.: On the efficiency of gradient based optimization algorithms for DNSbased optimal control in a turbulent channel flow, Comput. Fluids, 125, 11–24, 2016. a
Nocedal, J. and Wright, S. J.: Numerical Optimization, SpringerVerlag, New York, USA, 1999. a
Queney, P.: The problem of the airflow over mountains: a summary of theoretical studies, B. Am. Meteorol. Soc., 29, 16–26, 1948. a
Segalini, A. and Dahlberg, J. A.: Blockage effects in wind farms, Wind Energy, 23, 120–128, https://doi.org/10.1002/we.2413, 2019. a
Shamsoddin, S. and PortéAgel, F.: A model for the effect of pressure gradient on turbulent axisymmetric wakes, J. Fluid Mech., 837, R3, https://doi.org/10.1017/jfm.2017.864, 2018. a
Smith, R. B.: Linear theory of stratified hydrostatic flow past an isolated mountain, Tellus, 32, 348–364, 1980. a
Smith, R. B.: Interacting mountain waves and boundary layers, J. Atmos. Sci., 64, 594–607, 2007. a
Smith, R. B.: Gravity wave effects on wind farm efficiency, Wind Energy, 13, 449–458, 2010. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q
Smith, R. B., Jiang, Q., and Doyle, J. D.: A theory of gravity wave absorption by a boundary, J. Atmos. Sci., 63, 774–781, 2006. a
Stull, R. B.: An Introduction to Boundary Layer Meteorology, Kluwer Academic Publishers, Dordrecht, the Netherlands, 1988. a
Sutherland, B. R.: Internal gravity waves, Cambridge University Press, Cambridge, 2010. a, b
Volker, P. J. H.: Wake effects of large offshore wind farms – a study of the mesoscale atmosphere, PhD thesis, DTU Wind Energy, Roskilde, Denmark, 2014. a
Wolfe, P.: Convergence conditions for ascent methods, SIAM Rev., 11, 226–235, 1969. a
Wu, L. K. and PortéAgel, F.: Flow Adjustment Inside and Around Large FiniteSize Wind Farms, Energies, 10, 2164, https://doi.org/10.3390/en10122164, 2017. a, b
Wu, Y.T. and PortéAgel, F.: Simulation of turbulent flow inside and above wind farms: model validation and layout effects, Bound.Lay. Meteorol., 146, 181–205, 2013. a
 Abstract
 Introduction
 Methodology
 Numerical setup and case description
 Results and discussion
 Conclusions
 Appendix A: Derivation and verification of the adjoint equations and the adjoint gradient
 Appendix B: Optimal uniform thrust setpoint distribution
 Appendix C: Grid sensitivity
 Code availability
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Methodology
 Numerical setup and case description
 Results and discussion
 Conclusions
 Appendix A: Derivation and verification of the adjoint equations and the adjoint gradient
 Appendix B: Optimal uniform thrust setpoint distribution
 Appendix C: Grid sensitivity
 Code availability
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References