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 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. 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.

1 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 non-local effects such as gravity waves may also have strong implications on the wind-farm energy extraction (Allaerts and Meyers, 2018, 2019).

In a stable atmosphere, an air parcel which is vertically perturbed will have the tendency to fall back to its original position. In such a case, an oscillation is initiated that is driven by gravity and inertia; this is called a gravity wave. Mountains are examples of orographic obstacles that trigger vertical flow displacement and consequently gravity waves (Smith, 1980). The drag force exerted by the mountain is usually transported upward by these waves. At the point of breakdown, the drag force is released in the upper levels of the atmosphere, causing a slowdown of the 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).

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 (Allaerts and Meyers, 2018). The size of this region scales with the length of the farm. This phenomenon is one possible cause of flow blockage. Note that it differs from classical hydrodynamic blockage caused by the turbine induction, which typically scales with the turbine rotor diameter and which has also been studied recently in much detail (Bleeg et al., 2018; Segalini and Dahlberg, 2019). The goal of the current study is to determine a wind-farm thrust-coefficient distribution that minimizes the gravity-wave-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 wind-farm 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 % (Allaerts et al., 2018).

Gravity waves were also observed in mesoscale and large-eddy 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. Allaerts and Meyers (2017, 2018) have investigated the interaction between an “infinitely” wide wind farm and both a conventionally neutral and stable boundary layer in typical offshore conditions in a LES framework. They found that for low ABL heights, gravity waves induce strong pressure gradients and play an important role in the distribution of the kinetic energy within the farm. Wu and Porté-Agel (2017) considered a large 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 (2007, 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 set-point 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.

2 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 Sect. 2.1, and its optimization formulation is derived in Sect. 2.2.

2.1 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 Allaerts et al., 2018) and a second layer up to the top of the ABL. Finally, the third layer models the free atmosphere above the ABL following the approach of Smith (2010).

The 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.

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{1}}}{\partial t}}& +{\mathit{U}}_{\mathrm{1}}\cdot \mathrm{\nabla }{\mathit{u}}_{\mathrm{1}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }p+f_{\mathrm{c}}\mathbf{J}\cdot {\mathit{u}}_{\mathrm{1}}-{\mathit{\nu }}_{\mathrm{t},\mathrm{1}}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{u}}_{\mathrm{1}}\\ & -{\displaystyle \frac{{\mathbf{D}}^{\prime }}{H_{\mathrm{1}}}}\cdot \left({\mathit{u}}_{\mathrm{2}}-{\mathit{u}}_{\mathrm{1}}\right)+{\displaystyle \frac{{\mathbf{C}}^{\prime }}{H_{\mathrm{1}}}}\cdot {\mathit{u}}_{\mathrm{1}}={\displaystyle \frac{\mathit{f}}{H_{\mathrm{1}}}}\end{array}\text{(2)}& {\displaystyle }{\displaystyle \frac{\partial {\mathit{\eta }}_{\mathrm{1}}}{\partial t}}+{\mathit{U}}_{\mathrm{1}}\cdot \mathrm{\nabla }{\mathit{\eta }}_{\mathrm{1}}+H_{\mathrm{1}}\mathrm{\nabla }\cdot {\mathit{u}}_{\mathrm{1}}=\mathrm{0}\text{(3)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{2}}}{\partial t}}& +{\mathit{U}}_{\mathrm{2}}\cdot \mathrm{\nabla }{\mathit{u}}_{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }p+f_{\mathrm{c}}\mathbf{J}\cdot {\mathit{u}}_{\mathrm{2}}-{\mathit{\nu }}_{\mathrm{t},\mathrm{2}}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{u}}_{\mathrm{2}}\\ & +{\displaystyle \frac{{\mathbf{D}}^{\prime }}{H_{\mathrm{2}}}}\cdot \left({\mathit{u}}_{\mathrm{2}}-{\mathit{u}}_{\mathrm{1}}\right)=\mathrm{0}\end{array}\text{(4)}& {\displaystyle }{\displaystyle \frac{\partial {\mathit{\eta }}_{\mathrm{2}}}{\partial t}}+{\mathit{U}}_{\mathrm{2}}\cdot \mathrm{\nabla }{\mathit{\eta }}_{\mathrm{2}}+H_{\mathrm{2}}\mathrm{\nabla }\cdot {\mathit{u}}_{\mathrm{2}}=\mathrm{0}\end{array}$$

When the wind farm is not operating (i.e. the wind-farm drag force is zero), a horizontally invariant reference state of (U₁, H₁) and (U₂, H₂) characterizes the wind farm and upper layer, where U₁=(U₁, V₁) and U₂=(U₂, V₂) are the height-averaged horizontal components of the background velocity and H₁, H₂ represent the height of the two layers. Whenever the farm extracts power from the flow, small velocity and height perturbations (u₁, η₁) and (u₂, η₂) 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, ρ₀ 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 $\mathbf{J}={\mathit{e}}_x\otimes {\mathit{e}}_y-{\mathit{e}}_y\otimes {\mathit{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 wind-farm 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 wind-farm 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 $\mathit{f}={\mathit{f}}^{\left(\mathrm{0}\right)}+{\mathit{f}}^{\left(\mathrm{1}\right)}$ with

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }{\mathit{f}}^{\left(\mathrm{0}\right)}=-\mathit{\beta }{C}_{\mathrm{T}}B(x,y)\parallel {\mathit{U}}_{\mathrm{1}}\parallel {\mathit{U}}_{\mathrm{1}},\text{(6)}& {\displaystyle }& {\displaystyle }{\mathit{f}}^{\left(\mathrm{1}\right)}=-\mathit{\beta }{C}_{\mathrm{T}}B(x,y){\mathbf{U}}^{\prime }\cdot {\mathit{u}}_{\mathrm{1}},\end{array}$$

where

$$\begin{array}{}\text{(7)}& {\mathbf{U}}^{\prime }={\displaystyle \frac{\mathrm{1}}{\parallel {\mathit{U}}_{\mathrm{1}}\parallel }}\left({\mathit{U}}_{\mathrm{1}}\otimes {\mathit{U}}_{\mathrm{1}}+\mathbf{I}\parallel {\mathit{U}}_{\mathrm{1}}{\parallel }^{\mathrm{2}}\right)\end{array}$$

and with B(x, y) a box function equal to 1 within the wind-farm area and zero outside. The x and y axes denote the streamwise and spanwise directions, respectively. The wind-farm drag force magnitude in Eqs. (5) and (6) scales with

$$\begin{array}{}\text{(8)}& \mathit{\beta }={\displaystyle \frac{\mathit{\pi }{\mathit{\eta }}_{\mathrm{w}}\mathit{\gamma }}{\mathrm{8}{s}_x{s}_y}},\end{array}$$

where s_x and s_y are the streamwise and spanwise turbine spacings relative to the rotor diameter, η_w is the wake efficiency and $\mathit{\gamma }=u_{\mathrm{r}}^{\mathrm{2}}/\parallel {\mathit{U}}_{\mathrm{1}}{\parallel }^{\mathrm{2}}$ is a velocity shape factor with u_r the rotor-averaged wind speed (Allaerts and Meyers, 2018). Moreover, $\mathbf{I}={\mathit{e}}_x\otimes {\mathit{e}}_x+{\mathit{e}}_y\otimes {\mathit{e}}_y$ denotes the unit dyadic. Finally, C_T(x, y, t) represents the thrust-coefficient distribution. To compute the thrust coefficient ${\stackrel{\mathrm{̃}}{C}}_{\mathrm{T},k}\left(t\right)$ 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 ${\stackrel{\mathrm{̃}}{C}}_{\mathrm{T},k}\left(t\right)$ 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 gravity-wave 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₁). We decide to retain this term because it allows us to include gravity-wave feedback on wind-farm energy extraction. In fact, f⁽¹⁾≥0 so that it reduces the drag force that the farm exerts on the flow, thereby reducing effects of blockage in the model. We note that Allaerts and Meyers (2019) have shown that the flow perturbations computed with this simple drag force model have similar trends and orders of magnitude as the ones computed using a drag model that relies on the more detailed analytical wake model of Niayifar and Porté-Agel (2016). Therefore, we believe that the model adopted is a reasonable representation of reality.

The total vertical displacement of the inversion layer ${\mathit{\eta }}_{\mathrm{t}}={\mathit{\eta }}_{\mathrm{1}}+{\mathit{\eta }}_{\mathrm{2}}$ triggers gravity waves which induce pressure perturbations p. The relation between these two quantities is given by Smith (2010):

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{p}{{\mathit{\rho }}_{\mathrm{0}}}}={\mathcal{F}}^{-\mathrm{1}}\left(\widehat{\mathrm{\Phi }}\right)\ast {\mathit{\eta }}_t,\end{array}$$

where ℱ⁻¹ 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 $\partial p/\partial z=\mathrm{0}$). The complex stratification coefficient $\widehat{\mathrm{\Phi }}$ in Fourier components is expressed as

$$\begin{array}{}\text{(10)}& \widehat{\mathrm{\Phi }}=g^{\prime }+{\displaystyle \frac{i\left(N^{\mathrm{2}}-{\mathrm{\Omega }}^{\mathrm{2}}\right)}{m}}.\end{array}$$

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^{\prime }=g\mathrm{\Delta }\mathit{\theta }/{\mathit{\theta }}_{\mathrm{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 θ₀ 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

$$\begin{array}{}\text{(11)}& m^{\mathrm{2}}=\left(k^{\mathrm{2}}+l^{\mathrm{2}}\right)\left({\displaystyle \frac{N^{\mathrm{2}}}{{\mathrm{\Omega }}^{\mathrm{2}}}}-\mathrm{1}\right).\end{array}$$

According to the sign of m², we can have propagating or evanescent waves. Moreover, $\mathrm{\Omega }=\mathit{\omega }-\mathit{\kappa }\cdot {\mathit{U}}_{\mathrm{g}}$ denotes the intrinsic wave frequency with κ=(k, l) the horizontal wavenumber vector.

Finally, combining Eq. (9) with Eqs. (2) and (4), we can write the continuity equations for the wind farm and upper layer as

$$\begin{array}{}\text{(12)}& {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\partial p_{\mathrm{1}}}{\partial t}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\mathit{U}}_{\mathrm{1}}\cdot \mathrm{\nabla }{p}_{\mathrm{1}}+H_{\mathrm{1}}\mathrm{\nabla }\cdot \left[{\mathcal{F}}^{-\mathrm{1}}\left(\widehat{\mathrm{\Phi }}\right)\ast {\mathit{u}}_{\mathrm{1}}\right]=\mathrm{0},\text{(13)}& {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\partial p_{\mathrm{2}}}{\partial t}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\mathit{U}}_{\mathrm{2}}\cdot \mathrm{\nabla }{p}_{\mathrm{2}}+H_{\mathrm{2}}\mathrm{\nabla }\cdot \left[{\mathcal{F}}^{-\mathrm{1}}\left(\widehat{\mathrm{\Phi }}\right)\ast {\mathit{u}}_{\mathrm{2}}\right]=\mathrm{0},\end{array}$$

