Set-point optimization in wind farms to mitigate effects of flow blockage induced by gravity waves

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 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 wind-farm induced gravity waves and blockage. The analysis is performed using a fast boundary layer model which divides the vertical structure of the 5 atmosphere into three layers. The wind-farm drag force is applied over the whole wind-farm area in the lowest layer and is directly proportional to the wind-farm thrust set-point distribution. We implement an optimization model in order to derive the thrust-coefficient distribution which maximizes the wind-farm energy extraction. We use a continuous adjoint method to efficiently compute gradients for the optimization algorithm, which is based on a quasi-Newton method. Energy gains are evaluated with respect to a reference thrust-coefficient distribution based on the Betz–Joukowsky set point. We consider thrust 10 coefficients that can change in space, as well as in time, i.e. considering time-periodic signals. However, in all our optimization results, we find that optimal thrust-coefficient distributions are steady; any time-periodic distribution is less optimal. The (steady) optimal thrust-coefficient 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 U-shaped curve when the flow is supercritical. The sensitivity of the energy gain to the atmospheric state is studied using the developed 15 optimization tool for almost two thousand different atmospheric states. Energy gains of up to 14% are found for weakly stratified atmospheres in critical flow regimes.


Methodology
We now introduce the approach used for modelling wind-farm induced gravity waves and the method applied for maximizing the wind-farm energy output. The three-layer model is described in Section 2.1 and its optimization formulation is derived in Section 2.2. 80 In the work of Smith (2010), the atmospheric response to wind-farm 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 few hundreds meters from 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 85 model, the three-layer model divides the ABL into two layers: the wind-farm 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 ), 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).

Three-layer model
The model equations are derived starting from the incompressible three-dimensional Reynolds-Averaged Navier-Stokes 90 (RANS) equations for the ABL (Stull, 1988). The vertical component of the velocity field is a lot smaller than the horizontal using the following equations for the two layers in the ABL ∂u 1 ∂t + U 1 · ∇u 1 + 1 ∂u 2 ∂t + U 2 · ∇u 2 + 1 ∂η 2 ∂t + U 2 · ∇η 2 + H 2 ∇ · u 2 = 0.
(4) 100 When the wind farm is not operating (i.e., the wind-farm drag force is zero), a horizontally invariant reference state 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 height-averaged 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 105 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 depth-averaged turbulent viscosity, f c = 2Ω sin φ is the Coriolis frequency, with Ω the angular velocity of the earth and φ the latitude, J = e x ⊗e y −e y ⊗e x is the two-dimensional rotation dyadic with e x and e y two-dimensional unit vectors in the x-and y-direction, respectively. Finally, the perturbation of the friction at the ground and at the interface between both layers are described by the matrices C and D .

110
The right-hand side of Eq. 1 is characterized by the wind-farm drag force f . To limit the computational cost, we use a boxfunction wind-farm force model similar to Smith (2010). Therefore, we do not include the wake effects of each wind turbine individually, but instead we assume that the drag force is distributed over the whole wind-farm area. As for the flow equations, the wind-farm 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 thrust-coefficient distribution. Hence, the drag force is given by and with B(x, y) a box function equal to one within the wind-farm area and zero outside. The wind-farm drag force magnitude in Eqs. 5 and 6 scales with 120 where s x and s y are the streamwise and spanwise turbine spacing relative to the rotor diameter, η w is the wake efficiency and γ = u 2 r / U 1 2 is a velocity shape factor with u r the rotor-averaged wind speed . Moreover, unit dyadic. Relation 6 is nonlinear since it contains a product between time and space-dependent variables (i.e, C T and u 1 ).
We decide to retain this term because it allows us to include gravity-wave feedback on wind-farm 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.
The total vertical displacement of the inversion layer η t = η 1 + η 2 triggers gravity waves which induce pressure perturbations p. The relation between these two quantities is given by (Gill, 1982) 130 where F −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 boundary-layer approximation ∂p/∂z = 0). The complex stratification coefficient Φ in Fourier components is expressed as Relation 9 is obtained from linear three-dimensional, non-rotating, non-hydrostatic gravity-wave 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 = g∆θ/θ 0 accounts for two-dimensional trapped lee waves (further named as inversion waves) which corrugate the capping inversion layer. The potential-temperature 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 9, where m denotes the vertical 140 wavenumber which is given by According to the sign of m 2 we can have propagating or evanescent waves. Moreover, Ω = ω − κ · U g denotes the intrinsic wave frequency with κ = (k, l) the horizontal wavenumber vector.
Finally, combining Eq. 8 with Eq. 2 and Eq. 4, we can write the continuity equations for the wind-farm and upper layer as where p = p 1 + p 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.

