Set-point optimization in wind farms to mitigate effects of flow blockage induced by atmospheric 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 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 wind-farm-induced 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 wind-farm drag force is applied over the whole wind-farm area in the lowest layer and is directly proportional to the windfarm 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. Power gains are evaluated with respect to a reference thrust-coefficient 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 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 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.


Introduction
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 wind-farm 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 wind-farm energy extraction 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 large-scale 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).
L. Lanzilao and J. Meyers: Set-point optimization in wind farms to mitigate effects of flow blockage 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 wind-farm inflow velocity . 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 wind-farm thrust-coefficient distribution that minimizes the gravitywave-induced blockage effects, maximizing the flow wind speed and therefore the power production. Moreover, we investigate the impact of different atmospheric conditions on the optimal thrust-coefficient 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 wind-farm-induced gravity waves only in recent years. In the pioneering work of Smith (2010), a quasi-analytical model of atmospheric response to windfarm drag was used for modelling gravity-wave excitation due to diverging streamlines above the wind-farm area. Results have shown that gravity-wave excitation is strongly dependent upon the height of the boundary layer and the stability of the atmosphere aloft. Later, a fast boundary-layer model was proposed by Allaerts and Meyers (2019), who highlighted the crucial role of the inversion layer in determining gravity-wave patterns. The authors also used this model for an annual energy production study of the Belgian-Dutch offshore wind-farm cluster, showing that the annual energy loss due to the effect of self-induced gravity waves might be on the order of 4 % to 6 % .
Gravity waves were also observed in mesoscale and largeeddy simulation (LES) models. Fitch et al. (2012) and Volker (2014) proposed two different wind-farm parameterizations for the Weather Research and Forecasting model (WRF). Wind-farm-induced gravity waves were observed in both cases, causing flow deceleration several kilometres upstream of the farm. 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 finite-size wind farm operating in a conventionally neutral boundary layer (CNBL) with different free-atmosphere stratification, and they conclude that strongly stratified atmospheres decrease the turbine power output up to 35 % with respect to the weakly stratified cases. Wind-farm 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 wind-farm 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 wind-farm-induced gravity waves, to improve power extraction in waked turbines. Important control mechanisms include wake redirection (by yawing and tilting of the turbine), and turbine de-rating strategies. Control actions that influence wind-farm physics on a much larger scale, such as self-induced gravity waves, have not been explored to date.
In the current work, we concentrate on using wind-farm control to alter/improve the interaction between the wind farm and its self-induced gravity-wave system. To this end, we use the fast boundary-layer model proposed by Allaerts and Meyers (2019), which divides the vertical structure of the atmosphere into three layers (from here on named the three-layer model in the paper), and we reformulate it as an optimization problem. The objective function is defined as the wind-farm energy extracted over a time period T , while the constraints are the model equations plus a box constraint for the wind-farm thrust set-point distribution C T (x, y, t). Note that we do not use the tip-speed ratio and/or the pitch angle distribution as control parameters. Instead, we directly control the thrust set-point 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 (2007Smith ( , 2010. Consequently, the optimal thrust-coefficient distribution computed using the optimization formulation of the three-layer model takes into account the effects of self-induced gravity waves. Hence, we investigate whether it is possible to mitigate gravity-wave-induced blockage effects by varying the thrust set-point distribution within the wind-farm area.
The remainder of this paper is formulated as follows. The three-layer model and its optimization formulation are introduced in Sect. 2. Next, Sect. 3 describes the numerical setup, wind-farm 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. L. Lanzilao and J. Meyers: Set-point optimization in wind farms to mitigate effects of flow blockage 249

Methodology
We now introduce the approach used for modelling windfarm-induced gravity waves and the method applied for maximizing the wind-farm energy output. The three-layer model is described in Sect. 2.1, and its optimization formulation is derived in Sect. 2.2.

Three-layer model
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 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 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).
The three-layer model has been validated against LES results by Allaerts and Meyers (2019) (see Sect. 3 VAL2) on a two-dimensional (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 wind-farm area with a MAE of 5.6 %. Note that the three-layer 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 non-linear 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 three-dimensional Reynolds-averaged 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 wind-farm 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 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 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; and 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 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 right-hand side of Eq. (1) is characterized by the windfarm drag force f . We use a box-function wind-farm 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 wind-farm 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 set-point distribution (i.e. the C T value in every grid cell within the farm). 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 f = f (0) + f (1) with

250
L. Lanzilao and J. Meyers: Set-point optimization in wind farms to mitigate effects of flow blockage 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 γ = u 2 r / U 1 2 is a velocity shape factor with u r the rotor-averaged wind speed . Moreover, I = e x ⊗ e x + e y ⊗ e y denotes the unit dyadic. Finally, C T (x, y, t) represents the thrust-coefficient distribution. To compute the thrust coefficient C T,k (t) of a turbine at location (x k , y k ), it is possible to evaluate the thrust set-point distribution C T (x k , y k , t). A more accurate connection between C T,k (t) and the drag force f would, for example require the use of an analytical wind-farm 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 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. 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 η t = η 1 + η 2 triggers gravity waves which induce pressure perturbations p. The relation between these two quantities is given by Smith (2010): 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 com-plex stratification coefficientˆ in Fourier components is expressed aŝ Relation (Eq. 10) 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 (from here on named 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 (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, = ω − κ · U g denotes the intrinsic wave frequency with κ = (k, l) the horizontal wavenumber vector.
(2) and (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.

Optimization model
The goal of the optimization framework is to find a timeperiodic 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 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 time-periodic control (i.e. leading to a moving time average of the optimal control that is steady and does not lead to end-of-time effects). The wind-farm layout and the atmospheric state are inputs of the optimization model and are detailed in Sect. 3. Note that the relation between overall wind-farm drag and windfarm blockage is non-trivial. On the one hand, increased wind-farm drag leads to increased wind-farm blockage induced by gravity waves. This results from mass conservation and the upward displacement of the free atmosphere. On the other hand, increased wind-farm 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 non-linear 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 = D x × D y is the computational domain area. The non-linear relationship between C p and C T and the product between control and state variables in Eq. (15) imply 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 nonlinear time-periodic 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 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 N wf x N wf y N t , 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 directions, respectively.
It is common practice in a 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, defining an operator N (ψ, C T ) that denotes the state equations, we are enforcing 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 L. Lanzilao and J. Meyers: Set-point optimization in wind farms to mitigate effects of flow blockage where the only remaining constraints are the ones on the control parameters. The gradient of the reduced objective function ∇ 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 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 ψ * = [ζ 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 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 equations is roughly the same as for the forward equation, the computational cost for evaluating ∇ J is proportional to the cost of solving twice the state equations. To verify the approach, we compare the adjoint gradient to a standard finite-difference approximation in Appendix A4.

Numerical setup and case description
We define the model setup used to assess the potential of setpoint optimization in mitigating self-induced gravity-wave effects in this section. We discuss the numerical setup in Sect. 3.1. Next, the selected wind-farm layout is presented in Sect. 3.2. Finally, the atmospheric state is discussed in Sect. 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 setup. All terms in the equations are linear, except for the first-order term of the wind-farm 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 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 = 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 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 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 right-hand side needs gradient information to be computed, which is provided by the continuous adjoint method (see Appendix A for derivation and validation). ure 1 shows that the cost function decreases rapidly in the first two to three algorithm iterations, reaching a plateau afterwards. The use of a quasi-Newton 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 J (C T + αδC 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 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. 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 L-BFGS-B algorithm for solving the PDE-constrained optimization problem in the remainder of the article. To limit computational effort, a maximum of four L-BFGS-B 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 set-point 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.

Wind-farm layout
Allaerts and Meyers (2019) conducted a sensitivity study on the effects of wind-farm layout on gravity-wave-induced power losses. They show that these power losses monoton-ically increase with the size of the farm. Also, they mention that the losses are at a maximum when the wind-farm ratio L y /L x is close to 3/2, while being negligible for a very wide but short farm, and vice versa (i.e. negligible for L y /L x 1 and L x /L y 1). Since we are interested in optimal thrust-coefficient distributions in the presence of gravity waves, we have selected the "worst-case" wind-farm layout (i.e. a wind-farm 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  and Allaerts and Meyers (2019), which in size resembles the Belgian-Dutch wind-farm offshore cluster located in the North Sea, but simplified to a rectangular-shaped 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 = 5.61 (both non-dimensionalized 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 . 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.

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 free-atmosphere 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 non-dimensional surface roughness length z 0 = z 0 /H and the non-dimensional boundary-layer height h * = H f c /u * are the input parameters of Nieuwstadt 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 friction velocity is set to u * = 0.6 m s −1 . Finally, a surface roughness of z 0 = 10 −1 m is adopted. Using h * = 0.166 and z 0 = 10 −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 wind-farm 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 depth-averaged over the height H 1 and H 2 . Finally, the wind profile is oriented such that the wind in the wind-farm 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 F r = U B / g H and a non-dimensional group P N = U 2 B /N H U g , 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 (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 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, 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, the optimal thrust-coefficient 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.

Optimal thrust-coefficient distributions
The optimization model described in Sect. 2.2 is time and space dependent. Hence, the model is capable of finding a LGMRES  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 steady-state inflow conditions, meaning that we neglect meso-scale temporal variations in the velocity field (these could lead to time-dependent 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 steady-state results in the remainder of the manuscript, and we conclude for the time being that unsteady time-periodic 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 gravity-wave 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 large-scale 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 windfarm-induced gravity waves.
The steady-state optimal thrust-coefficient distributions obtained in sub-and supercritical conditions are analysed in the 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. Figure 2 illustrates a planform view of the perturbation flow patterns obtained with F r = 0.9 (top row) and F r = 1.1 (bottom row) using the reference model setup. The farm extracts energy from the flow, causing a momentum sink in the wind-farm layer. Due to the continuity constraint, an upward flow displacement above the wind-farm area takes place, which causes the boundary layer height to increase. Figure 2a shows that an inversion-layer vertical displacement of about 65 m takes place at the wind-farm 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 wind-farm centre. In both cases, the inversion-layer 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 inversion-layer 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 wind-farm induction region than the one in Fig. 2e. In fact, inversion waves travel upstream in subcritical conditions, which leads to more slow-down in the induction region. In both the suband 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 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 lower velocity deficits within the wind-farm 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 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 Fig. 2c and f. While maximizing the flow wind speed through the farm, we also maximize the wind-farm 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 thrust-coefficient 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 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 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 with P R = J R /T and P O = J O /T the power extracted using C T = C R T = 8/9 and C T = C O T , respectively. Further, we define where G denotes the relative power gain obtained using an optimal thrust-coefficient 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 non-uniform 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 set-point distributions displayed in Fig. 3 are related to the vertical displacement of the inversion layer over the wind-farm area. Figure 4 shows streamwise profiles of η t and C O T through the centre of the farm for F r = 0.9 and F r = 1.1. To reduce gravity-wave excitation, C O T 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 set-point distribution is sinusoidal as well, ex-  plaining the pattern displayed in Fig. 3a. On the other hand, η t assumes a U-shaped profile through the wind farm in supercritical conditions (see Fig. 4b), a profile that is also found in C O T (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 set-point distributions along the y direction. Figure 4 also shows that η O t,max < η R t,max in both suband 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 reduc-tion 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 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 . 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 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 cases, respectively. Figure 5b and d show streamwise profiles of velocity perturbations u R 1 and u O 1 through the centre of the farm for

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 (a, b, F r = 0.9) and supercritical (c, d, F r = 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). F r = 0.9 and F r = 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 F r = 1.1. The relative change in percentage between optimal and reference maximum flow perturbation values is summarized in Table 2.
The optimal thrust-coefficient 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.

Sensitivity study
Allaerts and Meyers (2019) pointed out that gravity-waveinduced 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 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. 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 as follows.
-The non-dimensional boundary layer height h * = H f 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 several orders of magnitude according to the sea state or land surface. We vary log 10 (z 0 ) between −4.2 and −2.8 in the current study.
-The non-dimensional Brunt-Väisälä frequency N/f c . The Brunt-Väisälä frequency is an important parameter in gravity-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 −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 10 ≤ N/f c ≤ 200.
-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.
Allaerts and Meyers (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 (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 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.
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 cases, respectively. The optimal thrust-coefficient 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 * , 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 4 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 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 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 state). Since the wind-farm 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 wind-farm layer in all simulations (V 1 = 0 m s −1 ).
To better understand the power 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 non-dispersive 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 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 corresponds 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 represents low-P N flows. 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 represents high-P N flows. The inversion-layer 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 shallow-water system for F r 1 (choking effect; Smith, 2010). Moreover, the perturbations' magnitude is strongly dependent upon the Froude number. (2019) defined a third regime where N = 0 and g = 0, which would correspond to F r, 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.

Smith (2010) and Allaerts and Meyers
Since we are interested in finding optimal thrust set-point distributions which allow the recovery of gravity-wave-induced power loss, we did not investigate this regime in the current study. Figure 6a-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 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 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. 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 above-mentioned 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 F r and P N values found over the sensitivity domain). Figure 7a-c and d-f illustrate the sensitivity of the optimal inversion-layer vertical-displacement reduction G η and power 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 configurations, respectively. As we discussed in Sect. 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 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 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 • for h * → 0, decreasing the dispersive impact of internal gravity waves. Figure 7d shows that the maximum power gain is indeed attained for h * = 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 log 10 (z 0 ) = −4.2, corresponding 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, power gains are close to zero in the case of high z 0 values. This is due to the additional frictional drag which dissipates perturbation energy, limiting gravity-wave excitation and consequently the potential of our optimization.
The sensitivity of G and G η to changes in free-atmosphere 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 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 (F r = 1), differently from the previous cases. The very high P N values (hence, the limited presence of internal waves) attained when 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 power recovery. The choking effect is not visible in Fig. 7d and 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 power gains are attained in critical and supercritical flow conditions 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 a farm's performance in such conditions (see Sect. 4.1 or Smith, 2010, andMeyers, 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).

Conclusions
In the current study, we for the first time investigated 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 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 the continuous adjoint method. To limit the computational cost, a box-function wind-farm force model was used, which assumes that the force is distributed over the whole wind-farm area. The wind-farm layout was inspired by the works of  and Allaerts 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 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 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 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 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 inversion-layer-displacement 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 gravity-wave-induced power loss. The strong gravity-wave feedback in high-P 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 (F r = 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 gravity-waveinduced power losses are also sensitive to the wind-farm 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 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 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 wind-farm 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, gravity-wave-induced 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 wind-farm 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 three-layer 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 higher-fidelity model (i.e. our in-house LES solver SP-Wind). However, this requires some work on the efficiency of non-reflecting boundary conditions in our LES solver 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.

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 Appendices A2 and A3, respectively. Finally, the comparison between a finite-difference 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 ∇ 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. Next, we define the state constraints of the optimization problem (i.e. the three-layer model equations) with shorthand 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 that J (C T ) = L(C T , ψ(C T ), ψ * ), where L 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 δψ = 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 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 ∂N /∂ψ. Starting with the velocity perturbations in the wind-farm layer, we have 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, ∂K/∂p 1 = ∂K/∂p 2 = 0. On the other hand, the adjoint momentum equations of the wind-farm layer are driven by the cost function. Using Eq. (A3), we obtain ∂K ∂u 1 = −3βC p B(x, y) U 1 U 1 . (A20) Figure A1 illustrates a planform view of the forward and adjoint solutions in subcritical flow conditions (F r = 0.9). Both solutions are derived assuming a steady-state 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 wind-farm 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 ∂N /∂C T , we need to evaluate the following inner product: which is easily rewritten as Moreover, we derive the first term on the right-hand 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 ∇ J = βB(x, y) H 1 U 1 U 1 · ζ 1 − H 1 U 1 dC p dC T U 1 2 + 3U 1 · u 1 + u 1 · U · ζ 1 . (A24)

A4 Verification of the adjoint gradient
The aim of this paragraph is to assess the quality of the gradient through 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. We define with ∇ J ADJ = ∇ J , δC T (A25) the directional derivative of ∇ J along δC T , where ∇ 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 round-off errors due to finiteprecision floating-point arithmetic. In other words, relation Eq. (A26) provides accurate gradient information only for a lower-and upper-bounded range of step length values. Next, we define where R and E 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 E 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 , 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 finite-difference 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 δC T (x, y) = cos (ak x x + π) + sin bk y y + π/5 (A30) 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 E is on the order of 10 −4 , showing the typical U-shaped curve (Nita et al., 2016). However, for 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 errors. We can appreciate that Fig. A2b displays a first-order truncation error in accordance with relation Eq. (A26).

Appendix C: Grid sensitivity
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 /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. Figure C1a and b display the cost function and power gain relative error, respectively, which are computed as where P F = J F /T and G F are the cost function (scaled with the time horizon T ) and power gain obtained with a H /4 grid resolution while P = J /T and G are the ones obtained with coarser grids. Note that the optimal distributions are steady state; therefore E P 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 first-order convergence. This is due to the two-dimensional 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.