where $p=p_{\mathrm{1}}+p_{\mathrm{2}}$ is intended as the sum of the pressure perturbations induced by the vertical displacements η₁ and η₂, respectively. This form will be used in the remainder of the manuscript.

2.2 Optimization model

The goal of the optimization framework is to find a time-periodic optimal thrust-coefficient distribution $C_{\mathrm{T}}^{\mathrm{O}}(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 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 wind-farm 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 wind-farm 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:

$$\begin{array}{}\text{(14)}& C_{\mathrm{p}}={\displaystyle \frac{C_{\mathrm{T}}}{\mathrm{2}}}\left(\mathrm{1}+\sqrt{\mathrm{1}-C_{\mathrm{T}}}\right).\end{array}$$

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

$$\begin{array}{}\text{(15)}& \begin{array}{rl}\mathcal{J}\left(\mathit{\psi },C_{\mathrm{T}}\right)=& -\mathit{\beta }\parallel {\mathit{U}}_{\mathrm{1}}\parallel \underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathrm{\Omega }}{\iint }{C}_{\mathrm{p}}B(x,y)\\ & \left(\parallel {\mathit{U}}_{\mathrm{1}}{\parallel }^{\mathrm{2}}+\mathrm{3}{\mathit{U}}_{\mathrm{1}}\cdot {\mathit{u}}_{\mathrm{1}}\right)\mathrm{d}\mathit{x}\mathrm{d}t,\end{array}\end{array}$$

where $\mathrm{\Omega }=D_x\times 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 𝒥 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 optimization problem constrained by a system of partial differential equations (PDEs):

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\min }_{\mathit{\psi },C_{\mathrm{T}}}& \mathcal{J}\left(\mathit{\psi },C_{\mathrm{T}}\right)\\ & \mathrm{s}.\mathrm{t}.\\ & {\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{1}}}{\partial t}}+{\mathit{U}}_{\mathrm{1}}\cdot \mathrm{\nabla }{\mathit{u}}_{\mathrm{1}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }{p}_{\mathrm{1}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }{p}_{\mathrm{2}}+f_{\mathrm{c}}\mathbf{J}\\ & \cdot {\mathit{u}}_{\mathrm{1}}-{\mathit{\nu }}_{\mathrm{t},\mathrm{1}}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{u}}_{\mathrm{1}}-{\displaystyle \frac{{\mathbf{D}}^{\prime }}{H_{\mathrm{1}}}}\cdot \left({\mathit{u}}_{\mathrm{2}}-{\mathit{u}}_{\mathrm{1}}\right)+{\displaystyle \frac{{\mathbf{C}}^{\prime }}{H_{\mathrm{1}}}}\cdot {\mathit{u}}_{\mathrm{1}}\\ & ={\displaystyle \frac{{\mathit{f}}^{\left(\mathrm{0}\right)}+{\mathit{f}}^{\left(\mathrm{1}\right)}}{H_{\mathrm{1}}}}\mathrm{in}\mathrm{\Omega }\times (\mathrm{0},T],\\ & {\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{2}}}{\partial t}}+{\mathit{U}}_{\mathrm{2}}\cdot \mathrm{\nabla }{\mathit{u}}_{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }{p}_{\mathrm{1}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }{p}_{\mathrm{2}}+f_{\mathrm{c}}\mathbf{J}\\ & \cdot {\mathit{u}}_{\mathrm{2}}-{\mathit{\nu }}_{\mathrm{t},\mathrm{2}}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{u}}_{\mathrm{2}}+{\displaystyle \frac{{\mathbf{D}}^{\prime }}{H_{\mathrm{2}}}}\\ & \cdot \left({\mathit{u}}_{\mathrm{2}}-{\mathit{u}}_{\mathrm{1}}\right)=\mathrm{0}\mathrm{in}\mathrm{\Omega }\times (\mathrm{0},T],\\ & {\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\partial p_{\mathrm{1}}}{\partial t}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\mathit{U}}_{\mathrm{1}}\cdot \mathrm{\nabla }{p}_{\mathrm{1}}+H_{\mathrm{1}}\mathrm{\nabla }\\ & \cdot \left[{\mathcal{F}}^{-\mathrm{1}}\left(\widehat{\mathrm{\Phi }}\right)\ast {\mathit{u}}_{\mathrm{1}}\right]=\mathrm{0}\mathrm{in}\mathrm{\Omega }\times (\mathrm{0},T],\\ & {\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\partial p_{\mathrm{2}}}{\partial t}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\mathit{U}}_{\mathrm{2}}\cdot \mathrm{\nabla }{p}_{\mathrm{2}}+H_{\mathrm{2}}\mathrm{\nabla }\\ & \cdot \left[{\mathcal{F}}^{-\mathrm{1}}\left(\widehat{\mathrm{\Phi }}\right)\ast {\mathit{u}}_{\mathrm{2}}\right]=\mathrm{0}\mathrm{in}\mathrm{\Omega }\times (\mathrm{0},T],\\ & \mathrm{0}\le C_{\mathrm{T}}<\mathrm{1}\mathrm{in}\mathrm{\Omega }\times (\mathrm{0},T],\\ & C_{\mathrm{T}}(x,y,\mathrm{0})=C_{\mathrm{T}}(x,y,T)\mathrm{in}\mathrm{\Omega }.\end{array}\end{array}$$

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₁, v₁, u₂, v₂, p₁, p₂] 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_x^{\mathrm{wf}}{N}_y^{\mathrm{wf}}{N}_t$, where N_t represents the number of time steps within the time horizon T, while $N_x^{\mathrm{wf}}$ and $N_y^{\mathrm{wf}}$ denote the number of grid points within the 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 𝒥(ψ, C_T) directly because such a problem would span both the state and control space. To avoid exploring the entire feasibility region, we require ψ(C_T) to be the solution of the state equations throughout the optimization process. In other words, defining an operator 𝒩(ψ, C_T) that denotes the state equations, we are enforcing 𝒩(ψ(C_T), C_T)=0 during optimization iterations. This technique leads us to a reduced optimization problem with feasibility region given by the control space (De Los Reyes, 2015). The reduced optimization problem is written as

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{\min }_{C_{\mathrm{T}}}\stackrel{\mathrm{̃}}{\mathcal{J}}\left(C_{\mathrm{T}}\right)& =\mathcal{J}\left(\mathit{\psi }\left(C_T\right),C_{\mathrm{T}}\right)\\ & \mathrm{s}.\mathrm{t}.\\ & \mathrm{0}\le C_{\mathrm{T}}<\mathrm{1}\mathrm{in}\mathrm{\Omega }\times (\mathrm{0},T],\\ & C_{\mathrm{T}}(x,y,\mathrm{0})=C_{\mathrm{T}}(x,y,T)\mathrm{in}\mathrm{\Omega },\end{array}\end{array}$$

where the only remaining constraints are the ones on the control parameters.

The gradient of the reduced objective function $\mathrm{\nabla }\stackrel{\mathrm{̃}}{\mathcal{J}}$ is needed for the solution of the reduced optimization problem. To this end, we use the continuous adjoint method. The adjoint (or backward) equations are given by (see Appendix A2 for detailed derivation)

$$\begin{array}{}\text{(18)}& \begin{array}{rl}-{\displaystyle \frac{\partial {\mathit{\zeta }}_{\mathrm{1}}}{\partial t}}& -{\mathit{U}}_{\mathrm{1}}\cdot \mathrm{\nabla }{\mathit{\zeta }}_{\mathrm{1}}+f_{\mathrm{c}}\mathbf{J}\cdot {\mathit{\zeta }}_{\mathrm{1}}-{\mathit{\nu }}_{\mathrm{t},\mathrm{1}}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\zeta }}_{\mathrm{1}}+{\displaystyle \frac{{\mathbf{D}}^{\prime }}{H_{\mathrm{1}}}}\cdot {\mathit{\zeta }}_{\mathrm{1}}\\ & +{\displaystyle \frac{{\mathbf{C}}^{\prime }}{H_{\mathrm{1}}}}\cdot {\mathit{\zeta }}_{\mathrm{1}}-{\displaystyle \frac{{\mathbf{D}}^{\prime }}{H_{\mathrm{2}}}}\cdot {\mathit{\zeta }}_{\mathrm{2}}-H_{\mathrm{1}}\left[{\mathcal{F}}^{-\mathrm{1}}\left(\widehat{\mathrm{\Phi }}\right)(-\mathit{x},-t)\right.\\ & \left.\ast \mathrm{\nabla }{\mathrm{\Pi }}_{\mathrm{1}}\right]+{\displaystyle \frac{\mathit{\beta }{C}_{\mathrm{T}}B(x,y)}{H_{\mathrm{1}}}}{\mathbf{U}}^{\prime }\cdot {\mathit{\zeta }}_{\mathrm{1}}=\mathrm{3}\mathit{\beta }{C}_{\mathrm{p}}B(x,y)\\ & \parallel {\mathit{U}}_{\mathrm{1}}\parallel {\mathit{U}}_{\mathrm{1}}\mathrm{in}\mathrm{\Omega }\times (\mathrm{0},T],\\ & -{\displaystyle \frac{\partial {\mathit{\zeta }}_{\mathrm{2}}}{\partial t}}-{\mathit{U}}_{\mathrm{2}}\cdot \mathrm{\nabla }{\mathit{\zeta }}_{\mathrm{2}}+f_{\mathrm{c}}\mathbf{J}\cdot {\mathit{\zeta }}_{\mathrm{2}}-{\mathit{\nu }}_{\mathrm{t},\mathrm{2}}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\zeta }}_{\mathrm{2}}+{\displaystyle \frac{{\mathbf{D}}^{\prime }}{H_{\mathrm{2}}}}\\ & \cdot {\mathit{\zeta }}_{\mathrm{2}}-{\displaystyle \frac{{\mathbf{D}}^{\prime }}{H_{\mathrm{1}}}}\cdot {\mathit{\zeta }}_{\mathrm{1}}-H_{\mathrm{2}}\left[{\mathcal{F}}^{-\mathrm{1}}\left(\widehat{\mathrm{\Phi }}\right)(-\mathit{x},-t)\ast \mathrm{\nabla }{\mathrm{\Pi }}_{\mathrm{2}}\right]\\ & =\mathrm{0}\mathrm{in}\mathrm{\Omega }\times (\mathrm{0},T],\\ & -{\displaystyle \frac{\partial {\mathrm{\Pi }}_{\mathrm{1}}}{\partial t}}-{\mathit{U}}_{\mathrm{1}}\cdot \mathrm{\nabla }{\mathrm{\Pi }}_{\mathrm{1}}-\mathrm{\nabla }\cdot {\mathit{\zeta }}_{\mathrm{1}}-\mathrm{\nabla }\cdot {\mathit{\zeta }}_{\mathrm{2}}\\ & =\mathrm{0}\mathrm{in}\mathrm{\Omega }\times (\mathrm{0},T],\\ & -{\displaystyle \frac{\partial {\mathrm{\Pi }}_{\mathrm{2}}}{\partial t}}-{\mathit{U}}_{\mathrm{2}}\cdot \mathrm{\nabla }{\mathrm{\Pi }}_{\mathrm{2}}-\mathrm{\nabla }\cdot {\mathit{\zeta }}_{\mathrm{1}}-\mathrm{\nabla }\cdot {\mathit{\zeta }}_{\mathrm{2}}\\ & =\mathrm{0}\mathrm{in}\mathrm{\Omega }\times (\mathrm{0},T].\end{array}\end{array}$$