150
The goal of the optimization framework is to find a time-periodic optimal thrust-coefficient distribution C O T (x, y, t) that minimizes the gravity-wave induced blockage effects, maximizing the flow wind speed and consequently the wind-farm energy extraction over a selected time period T . The background atmospheric state itself is presumed to be steady. The wind-farm layout and the atmospheric state are inputs of the optimization model and are detailed in Section 3.
The objective function of the optimization model consists in the energy extracted by the farm in the time interval [0, T ], and 155 it is defined as where Ω = D x × D y is the computational domain area and C p (x, y, t) is the power coefficient distribution. By using axial momentum theory (Burton et al., 2001), we find that the power coefficient depends upon the thrust coefficient according to the following non-linear relationship so that the objective function J is non-convex.
The wind-farm 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 non-linear time-periodic PDE-constrained optimization problem 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 175 wind-farm and upper layer, which also define the unknowns of the three-layer model. The control parameters consist of the value of the thrust set point in each grid cell within the wind-farm area. Hence, the size of the control space is proportional to where N t represents the number of time steps within the time horizon T , while N wf x and N wf y denote the number of grid points within the wind-farm area along the x and y-direction, respectively.
It is common practice in PDE-constrained optimization problem to not optimize the cost functional J (ψ, 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, if we use the shorthand notation N ψ, C T for the state constraints, we are enforcing that N ψ(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.

190
The gradient of the reduced objective function ∇J is needed for the solution of the reduced optimization problem. To this end, we use a continuous adjoint method. The adjoint (or backward) equations are given by (see Appendix A2 for detailed The adjoint variables are grouped in the vector ψ * = ζ 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 perturba-200 tions. Using the solution of the adjoint equations, the gradient of the cost function is expressed as (see Appendix A3 for details) where dC p /dC T is computed from Eq. 14. To compute the gradient ∇J we need to solve the forward and backward equations.
Since the cost for solving the adjoint equation is roughly the same as for the forward equation, the computational cost for 205 evaluating ∇J is proportional to the cost of solving twice the state equations. A verification of the approach, by comparing the adjoint gradient to a standard finite-difference approximation, is presented in Appendix A4.

Numerical setup and case description
We define the model setup used to asses the potential of set-point optimization in mitigating self-induced gravity-wave effects in this section. We discuss the numerical setup in Section 3.1. Next, the selected wind-farm layout is presented in Section 3.2.

210
Finally, the atmospheric state is discussed in Section 3.3.

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 set-up. All terms in the equations are linear, besides the first-order term of the wind-farm drag force (and its adjoint). These 215 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 3/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 × D y = 1000 × 400 km 2 , so that the perturbations die out before being recycled. The grid has N x × N y = 4000 × 1600 grid points which corresponds to a space resolution of ∆x = ∆y = 250 220 m or 6.4 × 10 6 DOF per layer. Finally, different time horizon values are used spanning from T = 10 minutes to T = 10 hours with number of time steps ranging from N t = 12 to N t = 120. The discretized forward and backward equations form two systems of 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 L-BFGS-B (limited-memory Broyden-Fletcher-Goldfarb-Shanno with box constraint) algorithm (Byrd et al., 1995) is an iterative quasi-Newton method. In the current appli-225 cation, 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 iteratively the Newton equation, 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 right-hand side needs gradient information to be computed, which are provided by the continuous adjoint method (see Appendix A for derivation and valida-230 tion). Fig. 1(a) shows that the cost functional is converged after six L-BFGS-B 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. 1(b) shows that the cost functional is mainly minimized within the first TNC iteration and convergence is reached after only four iterations. However, here, 63 function evaluations are needed. Hence, we will use the L-BFGS-B algorithm for solving the PDE-constrained optimization problem in 235 the remainder of the article. To limit computational effort, a maximum of four L-BFGS-B iterations will be performed.

Wind-farm layout
We choose a wind-farm layout that resembles the Belgian-Dutch wind-farm offshore cluster located in the North Sea in size, but simplified to a rectangular shaped wind-farm with length L x = 20 km and width L y = 30 km. This choice also leads to a widthto-depth ratio of 3/2 which is known to be the ratio for which self-induced gravity-wave effects are maximum (Allaerts and 240 Meyers, 2019). The wind turbine relative spacing along the x-and y-direction are s x = s y = 5.61 (both non-dimensionalized with respect to the turbine rotor diameter D). The turbine relative spacings are needed for computing the coefficient β (see Eq. 7). Both the wake efficiency η w and γ are assumed to be constant within the wind-farm area and are set to 0.9, similar to . 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.

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 250 (2019), we select two atmospheric states for initial testing in Section 4.1, corresponding to a subcritical and supercritical condition (see further below).
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 K and 3.7 K is used, which will lead to a sub-and supercritical flow, respectively (see below). Finally, a free atmosphere lapse rate Γ = 1 K/km is chosen. The lapse rate also 255 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 z 0 = z 0 /H and the non-dimensional boundary-layer height h * = Hf c /u * are the input parameters of Nieuwstadt's model, where f c is the Coriolis frequency and H = H 1 + H 2 is the ABL height. The wind-farm 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 fric-260 tion velocity is set to u * = 0.6 m/s. Finally, a surface roughness of z 0 = 10 −1 m is adopted. Using h * = 0.166 and z 0 = 10 −4 Table 1. Numerical setup, wind-farm layout and atmospheric state used in this manuscript.

Domain size
Dx × Dy = 1000 × 400 km 2 Grid size Nx × Ny = 4000 × 1600 The pressure gradient strengths induced by inversion and internal gravity waves are dependent upon the Froude number (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 (F r < 1), and a Froude number of 1.1 for ∆θ = 3.7 K, which leads to supercritical flow conditions (F r > 1). Further, P N = 1.92.
The numerical setup, wind-farm layout and background state variables are summarized in Table 1.

Results and discussion
We discuss the results of the optimization problem in the current section. Firstly, in Section 4.1 the optimal thrust-coefficient 275 distributions and relative energy gains are illustrated for two specific flow conditions introduced in Section 3.3. Thereafter, the sensitivity of the energy gain to the atmospheric state is studied in Section 4.2.

Optimal thrust-coefficient distributions
The optimization model described in Section 2.2 is time and space dependent. Hence, the model is capable of finding a timeperiodic optimal thrust-coefficient distribution over the wind-farm area in a fixed time interval [0, T ]. However, all optimal 280 thrust set-point 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 did not find any unsteady optimums. Therefore, we show only steady-state results in the remainder of the manuscript, and conclude that unsteady time-periodic excitation is less effective than a stationary spatially optimal distribution.
The steady-state optimal thrust-coefficient distributions obtained in sub-and supercritical conditions are analyzed in the 285 remainder of this section. To improve the understanding of such distributions, gravity-wave induced flow patterns obtained with C O T (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 set-point distribution over the wind-farm area is used, with C R T (x, y) = C Betz T = 8/9. in the wind-farm wake. The vertical displacement of air parcels triggers inversion waves on the 2D inversion layer surface and internal waves in the free atmosphere (3D waves). These waves induce pressure gradients, as visible in Fig. 2(b,e), where a region of high pressure builds up in correspondence with high η t values, leading to flow blockage. However, Fig. 2(b) shows a stronger adverse pressure gradient in the wind-farm induction region than the one in Fig. 2(e). In fact, inversion waves travel upstream in subcritical conditions, which leads to more slow-down in the induction region. In both sub-and supercritical 300 case, favourable pressure gradients develop within the wind-farm area which tend to accelerate the flow in the wind-farm exit region. Finally, Fig. 2(c,f) illustrate relative velocity reductions in the wind-farm layer. The stronger inversion strength found in the subcritical flow case transforms the inversion layer in a quasi-rigid lid, which limits vertical displacements. The lower streamlines divergence over the wind-farm area implies lower velocity reductions. Moreover, the favourable pressure gradient is stronger when F r = 0.9, allowing for higher flow wind speeds within the wind-farm area. This explains the higher velocity 305 reduction (up to 20%) seen in Fig. 2(f). Such a strong response could be on the limit of our small amplitude assumption. The planform view of pressure and velocity perturbations in the wind-farm 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 set-point distribution which reduces the velocity perturbations displayed in this end, we solve the optimization problem discussed in Section 2.2. The inputs of the optimization model are the wind-farm layout and the atmospheric conditions, which are detailed in Table 1. Moreover, an initial thrust-coefficient distribution needs to be specified. We have verified that for many different initial conditions the algorithm converges always to the same optimal solution. Therefore, a random initial thrust set-point distribution is chosen. The optimal configurations obtained for different Froude numbers are illustrated in Fig. 3. We find that the optimal thrust-coefficient distributions are non-uniform in space and 315 assume different spatial distributions according to the atmospheric state. In particular, when the flow is subcritical the optimal thrust set-point distribution assumes a sinusoidal behaviour in the streamwise direction while it becomes a U-shaped curve when the flow is supercritical. In both cases, C O T is almost invariant along the spanwise direction. We denote withJ R andJ O the energy extracted using C T = C R T = 8/9 and C T = C O T , respectively. Further, we define where G denotes the relative energy gain obtained using an optimal thrust-coefficient distribution instead of the reference one.
The energy gains attained in sub-and supercritical flow conditions are 5.3% and 7%, respectively. Clearly, energy gains are also strongly dependent on the atmospheric conditions. Therefore, a sensitivity study is carried out in Section 4.2.
The optimal set-point distributions displayed in Fig. 3 are related to the vertical displacement of the inversion layer over the wind-farm area. Fig. 4 shows streamwise profiles of η t and C O T through the center of the farm for F r = 0.9 and F r = 1.1. The 325 maximum value of η t is similar for both the sub-and supercritical case, although the profiles differ considerably. To reduce gravity-wave excitation, C O T is seen to be inversely related with η t . Moreover, Fig. 4(a) shows that the streamwise profile of η t has a sinusoidal behaviour. Hence, the optimal set-point distribution is sinusoidal as well, explaining the pattern displayed in Fig. 3(a). On the other hand, η t assumes a U-shaped profile through the wind farm in supercritical conditions (see Fig. 4(b)), a profile that is also found in C O T (see Fig. 3(b)). Moreover, Fig. 2(a,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 set-point distributions along the y-direction. Fig. 4 also shows that η R t,max > η O t,max in both sub-and supercritical conditions, meaning that the optimal thrust set-point distribution decreases the upward flow displacement over the wind-farm area. The maximum inversion-layer displacement is located at the entrance region of the farm. If we compare η R t and η O t in this region, a displacement reduction of 14.5% and 16.8% is attained with the optimal configuration for the sub-and supercritical case, respectively.

335
A lower vertical displacement of the inversion layer reduces gravity-wave 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 R T . Fig. 5(a,c) confirm this hypothesis, showing streamwise profiles of pressure perturbations p R and p O through the center of the farm for F r = 0.9 and F r = 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 case, respectively. Fig. 5(b,d) show streamwise profiles 340 of velocity perturbations u R 1 and u O 1 through the center of the farm for F r = 0.9 and F r = 1.1. The lower adverse pressure gradient strength in 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 case. Since the energy varies with the velocity cube, the higher velocity gain obtained in supercritical conditions 345 with respect to the subcritical ones explains the higher energy gain attained.
The optimal thrust-coefficient distributions and energy 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 energy gain to the atmospheric state is performed in the next section. (d) F r = 1.1 Figure 5. Streamwise profiles of (a,c) reference (p R ) and optimal (p O ) pressure perturbation and (b,d) reference (u R 1 ) and optimal (u O 1 ) velocity perturbation in subcritical (top row, Fr = 0.9) and supercritical (bottom row, Fr = 1.1) flow conditions. The wind-farm region is marked by vertical dashed lines, and the profiles have been obtained through the centre of the farm (y = 0).

Sensitivity study
Allaerts and Meyers (2019) pointed out that gravity-wave induced power loss is significant only for certain atmospheric states.
Since our aim is to recover its power loss, we also expect the energy gain to be sensitive to the atmospheric conditions. We note that gravity-wave patterns are also sensitive to the wind-farm layout. However, a sensitivity study over the wind-farm layout is beyond the scope of the article.

355
The nondimensionalization of the three-layer model equations with respect to the boundary layer height H and the friction velocity u * highlights four non-dimensional groups that govern the atmospheric state, which are: -The non-dimensional boundary layer height h * = Hf c /u * . Values of h * ≈ 0.1 denote shallow boundary layers typically found over sea, while h * ≈ 0.35 rather relates to a deep land-based boundary layer. We vary h * between 0.16 and 0.4; -The non-dimensional surface roughness length z 0 = z 0 /H. This number varies on several order of magnitude according 360 to the sea state or land surface. We vary log 10 (z 0 ) between −4.2 and −2.8 in the current study; wave theory which expresses the highest possible frequency for internal gravity waves (Gill, 1982). Typical values of free atmosphere lapse rate Γ range between 1 and 10 K/km. Low and high Γ values are associated with weakly and strongly stratified atmospheres, respectively. We vary Γ between 0.03 and 12 K/km corresponding to 10 ≤ N/f c ≤ 200; 365 -The inversion parameter g H/Au 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 H/Au 2 * ≈ 1, with A = 500 being an empirical constant. We vary the inversion parameter between 0.5 and 1.5. (2019) conducted a similar sensitivity study on the gravity-wave induced power loss on a wider range of non-dimensional numbers. However, since we are optimizing turbine thrust set points, we need to ensure that U 1 /U r < 1 370 (U r = 11 m/s is the rated wind speed of the DTU 10 MW IEA wind turbine), otherwise turbines would operate in above-rated wind speed regime and it would not make any sense to optimize their power production. The choice of the four ranges for the non-dimensional groups discussed above ensures that U 1 /U r ≤ 0.9 for all atmospheric states.

Allaerts and Meyers
Using the atmospheric state reported in Table 1, the non-dimensional numbers assume values h * = 0.166, z 0 = 10 −4 and N/f c = 58. The inversion parameter is equal to 1.046 and 0.691 in the sub-and supercritical case, respectively. The optimal 375 thrust-coefficient distributions discussed in Section 4.1 were obtained using these dimensionless group values. The sensitivity of the energy gain to atmospheric conditions is performed by varying h * , z 0 and N/f c against the inversion parameter g H/Au 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 four times bigger (∆x × ∆y = 1000 × 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 380 sensitivity study in Appendix B showing that the energy gain value changes of about 1% when the number of grid cells is increased of one order of magnitude (see Fig. B1). The high computational efficiency of the three-layer 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 condition). Since the wind-farm layout impact on energy gains is beyond the scope of our study, we impose the wind direction to be along the x-axis in the wind-farm layer in all simulations (V 1 = 0 m/s).

385
To better understand the energy gain sensitivity to atmospheric conditions, we examine how the non-dimensional parameters F r and P N impact the flow fields. The pressure gradients induced by inversion waves scale with g , therefore high inversion strengths correspond to strong inversion-wave feedback and low Froude number values. These two-dimensional waves are nondispervive with phase speed √ g H (Sutherland, 2010). Therefore, F r also represents the ratio of the bulk wind speed within the ABL to the velocity of the inversion waves. If F r < 1 (subcritical flow) the two-dimensional waves can affect the 390 upstream flow, while they can travel only downstream if F r > 1 (supercritical flow). The flow is said to be critical when F r = 1.
On the other hand, internal-wave induced pressure gradients are governed by the second non-dimensional group P N . Strong internal-wave feedback correspond to low P N values. In fact, strongly stratified atmospheres imply high N values, meaning that they account for higher internal-wave oscillation frequencies and phase speed (Sutherland, 2010).
Two different flow regimes can be identified: -Regime 1: low P N . The strongly stratified free atmosphere limits vertical displacement of air parcels, hence reduced streamline divergence over the wind-farm area is observed. This results in low velocity reductions and η t values. Moreover, the flow fields are F r -independent in these atmospheric states (Smith, 2010).
-Regime 2: high P N . The inversion-layer strength determines the flow fields properties since the influence of internal waves is negligible. The weakly stratified atmosphere makes the ABL to behave like an idealized shallow-water system 400 for F r 1 (choking effect (Smith, 2010)). Moreover, the perturbations magnitude are strongly dependent upon the Froude number. Smith (2010) and Allaerts and Meyers (2019) defined a third regime where N = 0 and g = 0, which would correspond to F r , P N → ∞ or to 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 set-point 405 distributions which allow to recover gravity-wave induced power loss, we did not investigate this regime in the current study. Fig. 6(a-c) illustrate the sensitivity of F r to changes in h * , z 0 and N/f 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 F r is invariant and quasi-invariant to N/f c and h * , respectively. On the other hand, changes in z 0 have a strong impact on the wind profile convexity and therefore on F r . The sensitivity of P N to 410 the atmospheric state is displayed in Fig. 6(d-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 internal-wave 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. 6(d,e). On the other hand, very high P N numbers (P N > 10) are attained in weakly stratified 415 conditions (see Fig. 6(f)). We will use the above mentioned regimes classification as a proxy for the interpretation of the energy gain sensitivity patterns (note that the term high and low in the regimes characterization are referred to the maximum and minimum F r and P N values found over the sensitivity domain). Fig. 7(a-c) and Fig. 7(d-f) illustrate the sensitivity of the optimal inversion-layer vertical-displacement reduction G η , and energy gain G to changes in h * , z 0 and N/f c , against the inversion parameter. The displacement reduction is defined as where η O t,max and η R t,max denote the maximum inversion-layer displacement attained with the optimal and reference model configuration, respectively. As we discussed in Section 4.1, the lowering of the inversion-layer vertical displacement reduces the strength of the adverse pressure gradient, increasing the flow wind speed and consequently the wind-farm energy output.  Allaerts and Meyers (2017 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 wind-farm area of 90 degrees for h * → 0, decreasing the dispersive impact of internal gravity waves. Fig. 7(d) displays that the maximum energy gain is indeed attained 430 for h * = 0.17 (i.e., for shallow boundary layer) in supercritical flow conditions, with gains of about 7.5% in correspondence to a displacement reduction of 17.5%. A similar pattern is seen in Fig. 7(e), where a maximum energy gain of 8.4% is attained again in supercritical conditions for log 10 (z 0 ) = −4.2, in correspondence to a displacement reduction of 19%. Both G and G η show higher sensitivity to changes in z 0 than in h * , decreasing rapidly for increasing value of surface roughness. Interestingly, energy gains are close to zero in case of high z 0 values. This is due to the additional frictional drag which dissipates perturbation 435 energy, limiting gravity-wave excitation and consequently the potential of our optimization.
The sensitivity of G and G η to changes in free atmosphere stability are shown in Fig. 7(c,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 c ) accounts for a clear distinction between regime 1 and regime 2. The former shows energy gains of about 5% while the latter attains gains of 14% in correspondence to inversion displacement reductions of 24%. Fig. 7(f)   differently from the previous cases. The very high P N values (hence, the limited presence of internal waves) attained in correspondence to F r = 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 energy recovery. The choking effect is not present in Fig. 7(d,e), since there P N ≈ 2 when F r = 1 (the flow perturbations are softened by internal waves).
Overall, higher inversion-layer displacement reductions and energy gains are attained in critical and supercritical flow con-445 ditions for high P N values (regime 2), that is for low h * , z 0 and N/f c . This is not surprising due to the strong impact that gravity waves have on farm's performance in such conditions (see Section 4.1 or Smith (2010) and Allaerts and Meyers (2019)).
Moreover, we observe strong gradients of G and G η along contours of P N in regime 2, and weak gradients in regime 1. This suggests that the flow properties are F r -independent for low P N values, confirming the observations of Smith (2010).
In the current study, we investigated for the first time the potential of thrust set-point optimization in large wind farms for mitigating gravity-wave induced blockage effects, with the aim of increasing the wind-farm energy extraction. Thus, a fast boundary layer model proposed by Allaerts and Meyers (2019) was adopted. The three-layer 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 455 first added the time-dependency to the model so that time-periodic gravity-wave patterns could be reproduced. Further, we reformulated the model as an optimization framework with the objective of maximizing the wind-farm energy output at all costs. Gradient information was derived using a continuous adjoint method. To limit the computational cost, a simple windfarm force model was used which assumes that the force is distributed over the whole wind-farm area. Hence, turbine-wake effects were not analytically resolved. The wind-farm layout was inspired by the works of  and Allaerts 460 and Meyers (2019), roughly representing the Belgian-Dutch offshore wind-farm cluster.
The optimization model was applied to two different atmospheric states representative of subcritical (F r = 0.9) and supercritical (F r = 1.1) flow conditions. The optimal configurations were then compared with a reference model setup which uses a uniform thrust-coefficient distribution. We did not observe dynamic behaviour in the optimal thrust set-point distributions for different choices of time horizon and time step, meaning that it is not necessary to excite non-stationary wave patterns to 465 further increase the wind-farm energy output. However, we observed interesting spatial patterns. The optimal thrust set-point distributions turned out to be inversely related with the inversion-layer vertical displacement η t . This has led to a sinusoidal and U-shaped C O T distribution along the streamwise direction in sub-and supercritical conditions, respectively. An inversion-layer displacement reduction of 14.5% and 16.8% was observed in sub-and supercritical conditions, which lowered the adverse pressure gradient strength in the wind-farm induction and entrance region. The reduced blockage effects allowed for higher 470 flow wind speeds through the farm. The optimal configurations showed energy 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 energy gain sensitivity to changes in atmospheric conditions was further studied. Thus, the developed thrust set-point optimization tool was applied for several wind profiles, inversion strengths and atmosphere stratifications for a total of 1960 different atmospheric states. Regions of high 475 inversion-layer-displacement reduction in the sensitivity domain strictly corresponded to regions of high energy gain. This has confirmed that it is essential to reduce the streamline divergence over the wind-farm area for limiting gravity-wave induced power loss. The strong gravity-wave feedback in high P N conditions made these atmospheric states the most suitable for energy recovery purposes. Energy gains up to 14% were found for weakly stratified atmospheres (P N ≈ 11) in correspondence of critical flow conditions (F r = 1). This is related to the large flow perturbations induced by the chocking effect (Smith, 2010). 480 Overall, energy gains above 4% were observed for 77% of the cases.
The results discussed in the current manuscript make wind-farm set-point optimization a promising tool for gravity-wave induced power loss recovery. However, many challenges remain before this can be translated to real wind-farm applications. First of all, the wind-farm force model needs to be improved. The analytical wake model developed by Niayifar and Porté-Agel (2016) and used by Allaerts and Meyers (2019) could be adopted for including turbine-wake effects in the optimization frame-485 work. Further, gravity-wave induced pressure gradient effects on turbine wakes recovery could be included using the model proposed by Shamsoddin and Porté-Agel (2018). This would allow us to optimize the turbine thrust set point of individual wind turbine placed in large wind farms. However, the complexity of the adjoint equations and the computational time will increase considerably. Moreover, the three-layer model has been validate with LES results only (Allaerts and Meyers (2019)). A more extensive validation of the model is topic for further research. Finally, we assumed that the free atmosphere is uniformly 490 stratified, and steady. The relaxation of these assumptions would extend the applicability of the model, e.g, to atmospheres with height-dependent Brunt-Väisälä frequency, and geostrophic wind, among others.

Appendix A: Derivation and verification of the adjoint equations and the adjoint gradient
The continuous adjoint method is briefly explained in Appendix A1. Next, the three-layer model adjoint equations and cost functional gradient are derived in Appendix A2 and A3, respectively. Finally, the comparison between a finite difference 495 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 500 variables and with C T = C T (x, y, t) the control parameter.
The reduced cost functional is defined as The gradient of the reduced cost functional ∇J is interpreted as the Riesz representation of the Gâteaux derivative operator at C T in any arbitrary direction δC T where H denotes the control Hilbert space.
notation N ψ, C T . The reduced formulation of the optimization problem implies by definition that N ψ(C T ), C T = 0, therefore we can write the reduced cost functional as where ψ * = ζ 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 thatJ C T = L C T , ψ(C T ), ψ * , where L is the Lagrangian of the optimization problem in 515 Eq. 15.
Using A4 and A5, the gradient of the reduced cost functional can be expressed as where δψ = dψ/dC T δC T . The adjoint of the operator ∂N /∂ψ is given by where the right-hand side is found using integration by parts. Similarly, the adjoint of ∂N /∂C 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 BT 1 = BT 2 = 0. Hence, substituting A7 and A8 into 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 530 We apply relation A7 for deriving the adjoint of the operator ∂N /∂ψ. Starting with the velocity perturbations in the wind-farm layer, we have Following the same procedure for the pressure perturbations p 1 and p 2 , we obtain Using A10, the resulting adjoint equations correspond to (A20)

550
The adjoint momentum equations of the upper layer are homogeneous, since the adjoint wind-farm drag force is felt only indirectly in this layer (∂K/∂u 2 = 0). Moreover, also ∂K/∂p 1 = ∂K/∂p 2 = 0. The adjoint momentum equations of the windfarm layer are driven by the cost function. Using A3, we obtain   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 wind-farm area and in the wake region.

A3 Derivation of the gradient
The adjoint gradient of the cost function is derived using relation A11. To compute the adjoint of the operator ∂N /∂C T , we need to evaluate the following inner product Moreover, we derive the first term on the right-hand side of A11 using A3, which results in Finally, we obtain the gradient expression by substituting A24 and A25 in A11, which gives 570 ∇J = βB(x, y) 1

A4 Verification of the adjoint gradient
The aim of this paragraph is to asses the quality of the gradient trough comparison with a finite difference approximation.
The comparison is done using a grid resolution of 500 m. All other parameters correspond to the ones listed in Table 1, with F r = 0.9.

575
We define with ∇J ADJ = ∇J , δC T (A27) the directional derivative of ∇J along δC T , where ∇J is the gradient computed with A26 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 A28 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 round-off errors due to finite-precision floatingpoint arithmetic. In other words, relation A28 provides accurate gradient information only for a lower an upper bounded 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 finite difference 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 C B T (x, y) = C Betz T 1 2 + 1 5 cos k x x + π + 1 5 sin k y y + π/5 (A31) where C Betz T = 8/9, k x = 2π/L x and k y = 2π/L y . Ideally, we should validate the adjoint-based gradient against the finitedifference one for all possible perturbations δC T . However, such validation would require 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 595 too expensive, therefore we select a limited class of perturbations given by δC T (x, y) = cos ak x x + π + sin bk y y + π/5 (A32) for different values of a and b.
Results of the comparison are shown in Fig. A2. We can appreciate that for 10 −11 ≤ α ≤ 10 −4 the ratio R is very close to unity and the relative error is in the order of 10 −4 , showing the typical U-shaped curve (Nita et al., 2016). However, for 600 smaller step length values the relative error increases due to the decreasing arithmetic accuracy of the finite-difference based gradient. The relative error also increases for α > 10 −4 due to discretization error. We can appreciate that Fig. A2b displays a first-order truncation error in accordance with relation A28.

Appendix B: Grid sensitivity
A grid sensitivity analysis is performed to determine the dependence of the optimization results on the grid cell size. To this 605 end, we fix the size of the numerical domain to 1000H × 400H and we vary the grid resolution spanning from 5H to H/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/4. Grid cells whereJ F and G F are the cost function and energy gain obtained with a H/4 grid resolution whileJ and G are the ones obtained with coarser grids. The cost function is evaluated using the reference case setup. The energy gain is obtained using the optimization model described in Section 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 ∞ ).

615
However, algebraic convergence is attained for functions f ∈ C p with p ≥ 0. Fig. B1 illustrates that we obtain a first-order convergence. This is due to the two-dimensional Heaviside function B(x, y) used for representing the wind-farm footprint, which is discontinuous with discontinuous derivatives. Fig. B1 also confirms that the results of the optimization model are grid-independent. In fact, the cost function and energy gain values change of about 1% and 4% when the number of grid cells is increased by two orders of magnitude (from 10 4 to 10 6 ). This justifies the use of a coarser grid in the sensitivity study 620 performed in section 4.2. Competing interests. The authors declare that they have no conflict of interest.