Note that the minus sign in the argument of ${\mathcal{F}}^{-\mathrm{1}}\left(\widehat{\mathrm{\Phi }}\right)(-\mathit{x}$, −t) is not a result of classical integration by parts, but rather arrives from applying Fubini's theorem to the convolution term in Eqs. (12) and (13) (see Appendix A2 for details). The adjoint variables are grouped in the vector ${\mathit{\psi }}^{*}=[{\mathit{\zeta }}_{\mathrm{1}}$, ζ₂, Π₁, Π₂], where ζ₁=(ζ₁, χ₁) and ζ₂=(ζ₂, χ₂) are the adjoint velocity perturbation fields in the wind farm and upper layer, respectively, while Π₁ and Π₂ 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)

$$\begin{array}{}\text{(19)}& \begin{array}{rl}\mathrm{\nabla }\stackrel{\mathrm{̃}}{\mathcal{J}}=& {\displaystyle \frac{\mathit{\beta }B(x,y)}{H_{\mathrm{1}}}}\left[\parallel {\mathit{U}}_{\mathrm{1}}\parallel {\mathit{U}}_{\mathrm{1}}\cdot {\mathit{\zeta }}_{\mathrm{1}}-H_{\mathrm{1}}\parallel {\mathit{U}}_{\mathrm{1}}\parallel {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{p}}}{\mathrm{d}{C}_{\mathrm{T}}}}\right.\\ & \left.\left(\parallel {\mathit{U}}_{\mathrm{1}}{\parallel }^{\mathrm{2}}+\mathrm{3}{\mathit{U}}_{\mathrm{1}}\cdot {\mathit{u}}_{\mathrm{1}}\right)+{\mathit{u}}_{\mathrm{1}}^{\top }\cdot {\mathbf{U}}^{\prime }\cdot {\mathit{\zeta }}_{\mathrm{1}}\right],\end{array}\end{array}$$

where $\mathrm{d}{C}_{\mathrm{p}}/\mathrm{d}{C}_{\mathrm{T}}$ is computed from Eq. (14). To compute the gradient $\mathrm{\nabla }\stackrel{\mathrm{̃}}{\mathcal{J}}$, we need to solve the forward and backward equations. Since the cost for solving the adjoint equations is roughly the same as for the forward equation, the computational cost for evaluating $\mathrm{\nabla }\stackrel{\mathrm{̃}}{\mathcal{J}}$ is proportional to the cost of solving twice the state equations. To verify the approach, we compare the adjoint gradient to a standard finite-difference approximation in Appendix A4.

3 Numerical setup and case description

We define the model setup used to assess the potential of set-point 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.

3.1 Numerical setup

Both the forward and adjoint equations are discretized using a Fourier–Galerkin spectral method in space and time. The use of Fourier modes in time automatically results in satisfying the periodicity conditions that we are aiming for in our optimization setup. All terms in the equations are linear, except for the 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 $\mathrm{3}/\mathrm{2}$ rule (Canuto et al., 1988). The use of Fourier modes in space forces periodic boundary conditions at the edges of the computational domain. Therefore, the domain has a sufficiently large dimension $D_x\times D_y=\mathrm{1000}\times \mathrm{400}$ km², so that the perturbations die out before being recycled. The grid has $N_x\times N_y=\mathrm{4000}\times \mathrm{1600}$ grid points, which corresponds to a space resolution of Δ=250 m or 6.4×10⁶ 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_xN_yN_t×6N_xN_yN_t, which are solved using the LGMRES algorithm (Baker et al., 2005).

Figure 1 Convergence of the cost functional over the (a) L-BFGS-B and (b) TNC iteration. The squares and triangles denote the cost functional value and the number of function evaluations, respectively.

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). Figure 1 shows that the cost function decreases rapidly in the first two to three algorithm iterations, reaching a plateau afterwards. The use of a 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 $\stackrel{\mathrm{̃}}{\mathcal{J}}(C_{\mathrm{T}}+\mathit{\alpha }\mathit{\delta }{C}_{\mathrm{T}})$ for all directions δC_T in the control space (at the expense of solving an auxiliary set of equations). In particular, Fig. 1a shows that the cost functional is converged after six 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⁶ 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.

3.2 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 monotonically 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 $\mathrm{3}/\mathrm{2}$, while being negligible for a very wide but short farm, and vice versa (i.e. negligible for $L_y/L_x\ll \mathrm{1}$ and $L_x/L_y\gg \mathrm{1}$). Since we are interested in optimal 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 Allaerts et al. (2018) 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=\mathrm{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 Allaerts and Meyers, 2018). Note that we do not define a specific layout or a number of turbines, but we only fix the density of turbines in the farm. The turbine dimensions are based on a DTU 10 MW IEA wind turbine (Bortolotti et al., 2019) with rotor diameter D=198 m and turbine hub height z_h=119 m.

3.3 Background state variables

The governing equations are linearized around a constant background state. To determine this state, we need vertical profiles of potential temperature, velocity, shear stress and eddy viscosity plus the surface roughness z₀ 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 θ₀=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⁻¹ 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 ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}}=z_{\mathrm{0}}/H$ and the non-dimensional boundary-layer height $h_{*}=H{f}_{\mathrm{c}}/u_{*}$ are the input parameters of Nieuwstadt model, where f_c is the Coriolis frequency and $H=H_{\mathrm{1}}+H_{\mathrm{2}}$ is the ABL height. The wind-farm layer height is assumed to be twice the turbine hub height, so H₁=2z_h. The ABL height is fixed to H=1000 m and the friction velocity is set to $u_{*}=\mathrm{0.6}$ m s⁻¹. Finally, a surface roughness of $z_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{1}}$ m is adopted. Using $h_{*}=\mathrm{0.166}$ and ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{4}}$ as input values for the Nieuwstadt model, we derive the velocity U₁, U₂, 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₁, all other physical quantities are depth-averaged over the height H₁ and H₂. 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₁=0 m s⁻¹).

The pressure gradient strengths induced by inversion and internal gravity waves are dependent upon the Froude number $Fr=U_{\mathrm{B}}/\sqrt{g^{\prime }H}$ and a non-dimensional group $P_{\mathrm{N}}=U_{\mathrm{B}}^{\mathrm{2}}/NH\parallel {\mathit{U}}_{\mathrm{g}}\parallel$, respectively (Smith, 2010; Allaerts and Meyers, 2019), where the velocity scale U_B is defined as

$$\begin{array}{}\text{(20)}& U_{\mathrm{B}}={\left({\displaystyle \frac{H_{\mathrm{1}}}{H}}{\displaystyle \frac{\mathrm{1}}{U_{\mathrm{1}}^{\mathrm{2}}}}+{\displaystyle \frac{H_{\mathrm{2}}}{H}}{\displaystyle \frac{\mathrm{1}}{U_{\mathrm{2}}^{\mathrm{2}}}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

The chosen background state defines a Froude number of 0.9 for Δθ=5.54 K, which implies subcritical flow conditions (Fr<1), and a Froude number of 1.1 for Δθ=3.7 K, which leads to supercritical flow conditions (Fr>1). Further, P_N expresses the impact of internal waves in the troposphere, which increases when P_N decreases. The background state defined in Table 1 leads to P_N=1.92. The numerical setup, wind-farm layout and background state variables are summarized in Table 1.

Table 1 Numerical setup, wind-farm layout and atmospheric state used in this manuscript.

4 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.

4.1 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 time-periodic optimal thrust-coefficient distribution over the wind-farm area in a fixed time interval [0, T]. However, all optimal 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 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 wind-farm-induced gravity waves.

Figure 2 Planform view of inversion-layer displacement (a, d), pressure perturbation (b, e) and relative velocity reduction (c, f) in the wind-farm layer in subcritical (a–c, Fr=0.9) and supercritical (d–f, Fr=1.1) flow conditions. The black rectangle indicates the wind-farm region.

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_{\mathrm{T}}^{\mathrm{O}}(x$, y) are compared with a reference case. The setup of the reference model is the one reported in Table 1, but instead a uniform thrust set-point distribution over the wind-farm area is used, with $C_{\mathrm{T}}^{\mathrm{R}}(x$, $y)=C_{\mathrm{T}}^{\mathrm{Betz}}=\mathrm{8}/\mathrm{9}$.

Figure 2 illustrates a planform view of the perturbation flow patterns obtained with Fr=0.9 (top row) and Fr=1.1 (bottom row) using the reference model setup. The farm extracts energy from the flow, causing a momentum sink in the 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 wind-farm exit region and assumes a wavy behaviour 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. 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 sub- and supercritical cases, favourable pressure gradients reduce the velocity deficits in the bulk of the farm. Finally, Fig. 2c and f illustrate relative velocity reductions in the 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 Fr=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).

Figure 3 Planform view of (a) optimal thrust-coefficient distribution in subcritical (Fr=0.9) and (b) supercritical flow conditions (Fr=1.1). The length and width of the wind farm are 20 and 30 km, respectively.

Figure 4 Streamwise profiles of optimal thrust set-point distribution ($C_{\mathrm{T}}^{\mathrm{O}}$), reference (${\mathit{\eta }}_{\mathrm{t}}^{\mathrm{R}}$) and optimal (${\mathit{\eta }}_{\mathrm{t}}^{\mathrm{O}}$) inversion-layer displacement in (a) subcritical (Fr=0.9), and (b) supercritical flow conditions (Fr=1.1). The wind-farm region is marked by vertical dashed lines, and the profiles have been obtained through the centre of the farm (y=0).

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 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 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_{\mathrm{T}}^{\mathrm{O}}$ is almost invariant along the spanwise direction.

We denote with ${\mathcal{P}}^{\mathrm{R}}=\stackrel{\mathrm{̃}}{{\mathcal{J}}^{\mathrm{R}}}/T$ and ${\mathcal{P}}^{\mathrm{O}}={\stackrel{\mathrm{̃}}{\mathcal{J}}}^{\mathrm{O}}/T$ the power extracted using $C_{\mathrm{T}}=C_{\mathrm{T}}^{\mathrm{R}}=\mathrm{8}/\mathrm{9}$ and $C_{\mathrm{T}}=C_{\mathrm{T}}^{\mathrm{O}}$, respectively. Further, we define

$$\begin{array}{}\text{(21)}& \mathcal{G}={\displaystyle \frac{{\mathcal{P}}^{\mathrm{O}}-{\mathcal{P}}^{\mathrm{R}}}{{\mathcal{P}}^{\mathrm{R}}}},\end{array}$$

where 𝒢 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.

Figure 5 Streamwise profiles of (a, c) reference (p^R) and optimal (p^O) pressure perturbation and (b, d) reference ($u_{\mathrm{1}}^{\mathrm{R}}$) and optimal ($u_{\mathrm{1}}^{\mathrm{O}}$) velocity perturbation in subcritical (a, b, Fr=0.9) and supercritical (c, d, 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).

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_{\mathrm{T}}^{\mathrm{O}}$ through the centre of the farm for Fr=0.9 and Fr=1.1. To reduce gravity-wave excitation, $C_{\mathrm{T}}^{\mathrm{O}}$ is seen to be inversely related with η_t. In fact, Fig. 4a shows that the streamwise profile of η_t has a sinusoidal behaviour. Hence, the optimal set-point distribution is sinusoidal as well, explaining 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_{\mathrm{T}}^{\mathrm{O}}$ (see Fig. 3b). Moreover, Fig. 2a and d show that the gradient of η_t along the spanwise direction is much smaller than the one along the streamwise direction, explaining the almost constant thrust set-point distributions along the y direction. Figure 4 also shows that ${\mathit{\eta }}_{\mathrm{t},\max }^{\mathrm{O}}<{\mathit{\eta }}_{\mathrm{t},\max }^{\mathrm{R}}$ 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 ${\mathit{\eta }}_{\mathrm{t}}^{\mathrm{R}}$ and ${\mathit{\eta }}_{\mathrm{t}}^{\mathrm{O}}$ in this region, a displacement reduction of 14.5 % and 16.8 % is attained with the optimal configuration for the sub- and supercritical cases, respectively.

A lower vertical displacement of the inversion layer reduces 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_{\mathrm{T}}^{\mathrm{R}}$. Figure 5a and c confirm this hypothesis, showing streamwise profiles of pressure perturbations p^R and p^O through the centre of the farm for Fr=0.9 and Fr=1.1. The pressure peak is located at the entrance of the farm and a pressure peak reduction of 14.3 % and 16.2 % is attained with the optimal configuration for the sub- and supercritical cases, respectively. Figure 5b and d show streamwise profiles of velocity perturbations $u_{\mathrm{1}}^{\mathrm{R}}$ and $u_{\mathrm{1}}^{\mathrm{O}}$ through the centre of the farm for Fr=0.9 and Fr=1.1. The lower adverse pressure gradient strength attained with the optimal configuration allows for a lower velocity perturbation u₁ in the induction region with respect to the reference case. Moreover, the optimal configuration also reduces the streamline divergence, accounting for higher flow wind speeds through the farm. Consequently, a velocity perturbation reduction of 13.4 % and 15.5 % is attained for the sub- and supercritical cases, which explains the higher power gain obtained for Fr=1.1. The relative change in percentage between optimal and reference maximum flow perturbation values is summarized in Table 2.

Table 2 Relative change in percentage between optimal and reference maximum flow perturbation values. Power gains are also included.

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.

4.2 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 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}_{\mathrm{c}}/u_{*}$. Values of $h_{*}\approx \mathrm{0.1}$ denote shallow boundary layers typically found over sea, while $h_{*}\approx \mathrm{0.35}$ rather relates to a deep land-based boundary layer. We vary h_* between 0.16 and 0.4.

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

• The non-dimensional Brunt–Väisälä frequency $N/f_{\mathrm{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⁻¹. Low and high Γ values are associated with weakly and strongly stratified atmospheres, respectively. We vary Γ between 0.03 and 12 K km⁻¹ corresponding to $\mathrm{10}\le N/f_{\mathrm{c}}\le \mathrm{200}$.

• The inversion parameter $g^{\prime }H/A{u}_{*}^{\mathrm{2}}$. According to Csanady (1974), the height of the inversion layer is determined by a balance of surface stress and buoyancy. Equilibrium conditions are reached when $g^{\prime }H/A{u}_{*}^{\mathrm{2}}\approx \mathrm{1}$, with A=500 being an empirical constant. We vary the inversion parameter between 0.5 and 1.5.

Allaerts and Meyers (2019) conducted a similar sensitivity study on the 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_{\mathrm{1}}/U_{\mathrm{r}}<\mathrm{1}$ (U_r=11 m s⁻¹ 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_{\mathrm{1}}/U_{\mathrm{r}}\le \mathrm{0.9}$ for all atmospheric states.

Using the atmospheric state reported in Table 1, the non-dimensional numbers assume values $h_{*}=\mathrm{0.166}$, ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{4}}$ and $N/f_{\mathrm{c}}=\mathrm{58}$. The inversion parameter is equal to 1.046 and 0.691 in the sub- and supercritical cases, respectively. The optimal 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_*, ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}}$ and $N/f_{\mathrm{c}}$ against the inversion parameter $g^{\prime }H/A{u}_{*}^{\mathrm{2}}$, similarly to Allaerts and Meyers (2019). The numerical setup is the one detailed in Table 1. However, we use a grid cell size which is 4 times bigger ($\mathrm{\Delta }x\times \mathrm{\Delta }y=\mathrm{1000}\times \mathrm{1000}$ m²), meaning that we use 4×10⁵ cells instead of 6.4×10⁶, 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₁=0 m s⁻¹).

Figure 6 Sensitivity of (a–c) Froude number $Fr=U_{\mathrm{B}}/\sqrt{g^{\prime }H}$ and (d–f) $P_{\mathrm{N}}=U_{\mathrm{B}}^{\mathrm{2}}/NH\parallel {\mathit{U}}_{\mathrm{g}}\parallel$ to atmospheric conditions. (a, d) Non-dimensional boundary layer height h_* (with ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{4}}$ and $N/f_{\mathrm{c}}=\mathrm{58}$), (b, d) logarithm of the non-dimensional surface roughness length ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}}$ (with $h_{*}=\mathrm{0.166}$ and $N/f_{\mathrm{c}}=\mathrm{58}$) and (c, f) ratio of Brunt–Väisälä frequency to Coriolis parameter $N/f_{\mathrm{c}}$ (with $h_{*}=\mathrm{0.166}$ and ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{4}}$) against the inversion parameter $g^{\prime }H/A{u}_{*}^{\mathrm{2}}$. The black solid lines in (a–c) correspond to critical flow conditions (Fr=1) while the dashed black ones in (d–f) correspond to flow conditions of P_N=1.5. The markers △ and ∘ represent the sub- and supercritical flow cases studied in Section 4.1, respectively. Note that (f) uses a different scale than (d, e).

To better understand the power gain sensitivity to atmospheric conditions, we examine how the non-dimensional parameters Fr and P_N impact the flow fields. The pressure gradients induced by inversion waves scale with g^′; therefore high inversion strengths correspond to strong inversion-wave feedback and low Froude number values. These two-dimensional waves are non-dispersive with phase speed $\sqrt{g^{\prime }H}$ (Sutherland, 2010). Therefore, Fr also represents the ratio of the bulk wind speed within the ABL to the velocity of the inversion waves. If Fr<1 (subcritical flow) the two-dimensional waves can affect the upstream flow, while they can travel only downstream if Fr>1 (supercritical flow). The flow is said to be critical when Fr=1. On the other hand, 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 Fr 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 Fr≃1 (choking effect; Smith, 2010). Moreover, the perturbations' magnitude is strongly dependent upon the Froude number.

Smith (2010) and Allaerts and Meyers (2019) defined a third regime where N=0 and $g^{\prime }=\mathrm{0}$, which would correspond to Fr, P_N→∞ or a purely neutral atmosphere. Gravity waves are not excited in this particular flow condition, and only drag forces and frictional effects play a role in the flow behaviour. Since we are interested in finding optimal thrust 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 Fr to changes in h_*, ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}}$ and $N/f_{\mathrm{c}}$ against the inversion parameter. In all cases, the Froude number ranges from approximatively 0.5 to 1.4. The black line denotes critical flow conditions. Lines of constant Froude number run parallel to this line, meaning that Fr is invariant and quasi-invariant to $N/f_{\mathrm{c}}$ and h_*, respectively. On the other hand, changes in ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}}$ have a strong impact on the wind profile convexity and therefore on Fr. The sensitivity of P_N to the atmospheric state is displayed in Fig. 6d–f. P_N is not dependent on the inversion parameter. Hence, lines of constant P_N values are vertical and parallel to the dashed black line, which denotes atmospheric conditions for which P_N=1.5. This line divides the domain in regions where the 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 Fr and P_N values found over the sensitivity domain).

Figure 7 Sensitivity of (a–c) inversion-layer vertical-displacement reduction 𝒢_η and (d–f) power gain 𝒢 to atmospheric conditions. (a, d) Non-dimensional boundary layer height h_* (with ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{4}}$ and $N/f_{\mathrm{c}}=\mathrm{58}$), (b, d) logarithm of the non-dimensional surface roughness length ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}}$ (with $h_{*}=\mathrm{0.166}$ and $N/f_{\mathrm{c}}=\mathrm{58}$) and (c, f) ratio of Brunt–Väisälä frequency to Coriolis parameter $N/f_{\mathrm{c}}$ (with $h_{*}=\mathrm{0.166}$ and ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{4}}$) against the inversion parameter $g^{\prime }H/A{u}_{*}^{\mathrm{2}}$. The black solid line corresponds to critical flow conditions (Fr=1) while the black dashed line corresponds to flow conditions of P_N=1.5. The markers △ and ∘ represent the sub- and supercritical flow case studied in Sect. 4.1, respectively. Note that (c) and (f) use a different scale than (a, b) and (d, e).

Figure 7a–c and d–f illustrate the sensitivity of the optimal inversion-layer vertical-displacement reduction 𝒢_η and power gain 𝒢 to changes in h_*, ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}}$ and $N/f_{\mathrm{c}}$, against the inversion parameter. The displacement reduction is defined as

$$\begin{array}{}\text{(22)}& {\mathcal{G}}_{\mathit{\eta }}=\left|{\displaystyle \frac{{\mathit{\eta }}_{\mathrm{t},\max }^{\mathrm{O}}-{\mathit{\eta }}_{\mathrm{t},\max }^{\mathrm{R}}}{{\mathit{\eta }}_{\mathrm{t},\max }^{\mathrm{R}}}}\right|,\end{array}$$

where ${\mathit{\eta }}_{\mathrm{t},\max }^{\mathrm{O}}$ and ${\mathit{\eta }}_{\mathrm{t},\max }^{\mathrm{R}}$ 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, 2018, 2019) found that for low ABL heights, gravity waves induce strong pressure gradients and play an important role in the distribution of the kinetic energy within the farm. Indeed, the large geostrophic wind angle found in shallow boundary layers redirects the favourable pressure gradient seen over the wind-farm area of 90^∘ for $h_{*}\to \mathrm{0}$, decreasing the dispersive impact of internal gravity waves. Figure 7d shows that the maximum power gain is indeed attained for $h_{*}=\mathrm{0.17}$ (i.e. for shallow boundary layer) in supercritical flow conditions, with gains of about 7.5 % corresponding to a displacement reduction of 17.5 %. A similar pattern is seen in Fig. 7e, where a maximum power gain of 8.4 % is attained again in supercritical conditions for ${\log }_{\mathrm{10}}\left({\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}}\right)=-\mathrm{4.2}$, corresponding to a displacement reduction of 19%. Both 𝒢 and 𝒢_η show higher sensitivity to changes in ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}}$ than in h_*, decreasing rapidly for increasing value of surface roughness. Interestingly, power gains are close to zero in the case of high ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{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 𝒢 and 𝒢_η 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_{\mathrm{c}}$) accounts for a clear distinction between regime 1 and regime 2. The former shows power gains of about 5 % while the latter attains gains of 14 % corresponding to inversion displacement reductions of 24 %. Figure 7f illustrates that the power gain peak is obtained in critical flow conditions (Fr=1), differently from the previous cases. The very high P_N values (hence, the limited presence of internal waves) attained when Fr=1 allow for the choking effect to take place (Smith, 2010; Allaerts and Meyers, 2019). Very large flow perturbations are triggered in these atmospheric conditions, leaving greater potential for power recovery. The choking effect is not visible in Fig. 7d and e since there P_N≈2 when Fr=1 (the flow perturbations are softened by internal waves).

Overall, higher 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_*, ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}}$ and $N/f_{\mathrm{c}}$. This is not surprising due to the strong impact that gravity waves have on a farm's performance in such conditions (see Sect. 4.1 or Smith, 2010, and Allaerts and Meyers, 2019). Moreover, we observe strong gradients of 𝒢 and 𝒢_η along contours of P_N in regime 2 and weak gradients in regime 1. This suggests that the flow properties are Fr independent for low P_N values, confirming the observations of Smith (2010).

5 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 Allaerts et al. (2018) 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 (Fr=0.9) and supercritical (Fr=1.1) flow conditions. The optimal configurations were then compared with a reference model setup which uses a uniform 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 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_{\mathrm{T}}^{\mathrm{O}}$ 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 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 power recovery purposes. Power gains of up to 14 % were found for weakly stratified atmospheres (P_N≈11) in correspondence with critical flow conditions (Fr=1). This is related to the large flow perturbations induced by the choking effect (Smith, 2010). Overall, power gains above 4 % were observed for 77 % of the cases. We note that the gravity-wave-induced power losses are also sensitive to the wind-farm layout. Optimization of layout (including, e.g. relevant techno-economical 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 wind-farm 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 (Allaerts and Meyers, 2017, 2018). Finally, we assumed that the free atmosphere is uniformly stratified and steady. The relaxation of these assumptions would extend the applicability of the model, e.g. to atmospheres with 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 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² inner product over the time interval [0, T] and simulation domain Ω:

$$\begin{array}{}\text{(A1)}& (\mathit{a},\mathit{b})=\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathrm{\Omega }}{\iint }\mathit{a}\cdot \mathit{b}\mathrm{d}\mathit{x}\mathrm{d}t,\end{array}$$

where a and b are two generic vectors. Moreover, we denote with ψ=[u₁, v₁, u₂, v₂, p₁, p₂] 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

$$\begin{array}{}\text{(A2)}& \stackrel{\mathrm{̃}}{\mathcal{J}}\left(C_{\mathrm{T}}\right)=\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathrm{\Omega }}{\iint }\mathcal{K}\left(\mathit{\psi }\left(C_{\mathrm{T}}\right),C_{\mathrm{T}}\right)\mathrm{d}\mathit{x}\mathrm{d}t,\end{array}$$

where

$$\begin{array}{}\text{(A3)}& \mathcal{K}\left(\mathit{\psi }\left(C_{\mathrm{T}}\right),C_{\mathrm{T}}\right)=-\mathit{\beta }\parallel {\mathit{U}}_{\mathrm{1}}\parallel C_{\mathrm{p}}B(x,y)\left(\parallel {\mathit{U}}_{\mathrm{1}}{\parallel }^{\mathrm{2}}+\mathrm{3}{\mathit{U}}_{\mathrm{1}}\cdot {\mathit{u}}_{\mathrm{1}}\right).\end{array}$$

The gradient of the reduced cost functional $\mathrm{\nabla }\stackrel{\mathrm{̃}}{\mathcal{J}}$ is interpreted as the Riesz representation of the Gâteaux derivative operator at C_T in any arbitrary direction δC_T:

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{\mathcal{J}}}_{C_{\mathrm{T}}}\left(\mathit{\delta }{C}_{\mathrm{T}}\right)& \equiv {\displaystyle \frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}}\stackrel{\mathrm{̃}}{\mathcal{J}}\left(C_{\mathrm{T}}+\mathit{\alpha }\mathit{\delta }{C}_{\mathrm{T}}\right){|}_{\mathit{\alpha }=\mathrm{0}}\\ & =\left(\mathrm{\nabla }\stackrel{\mathrm{̃}}{\mathcal{J}},\mathit{\delta }{C}_{\mathrm{T}}\right)\forall \mathit{\delta }{C}_{\mathrm{T}}\in \mathcal{H},\end{array}\end{array}$$

where ℋ 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 𝒩(ψ, C_T). The reduced formulation of the optimization problem implies by definition that 𝒩(ψ(C_T), C_T)=0; therefore we can write the reduced cost functional as

$$\begin{array}{}\text{(A5)}& \stackrel{\mathrm{̃}}{\mathcal{J}}\left(C_{\mathrm{T}}\right)=\mathcal{J}\left(\mathit{\psi }\left(C_{\mathrm{T}}\right),C_{\mathrm{T}}\right)+\left({\mathit{\psi }}^{*},\mathcal{N}\left(\mathit{\psi }\left(C_{\mathrm{T}}\right),C_{\mathrm{T}}\right)\right),\end{array}$$

where ${\mathit{\psi }}^{*}=[{\mathit{\zeta }}_{\mathrm{1}}$, χ₁, ζ₂, χ₂, Π₁, Π₂] denotes the vector containing the adjoint variables which play the role of Lagrange multipliers. In fact, it is easy to notice that $\stackrel{\mathrm{̃}}{\mathcal{J}}\left(C_{\mathrm{T}}\right)=\mathcal{L}(C_{\mathrm{T}}$, ψ(C_T), ψ^*), where ℒ is the Lagrangian of the optimization problem in Eq. (16).

Using Eqs. (A4) and (A5), the gradient of the reduced cost functional can be expressed as

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}\left(\mathrm{\nabla }\stackrel{\mathrm{̃}}{\mathcal{J}},\mathit{\delta }{C}_{\mathrm{T}}\right)& =\left({\displaystyle \frac{\partial \mathcal{K}}{\partial C_{\mathrm{T}}}},\mathit{\delta }{C}_{\mathrm{T}}\right)+\left({\mathit{\psi }}^{*},{\displaystyle \frac{\partial \mathcal{N}}{\partial C_{\mathrm{T}}}}\mathit{\delta }{C}_{\mathrm{T}}\right)\\ & +\left({\displaystyle \frac{\partial \mathcal{K}}{\partial \mathit{\psi }}},\mathit{\delta }\mathit{\psi }\right)+\left({\mathit{\psi }}^{*},{\displaystyle \frac{\partial \mathcal{N}}{\partial \mathit{\psi }}}\mathit{\delta }\mathit{\psi }\right),\end{array}\end{array}$$

where $\mathit{\delta }\mathit{\psi }=\mathrm{d}\mathit{\psi }/\mathrm{d}{C}_{\mathrm{T}}\mathit{\delta }{C}_{\mathrm{T}}$. The adjoint of the operator $\partial \mathcal{N}/\partial \mathit{\psi }$ is given by

$$\begin{array}{}\text{(A7)}& \left({\mathit{\psi }}^{*},{\displaystyle \frac{\partial \mathcal{N}}{\partial \mathit{\psi }}}\mathit{\delta }\mathit{\psi }\right)=\left({\left[{\displaystyle \frac{\partial \mathcal{N}}{\partial \mathit{\psi }}}\right]}^{*}{\mathit{\psi }}^{*},\mathit{\delta }\mathit{\psi }\right)+{\mathrm{BT}}_{\mathrm{1}},\end{array}$$

where the right-hand side is found using integration by parts. Similarly, the adjoint of $\partial \mathcal{N}/\partial C_{\mathrm{T}}$ is expressed as

$$\begin{array}{}\text{(A8)}& \left({\mathit{\psi }}^{*},{\displaystyle \frac{\partial \mathcal{N}}{\partial C_{\mathrm{T}}}}\mathit{\delta }{C}_{\mathrm{T}}\right)=\left({\left[{\displaystyle \frac{\partial \mathcal{N}}{\partial C_{\mathrm{T}}}}\right]}^{*}{\mathit{\psi }}^{*},\mathit{\delta }{C}_{\mathrm{T}}\right)+{\mathrm{BT}}_{\mathrm{2}}.\end{array}$$

The boundary terms BT₁ and BT₂ arise as a result of the integration by parts. Due to spatial and time periodicity constraints, it is easy to show that ${\mathrm{BT}}_{\mathrm{1}}={\mathrm{BT}}_{\mathrm{2}}=\mathrm{0}$. Hence, substituting Eqs. (A7) and (A8) into Eq. (A6), we obtain

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}\left(\mathrm{\nabla }\stackrel{\mathrm{̃}}{\mathcal{J}},\mathit{\delta }{C}_{\mathrm{T}}\right)& =\left({\displaystyle \frac{\partial \mathcal{K}}{\partial C_{\mathrm{T}}}}+{\left[{\displaystyle \frac{\partial \mathcal{N}}{\partial C_{\mathrm{T}}}}\right]}^{*}{\mathit{\psi }}^{*},\mathit{\delta }{C}_{\mathrm{T}}\right)\\ & +\left({\displaystyle \frac{\partial \mathcal{K}}{\partial \mathit{\psi }}}+{\left[{\displaystyle \frac{\partial \mathcal{N}}{\partial \mathit{\psi }}}\right]}^{*}{\mathit{\psi }}^{*},\mathit{\delta }\mathit{\psi }\right).\end{array}\end{array}$$

Further, we assume that the adjoint variables satisfy the following relation:

$$\begin{array}{}\text{(A10)}& \left({\displaystyle \frac{\partial \mathcal{K}}{\partial \mathit{\psi }}}+{\left[{\displaystyle \frac{\partial \mathcal{N}}{\partial \mathit{\psi }}}\right]}^{*}{\mathit{\psi }}^{*},\mathit{\delta }\mathit{\psi }\right)=\mathrm{0},\end{array}$$

which defines the adjoint equations. Therefore, the adjoint gradient is given by

$$\begin{array}{}\text{(A11)}& \mathrm{\nabla }\stackrel{\mathrm{̃}}{\mathcal{J}}={\displaystyle \frac{\partial \mathcal{K}}{\partial C_{\mathrm{T}}}}+{\left[{\displaystyle \frac{\partial \mathcal{N}}{\partial C_{\mathrm{T}}}}\right]}^{*}{\mathit{\psi }}^{*}.\end{array}$$

A2 Derivation of the adjoint equations

We apply relation Eq. (A7) for deriving the adjoint of the operator $\partial \mathcal{N}/\partial \mathit{\psi }$. Starting with the velocity perturbations in the wind-farm layer, we have

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}\left({\mathit{\psi }}^{*},{\displaystyle \frac{\partial \mathcal{N}}{\partial {\mathit{u}}_{\mathrm{1}}}}\mathit{\delta }{\mathit{u}}_{\mathrm{1}}\right)& =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathrm{\Omega }}{\iint }\left[{\displaystyle \frac{\partial \mathit{\delta }{\mathit{u}}_{\mathrm{1}}}{\partial t}}+{\mathit{U}}_{\mathrm{1}}\cdot \mathrm{\nabla }\mathit{\delta }{\mathit{u}}_{\mathrm{1}}\right.\\ & +f_{\mathrm{c}}\mathbf{J}\cdot \mathit{\delta }{\mathit{u}}_{\mathrm{1}}-{\mathit{\nu }}_{t,\mathrm{1}}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{\delta }{\mathit{u}}_{\mathrm{1}}+{\displaystyle \frac{{\mathbf{D}}^{\prime }}{H_{\mathrm{1}}}}\cdot \mathit{\delta }{\mathit{u}}_{\mathrm{1}}\\ & \left.+{\displaystyle \frac{{\mathbf{C}}^{\prime }}{H_{\mathrm{1}}}}\cdot \mathit{\delta }{\mathit{u}}_{\mathrm{1}}+-{\displaystyle \frac{\mathrm{1}}{H_{\mathrm{1}}}}{\displaystyle \frac{\partial {\mathit{f}}^{\left(\mathrm{1}\right)}}{\partial {\mathit{u}}_{\mathrm{1}}}}{|}_{\mathit{\delta }{\mathit{u}}_{\mathrm{1}}}\right]\\ & \cdot {\mathit{\zeta }}_{\mathrm{1}}\mathrm{d}\mathit{x}\mathrm{d}t+\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathrm{\Omega }}{\iint }\left[-{\displaystyle \frac{{\mathbf{D}}^{\prime }}{H_{\mathrm{2}}}}\cdot \mathit{\delta }{\mathit{u}}_{\mathrm{1}}\right]\\ & \cdot {\mathit{\zeta }}_{\mathrm{2}}\mathrm{d}\mathit{x}\mathrm{d}t+\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathrm{\Omega }}{\iint }\left[H_{\mathrm{1}}{\mathcal{F}}^{-\mathrm{1}}\left(\widehat{\mathrm{\Phi }}\right)\right.\\ & \left.\ast \mathrm{\nabla }\cdot \mathit{\delta }{\mathit{u}}_{\mathrm{1}}\right]{\mathrm{\Pi }}_{\mathrm{1}}\mathrm{d}\mathit{x}\mathrm{d}t,\end{array}\end{array}$$

and by computing an integration by parts we obtain

$$\begin{array}{}\text{(A13)}& \begin{array}{rl}\left({\left[{\displaystyle \frac{\partial \mathcal{N}}{\partial {\mathit{u}}_{\mathrm{1}}}}\right]}^{*}{\mathit{\psi }}^{*},\mathit{\delta }{\mathit{u}}_{\mathrm{1}}\right)& =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathrm{\Omega }}{\iint }\left[-{\displaystyle \frac{\partial {\mathit{\zeta }}_{\mathbf{1}}}{\partial t}}-{\mathit{U}}_{\mathrm{1}}\right.\\ & \cdot \mathrm{\nabla }{\mathit{\zeta }}_{\mathrm{1}}+f_{\mathrm{c}}\mathbf{J}\cdot {\mathit{\zeta }}_{\mathrm{1}}-{\mathit{\nu }}_{t,\mathrm{1}}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\zeta }}_{\mathrm{1}}\\ & +{\displaystyle \frac{{\mathbf{D}}^{\prime }}{H_{\mathrm{1}}}}\cdot {\mathit{\zeta }}_{\mathrm{1}}+{\displaystyle \frac{{\mathbf{C}}^{\prime }}{H_{\mathrm{1}}}}\cdot {\mathit{\zeta }}_{\mathrm{1}}\\ & +{\displaystyle \frac{\mathit{\beta }{C}_{\mathrm{T}}B(x,y)}{H_{\mathrm{1}}}}{\mathbf{U}}^{\prime }\cdot {\mathit{\zeta }}_{\mathrm{1}}-{\displaystyle \frac{{\mathbf{D}}^{\prime }}{H_{\mathrm{2}}}}\\ & \cdot {\mathit{\zeta }}_{\mathrm{2}}-H_{\mathrm{1}}\left[{\mathcal{F}}^{-\mathrm{1}}\left(\widehat{\mathrm{\Phi }}\right)(-\mathit{x},-t)\right.\\ & \left.\left.\ast \mathrm{\nabla }{\mathrm{\Pi }}_{\mathrm{1}}\right]\right]\cdot \mathit{\delta }{\mathit{u}}_{\mathrm{1}}\mathrm{d}\mathit{x}\mathrm{d}t.\end{array}\end{array}$$

Note that the minus sign in the argument of ${\mathcal{F}}^{-\mathrm{1}}\left(\widehat{\mathrm{\Phi }}\right)(-\mathit{x}$, −t) does not come from classical integration by parts. In fact, given three functions f, g, h∈L¹(Ω), it can be shown that

$$\begin{array}{}\text{(A14)}& \begin{array}{rl}\underset{\mathrm{\Omega }}{\int }\left[f\right(x)& \ast g\left(x\right)\left]h\right(x)\mathrm{d}x\\ & =\underset{\mathrm{\Omega }}{\int }\underset{{\mathrm{\Omega }}^{\prime }}{\int }\left[f\right(x-x^{\prime }\left)g\right(x^{\prime }\left)\mathrm{d}{x}^{\prime }\right]h\left(x\right)\mathrm{d}x\\ & =\underset{{\mathrm{\Omega }}^{\prime }}{\int }\underset{\mathrm{\Omega }}{\int }f(-(x^{\prime }-x\left)\right)h\left(x\right)\mathrm{d}xg\left(x^{\prime }\right)\mathrm{d}{x}^{\prime }\\ & =\underset{\mathrm{\Omega }}{\int }\left[f\right(-x)\ast h(x\left)\right]g\left(x\right)\mathrm{d}x,\end{array}\end{array}$$

where in the second passage we have changed the order of integration (Fubini's theorem). This property allows us to write

$$\begin{array}{}\text{(A15)}& \begin{array}{rl}-H_{\mathrm{1}}& \underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathrm{\Omega }}{\iint }\left[{\mathcal{F}}^{-\mathrm{1}}\left(\widehat{\mathrm{\Phi }}\right)\ast \mathit{\delta }{\mathit{u}}_{\mathrm{1}}\right]\cdot \mathrm{\nabla }{\mathrm{\Pi }}_{\mathrm{1}}\mathrm{d}\mathit{x}\mathrm{d}t\\ & =-H_{\mathrm{1}}\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathrm{\Omega }}{\iint }\left[{\mathcal{F}}^{-\mathrm{1}}\left(\widehat{\mathrm{\Phi }}\right)(-\mathit{x},-t)\ast \mathrm{\nabla }{\mathrm{\Pi }}_{\mathrm{1}}\right]\\ & \cdot \mathit{\delta }{\mathit{u}}_{\mathrm{1}}\mathrm{d}\mathit{x}\mathrm{d}t.\end{array}\end{array}$$

Similarly, for the velocity perturbations in the upper layer, we have that

$$\begin{array}{}\text{(A16)}& \begin{array}{rl}\left({\left[{\displaystyle \frac{\partial \mathcal{N}}{\partial {\mathit{u}}_{\mathrm{2}}}}\right]}^{*}{\mathit{\psi }}^{*},\mathit{\delta }{\mathit{u}}_{\mathrm{2}}\right)& =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathrm{\Omega }}{\iint }\left[-{\displaystyle \frac{{\mathbf{D}}^{\prime }}{H_{\mathrm{1}}}}\cdot {\mathit{\zeta }}_{\mathbf{1}}-{\displaystyle \frac{\partial {\mathit{\zeta }}_{\mathrm{2}}}{\partial t}}\right.\\ & -{\mathit{U}}_{\mathrm{2}}\cdot \mathrm{\nabla }{\mathit{\zeta }}_{\mathrm{2}}+f_{\mathrm{c}}\mathbf{J}\cdot {\mathit{\zeta }}_{\mathrm{2}}\\ & -{\mathit{\nu }}_{\mathrm{t},\mathrm{2}}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\zeta }}_{\mathrm{2}}+{\displaystyle \frac{{\mathbf{D}}^{\prime }}{H_{\mathrm{2}}}}\cdot {\mathit{\zeta }}_{\mathrm{2}}\\ & +-H_{\mathrm{2}}\left[{\mathcal{F}}^{-\mathrm{1}}\left(\widehat{\mathrm{\Phi }}\right)(-\mathit{x},-t)\right.\\ & \left.\left.\ast \mathrm{\nabla }{\mathrm{\Pi }}_{\mathrm{2}}\right]\right]\cdot \mathit{\delta }{\mathit{u}}_{\mathrm{2}}\mathrm{d}\mathit{x}\mathrm{d}t.\end{array}\end{array}$$

Following the same procedure for the pressure perturbations p₁ and p₂, we obtain

$$\begin{array}{}\text{(A17)}& \begin{array}{rl}\left({\left[{\displaystyle \frac{\partial \mathcal{N}}{\partial p_{\mathrm{1}}}}\right]}^{*}{\mathit{\psi }}^{*},\mathit{\delta }{p}_{\mathrm{1}}\right)& =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathrm{\Omega }}{\iint }\left[-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }\cdot {\mathit{\zeta }}_{\mathrm{1}}-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }\right.\\ & \left.\cdot {\mathit{\zeta }}_{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathrm{\Pi }}_{\mathrm{1}}}{\partial t}}-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\mathit{U}}_{\mathrm{1}}\cdot \mathrm{\nabla }{\mathrm{\Pi }}_{\mathrm{1}}\right]\\ & \mathit{\delta }{p}_{\mathrm{1}}\mathrm{d}\mathit{x}\mathrm{d}t\end{array}\end{array}$$

and

$$\begin{array}{}\text{(A18)}& \begin{array}{rl}\left({\left[{\displaystyle \frac{\partial \mathcal{N}}{\partial p_{\mathrm{2}}}}\right]}^{*}{\mathit{\psi }}^{*},\mathit{\delta }{p}_{\mathrm{2}}\right)& =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathrm{\Omega }}{\iint }\left[-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }\cdot {\mathit{\zeta }}_{\mathrm{1}}-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }\right.\\ & \left.\cdot {\mathit{\zeta }}_{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathrm{\Pi }}_{\mathrm{2}}}{\partial t}}-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\mathit{U}}_{\mathrm{2}}\cdot \mathrm{\nabla }{\mathrm{\Pi }}_{\mathrm{2}}\right]\\ & \mathit{\delta }{p}_{\mathrm{2}}\mathrm{d}\mathit{x}\mathrm{d}t.\end{array}\end{array}$$

Using Eq. (A10), the resulting adjoint equations correspond to

$$\begin{array}{}\text{(A19)}& \begin{array}{rl}-{\displaystyle \frac{\partial {\mathit{\zeta }}_{\mathrm{1}}}{\partial t}}& -{\mathit{U}}_{\mathrm{1}}\cdot \mathrm{\nabla }{\mathit{\zeta }}_{\mathrm{1}}+f_{\mathrm{c}}\mathbf{J}\cdot {\mathit{\zeta }}_{\mathrm{1}}-{\mathit{\nu }}_{t,\mathrm{1}}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\zeta }}_{\mathrm{1}}+{\displaystyle \frac{{\mathbf{D}}^{\prime }}{H_{\mathrm{1}}}}\\ & \cdot {\mathit{\zeta }}_{\mathrm{1}}+{\displaystyle \frac{{\mathbf{C}}^{\prime }}{H_{\mathrm{1}}}}\cdot {\mathit{\zeta }}_{\mathrm{1}}-{\displaystyle \frac{{\mathbf{D}}^{\prime }}{H_{\mathrm{2}}}}\cdot {\mathit{\zeta }}_{\mathrm{2}}-H_{\mathrm{1}}\\ & \left[{\mathcal{F}}^{-\mathrm{1}}\left(\widehat{\mathrm{\Phi }}\right)(-\mathit{x},-t)\ast \mathrm{\nabla }{\mathrm{\Pi }}_{\mathrm{1}}\right]+{\displaystyle \frac{\mathit{\beta }{C}_{\mathrm{T}}B(x,y)}{H_{\mathrm{1}}}}{\mathbf{U}}^{\prime }\\ & \cdot {\mathit{\zeta }}_{\mathrm{1}}=-{\displaystyle \frac{\partial \mathcal{K}}{\partial {\mathit{u}}_{\mathrm{1}}}}\mathrm{in}\mathrm{\Omega }\times (\mathrm{0},T],\\ & -{\displaystyle \frac{\partial {\mathit{\zeta }}_{\mathrm{2}}}{\partial t}}-{\mathit{U}}_{\mathrm{2}}\cdot \mathrm{\nabla }{\mathit{\zeta }}_{\mathrm{2}}+f_{\mathrm{c}}\mathbf{J}\cdot {\mathit{\zeta }}_{\mathrm{2}}-{\mathit{\nu }}_{t,\mathrm{2}}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\zeta }}_{\mathrm{2}}+{\displaystyle \frac{{\mathbf{D}}^{\prime }}{H_{\mathrm{2}}}}\\ & \cdot {\mathit{\zeta }}_{\mathrm{2}}-{\displaystyle \frac{{\mathbf{D}}^{\prime }}{H_{\mathrm{1}}}}\cdot {\mathit{\zeta }}_{\mathrm{1}}-H_{\mathrm{2}}\left[{\mathcal{F}}^{-\mathrm{1}}\left(\widehat{\mathrm{\Phi }}\right)(-\mathit{x},-t)\ast \mathrm{\nabla }{\mathrm{\Pi }}_{\mathrm{2}}\right]\\ & =-{\displaystyle \frac{\partial \mathcal{K}}{\partial {\mathit{u}}_{\mathrm{2}}}}\mathrm{in}\mathrm{\Omega }\times (\mathrm{0},T],\\ & -{\displaystyle \frac{\partial {\mathrm{\Pi }}_{\mathrm{1}}}{\partial t}}-{\mathit{U}}_{\mathrm{1}}\cdot \mathrm{\nabla }{\mathrm{\Pi }}_{\mathrm{1}}-\mathrm{\nabla }\cdot {\mathit{\zeta }}_{\mathrm{1}}-\mathrm{\nabla }\cdot {\mathit{\zeta }}_{\mathrm{2}}\\ & =-{\mathit{\rho }}_{\mathrm{0}}{\displaystyle \frac{\partial \mathcal{K}}{\partial p_{\mathrm{1}}}}\mathrm{in}\mathrm{\Omega }\times (\mathrm{0},T],\\ & -{\displaystyle \frac{\partial {\mathrm{\Pi }}_{\mathrm{2}}}{\partial t}}-{\mathit{U}}_{\mathrm{2}}\cdot \mathrm{\nabla }{\mathrm{\Pi }}_{\mathrm{2}}-\mathrm{\nabla }\cdot {\mathit{\zeta }}_{\mathrm{1}}-\mathrm{\nabla }\cdot {\mathit{\zeta }}_{\mathrm{2}}\\ & =-{\mathit{\rho }}_{\mathrm{0}}{\displaystyle \frac{\partial \mathcal{K}}{\partial p_{\mathrm{2}}}}\mathrm{in}\mathrm{\Omega }\times (\mathrm{0},T].\end{array}\end{array}$$

The adjoint momentum equations of the upper layer are homogeneous since the adjoint wind-farm drag force is felt only indirectly in this layer ($\partial \mathcal{K}/\partial {\mathit{u}}_{\mathrm{2}}=\mathrm{0}$). Moreover, $\partial \mathcal{K}/\partial p_{\mathrm{1}}=\partial \mathcal{K}/\partial p_{\mathrm{2}}=\mathrm{0}$. On the other hand, the adjoint momentum equations of the wind-farm layer are driven by the cost function. Using Eq. (A3), we obtain

$$\begin{array}{}\text{(A20)}& {\displaystyle \frac{\partial \mathcal{K}}{\partial {\mathit{u}}_{\mathrm{1}}}}=-\mathrm{3}\mathit{\beta }{C}_{\mathrm{p}}B(x,y)\parallel {\mathit{U}}_{\mathrm{1}}\parallel {\mathit{U}}_{\mathrm{1}}.\end{array}$$

Figure A1 illustrates a planform view of the forward and adjoint solutions in subcritical flow conditions (Fr=0.9). Both solutions are derived assuming a steady-state formulation of the optimization problem. The numerical setup, wind-farm 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 $\partial \mathcal{N}/\partial C_{\mathrm{T}}$, we need to evaluate the following inner product:

$$\begin{array}{}\text{(A21)}& \begin{array}{rl}\left({\mathit{\psi }}^{*},{\displaystyle \frac{\partial \mathcal{N}}{\partial C_{\mathrm{T}}}}\mathit{\delta }{C}_{\mathrm{T}}\right)& =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathrm{\Omega }}{\iint }\left[-{\displaystyle \frac{\mathrm{1}}{H_{\mathrm{1}}}}{\displaystyle \frac{\partial {\mathit{f}}^{\left(\mathrm{0}\right)}}{\partial C_{\mathrm{T}}}}{|}_{\mathit{\delta }{C}_{\mathrm{T}}}\right.\\ & \left.-{\displaystyle \frac{\mathrm{1}}{H_{\mathrm{1}}}}{\displaystyle \frac{\partial {\mathit{f}}^{\left(\mathrm{1}\right)}}{\partial C_{\mathrm{T}}}}{|}_{\mathit{\delta }{C}_{\mathrm{T}}}\right]\cdot {\mathit{\zeta }}_{\mathbf{1}}\mathrm{d}\mathit{x}\mathrm{d}t\\ & =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathrm{\Omega }}{\iint }\left[{\displaystyle \frac{\mathit{\beta }B(x,y)}{H_{\mathrm{1}}}}\mathit{\delta }{C}_{\mathrm{T}}\parallel {\mathit{U}}_{\mathrm{1}}\parallel {\mathit{U}}_{\mathrm{1}}\right.\\ & \left.+{\displaystyle \frac{\mathit{\beta }B(x,y)}{H_{\mathrm{1}}}}\mathit{\delta }{C}_{\mathrm{T}}{\mathbf{U}}^{\prime }\cdot {\mathit{u}}_{\mathrm{1}}\right]\cdot {\mathit{\zeta }}_{\mathrm{1}}\mathrm{d}\mathit{x}\mathrm{d}t,\end{array}\end{array}$$

which is easily rewritten as

$$\begin{array}{}\text{(A22)}& \begin{array}{rl}\left({\left[{\displaystyle \frac{\partial \mathcal{N}}{\partial C_{\mathrm{T}}}}\right]}^{*}{\mathit{\psi }}^{*},\mathit{\delta }{C}_{\mathrm{T}}\right)& =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathrm{\Omega }}{\iint }\left[{\displaystyle \frac{\mathit{\beta }B(x,y)}{H_{\mathrm{1}}}}\left(\parallel {\mathit{U}}_{\mathrm{1}}\parallel {\mathit{U}}_{\mathrm{1}}\right.\right.\\ & \left.\left.\cdot {\mathit{\zeta }}_{\mathrm{1}}+{\mathit{u}}_{\mathrm{1}}^{\top }\cdot {\mathbf{U}}^{\prime }\cdot {\mathit{\zeta }}_{\mathbf{1}}\right)\right]\mathit{\delta }{C}_{\mathrm{T}}\mathrm{d}\mathit{x}\mathrm{d}t.\end{array}\end{array}$$

Moreover, we derive the first term on the right-hand side of Eq. (A11) using Eq. (A3), which results in

$$\begin{array}{}\text{(A23)}& {\displaystyle \frac{\partial \mathcal{K}}{\partial C_{\mathrm{T}}}}=-\mathit{\beta }B(x,y)\parallel {\mathit{U}}_{\mathrm{1}}\parallel {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{p}}}{\mathrm{d}{C}_{\mathrm{T}}}}\left(\parallel {\mathit{U}}_{\mathrm{1}}{\parallel }^{\mathrm{2}}+\mathrm{3}{\mathit{U}}_{\mathrm{1}}\cdot {\mathit{u}}_{\mathrm{1}}\right).\end{array}$$

Finally, we obtain the gradient expression by substituting Eqs. (A22) and (A23) in Eq. (A11), which gives

$$\begin{array}{}\text{(A24)}& \begin{array}{rl}\mathrm{\nabla }\stackrel{\mathrm{̃}}{\mathcal{J}}=& {\displaystyle \frac{\mathit{\beta }B(x,y)}{H_{\mathrm{1}}}}\left[\parallel {\mathit{U}}_{\mathrm{1}}\parallel {\mathit{U}}_{\mathrm{1}}\cdot {\mathit{\zeta }}_{\mathrm{1}}-H_{\mathrm{1}}\parallel {\mathit{U}}_{\mathrm{1}}\parallel {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{p}}}{\mathrm{d}{C}_{\mathrm{T}}}}\right.\\ & \left.\left(\parallel {\mathit{U}}_{\mathrm{1}}{\parallel }^{\mathrm{2}}+\mathrm{3}{\mathit{U}}_{\mathrm{1}}\cdot {\mathit{u}}_{\mathrm{1}}\right)+{\mathit{u}}_{\mathrm{1}}^{\top }\cdot {\mathbf{U}}^{\prime }\cdot {\mathit{\zeta }}_{\mathrm{1}}\right].\end{array}\end{array}$$

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 Fr=0.9.

We define with

$$\begin{array}{}\text{(A25)}& \mathrm{\nabla }{\stackrel{\mathrm{̃}}{\mathcal{J}}}_{\mathrm{ADJ}}=\left(\mathrm{\nabla }\stackrel{\mathrm{̃}}{\mathcal{J}},\mathit{\delta }{C}_{\mathrm{T}}\right)\end{array}$$

the directional derivative of $\mathrm{\nabla }\stackrel{\mathrm{̃}}{\mathcal{J}}$ along δC_T, where $\mathrm{\nabla }\stackrel{\mathrm{̃}}{\mathcal{J}}$ is the gradient computed with Eq. (A24) and δC_T is a perturbation of the baseline control C_T. Using finite difference, the same directional derivative can be approximated as

$$\begin{array}{}\text{(A26)}& \mathrm{\nabla }{\stackrel{\mathrm{̃}}{\mathcal{J}}}_{\mathrm{FD}}={\displaystyle \frac{\stackrel{\mathrm{̃}}{\mathcal{J}}\left(C_{\mathrm{T}}+\mathit{\alpha }\mathit{\delta }{C}_{\mathrm{T}}\right)-\stackrel{\mathrm{̃}}{\mathcal{J}}\left(C_{\mathrm{T}}\right)}{\mathit{\alpha }}}+\mathcal{O}\left(\mathit{\alpha }\right).\end{array}$$

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 finite-precision 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

$$\begin{array}{}\text{(A27)}& {\displaystyle }& {\displaystyle }R={\displaystyle \frac{\mathrm{\nabla }{\stackrel{\mathrm{̃}}{\mathcal{J}}}_{\mathrm{ADJ}}}{\mathrm{\nabla }{\stackrel{\mathrm{̃}}{\mathcal{J}}}_{\mathrm{FD}}}},\text{(A28)}& {\displaystyle }& {\displaystyle }\mathcal{E}=\left|{\displaystyle \frac{\mathrm{\nabla }{\stackrel{\mathrm{̃}}{\mathcal{J}}}_{\mathrm{ADJ}}-\mathrm{\nabla }{\stackrel{\mathrm{̃}}{\mathcal{J}}}_{\mathrm{FD}}}{\mathrm{\nabla }{\stackrel{\mathrm{̃}}{\mathcal{J}}}_{\mathrm{FD}}}}\right|,\end{array}$$

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:

$$\begin{array}{}\text{(A29)}& \begin{array}{rl}{C}_{\mathrm{T}}^{\mathrm{B}}(x,y)& =C_{\mathrm{T}}^{\mathrm{Betz}}\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{5}}}\cos \left(k_{x}x+\mathit{\pi }\right)\right.\\ & \left.+{\displaystyle \frac{\mathrm{1}}{\mathrm{5}}}\sin \left(k_{y}y+\mathit{\pi }/\mathrm{5}\right)\right],\end{array}\end{array}$$

where $C_{\mathrm{T}}^{\mathrm{Betz}}=\mathrm{8}/\mathrm{9}$, $k_x=\mathrm{2}\mathit{\pi }/L_x$ and $k_y=\mathrm{2}\mathit{\pi }/L_y$. Ideally, we should validate the 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³ 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

$$\begin{array}{}\text{(A30)}& \mathit{\delta }{C}_{\mathrm{T}}(x,y)=\cos \left(a{k}_{x}x+\mathit{\pi }\right)+\sin \left(b{k}_{y}y+\mathit{\pi }/\mathrm{5}\right)\end{array}$$

for different values of a and b.

Results of the comparison are shown in Fig. A2. We can appreciate that for ${\mathrm{10}}^{-\mathrm{11}}\le \mathit{\alpha }\le {\mathrm{10}}^{-\mathrm{4}}$ the ratio R is very close to unity, and the relative error ℰ is on the order of 10⁻⁴, 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 $\mathit{\alpha }>{\mathrm{10}}^{-\mathrm{4}}$ due to discretization errors. We can appreciate that Fig. A2b displays a first-order truncation error in accordance with relation Eq. (A26).

Figure A1 Planform view of (a) pressure perturbation, (b) velocity perturbation in the wind-farm layer, (c) velocity perturbation in the upper layer, (d) adjoint pressure $\mathrm{\Pi }={\mathrm{\Pi }}_{\mathrm{1}}+{\mathrm{\Pi }}_{\mathrm{2}}$, (e) adjoint velocity field in the wind-farm layer and (f) adjoint velocity field in the upper layer in subcritical (Fr=0.9) flow conditions. The black rectangle indicates the wind-farm region.

Figure A2 (a) Ratio and (b) relative error between adjoint and finite-difference-based gradients.

Appendix B: Optimal uniform thrust set-point distribution

In the current section, we use the optimization framework derived in Sect. 2.2 to find an optimal uniform and steady thrust-coefficient distribution that minimizes the gravity-wave-induced blockage effects. To avoid confusion, we will denote with $C_{\mathrm{T}}^{\mathrm{O}}$ and $C_{\mathrm{T}}^{\mathrm{O},\mathrm{u}}$ the optimal non-uniform and uniform distributions, respectively. The wind-farm layout and the atmospheric state are the ones detailed in Sect. 3.

Figure B1a and b display the optimal spatially invariant $C_{\mathrm{T}}^{\mathrm{O},\mathrm{u}}$ together with the streamwise profile of $C_{\mathrm{T}}^{\mathrm{O}}$ through the centre of the farm and its averaged value over the wind-farm area $\langle C_{\mathrm{T}}^{\mathrm{O}}\rangle$ for the sub- and supercritical cases, respectively. Moreover, $C_{\mathrm{T}}^{\mathrm{R}}$ denotes the thrust distribution used in the reference model. Interestingly, $C_{\mathrm{T}}^{\mathrm{O},\mathrm{u}}$ corresponds to the average of the non-uniform distribution in both cases. Since $C_{\mathrm{T}}^{\mathrm{O}}$ is sensitive to the atmospheric conditions, we expect $C_{\mathrm{T}}^{\mathrm{O},\mathrm{u}}$ to also depend on the atmospheric state (in fact, we observe a different value of $C_{\mathrm{T}}^{\mathrm{O},\mathrm{u}}$ in sub- and supercritical conditions).

In the current example, the power gain 𝒢 (see Eq. 21) over the reference model configuration obtained with the non-uniform distributions $C_{\mathrm{T}}^{\mathrm{O}}$ is 5.3 % and 7 % for the sub- and supercritical cases, respectively. For the optimal uniform distributions, we obtain a power gain of 5 % and 6.6 %.

Figure B1 Reference thrust set point ($C_{\mathrm{T}}^{\mathrm{R}}$), optimal non-uniform thrust set point ($C_{\mathrm{T}}^{\mathrm{O}}$) and its averaged value over the wind-farm area ($\langle C_{\mathrm{T}}^{\mathrm{O}}\rangle$) and optimal uniform thrust-coefficient distribution ($C_{\mathrm{T}}^{\mathrm{O},\mathrm{u}}$) in (a) subcritical and (b) supercritical flow conditions. The $C_{\mathrm{T}}^{\mathrm{O}}$ profiles are taken through the centre of the farm (y=0).

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/\mathrm{3}$, or equivalently from 1.6×10⁴ to 3.6×10⁶ DOF per layer. The results obtained are compared with the ones derived on a finer grid with resolution equal to $H/\mathrm{4}$.

Figure C1a and b display the cost function and power gain relative error, respectively, which are computed as

$$\begin{array}{}\text{(C1)}& {\displaystyle }& {\displaystyle }{\mathcal{E}}_{\mathcal{P}}=\left|{\displaystyle \frac{{\mathcal{P}}^{\mathrm{F}}-\mathcal{P}}{\mathcal{P}}}\right|,\text{(C2)}& {\displaystyle }& {\displaystyle }{\mathcal{E}}_{\mathcal{G}}=\left|{\displaystyle \frac{{\mathcal{G}}^{\mathrm{F}}-\mathcal{G}}{\mathcal{G}}}\right|,\end{array}$$

where ${\mathcal{P}}^{\mathrm{F}}={\stackrel{\mathrm{̃}}{\mathcal{J}}}^{\mathrm{F}}/T$ and 𝒢^F are the cost function (scaled with the time horizon T) and power gain obtained with a $H/\mathrm{4}$ grid resolution while $\mathcal{P}=\stackrel{\mathrm{̃}}{\mathcal{J}}/T$ and 𝒢 are the ones obtained with coarser grids. Note that the optimal distributions are steady state; therefore ℰ_𝒫 is not dependent on the choice of the time horizon T. The cost function is evaluated using the reference case setup. The power gain is obtained using the optimization model described in Sect. 2. The model setup is reported in Table 1.

Spectral methods are known to have exponential convergence when used for discretizing smooth functions (i.e. f∈C^∞). However, algebraic convergence is attained for functions f∈C^p with p≥0. Figure C1 illustrates that we obtain a 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. 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⁴ to 10⁶). This justifies the use of a coarser grid in the sensitivity study performed in Sect. 4.2.

Figure C1 (a) Cost function and (b) power gain relative error between a grid with resolution $H/\mathrm{4}$ and coarser grids in sub- and supercritical flow conditions.