Including realistic upper atmospheres in a wind-farm gravity-wave model

Recent research suggests that atmospheric gravity waves can affect offshore wind-farm performance. A fast wind-farm boundary layer model has been proposed to simulate the effects of these gravity waves on wind-farm operation by Allaerts and Meyers (2019). The current work extends the applicability of that model to free atmospheres in which wind and stability vary with altitude. We validate the model using reference cases from literature on mountain waves. Analysis of a reference flow shows that internal gravity-wave resonance caused by the atmospheric non-uniformity can prohibit perturbations in the atmospheric boundary layer (ABL) at the wavelengths where it occurs. To determine the overall impact of the vertical variations in the atmospheric conditions on wind-farm operation, we consider 1 year of operation of the Belgian–Dutch wind-farm cluster with the extended model. We find that this impact on individual flow cases is often of the same order of magnitude as the total flow perturbation. In 16.6 % of the analyzed flows, the relative difference in upstream velocity reduction between uniform and non-uniform free atmospheres is more than 30 %. However, this impact is small when averaged over all cases. This suggests that variations in the atmospheric conditions should be taken into account when simulating wind-farm operation in specific atmospheric conditions.

1 Introduction

In recent years, it has been well documented that wind farms form a blockage to the flow in and around them (Bleeg et al., 2018), thereby displacing the atmospheric boundary layer (ABL). Such displacements can propagate through the overlying inversion layer and free atmosphere as waves in stably stratified atmospheres, conditions which frequently occur at sea. As offshore wind farms in Europe increase in size and installed capacity (WindEurope, 2018), improving the understanding and simulation of gravity-wave–wind-farm interaction becomes crucial to optimizing turbine control and wind-farm layout (see, e.g., Lanzilao and Meyers, 2021b).

Previous work on the interaction between gravity waves and wind farms has assumed the free atmosphere to be uniformly stratified, with a constant background wind (Smith, 2010; Allaerts and Meyers, 2017, 2019). However, for waves with a horizontal scale of a few kilometers, such as those triggered by wind-farm blockage, vertical variations in the atmospheric conditions can drastically influence the pressure feedback they induce (Teixeira et al., 2013). Therefore, the goal of this work is to determine how variations in the free atmosphere's properties with altitude affect the interaction between gravity waves and the flow in the ABL and overall wind-farm performance.

Allaerts and Meyers (2019) proposed a model focused on the interaction between wind turbines and gravity waves. Based on earlier work by Smith (2010), the model explicitly simulates the flow in the ABL as two height-averaged horizontal layers and incorporates the free atmosphere as a boundary condition. This boundary condition links variations in ABL height to the pressure gradients induced by the corresponding gravity waves. In Fourier space, the relation between these height and pressure perturbations is defined by the stratification coefficient $\widehat{\mathrm{\Phi }}$. As the model is based on a three-layer representation of the atmosphere, it is called the three-layer model (TLM).

Currently, the TLM can only describe uniformly stratified free atmospheres, which places a strong restriction on the atmospheric conditions that can be represented. This work adapts the TLM for flow profiles that vary with altitude and studies how these variations change the interaction between the ABL flow and gravity waves. A common approach in gravity-wave theory is to use a piecewise representation of the upper atmosphere, where the profiles of the stratification and the wind speed are split up in a discrete number of layers (Tolstoy, 1973; Gossard and Hooke, 1975; Baines, 1998; Smith et al., 2002; Teixeira et al., 2013; Yu and Teixeira, 2015). We generalize this approach by allowing for an arbitrarily large number of layers so that any atmospheric profile can be accurately analyzed.

It is well known that variations in the atmospheric state can cause wave reflection, which might lead to internal gravity-wave resonance (Gill, 1982). Using the extended TLM, we examine how this resonance affects the pressure feedback on the ABL. Extensive research has been performed on internal gravity-wave resonance, with most of it focusing on flow around topographies, mountain waves, and the effect on the free atmosphere (Teixeira, 2014). As a result, the vertical displacement is simply given by the shape of the topography, not by a change in ABL height. In wind farms, this is different: the displacement is given by a change in ABL depth due to wind-farm blockage, in which the induced gravity waves can play an active role. We therefore set up a qualitative analysis of the feedback loop between ABL flow perturbations and the associated resonant gravity-wave feedback. Finally, in order to estimate the overall effect of vertical variations in the atmospheric conditions, 1 year of operation of the Belgian–Dutch wind-farm cluster is simulated.

The remainder of this paper is organized as follows. Section 2 reviews how gravity waves and the ABL flow interact and how the TLM models ABL flow. In Sect. 3, the analytical expressions for the stratification coefficients for non-uniform atmospheres are derived. To evaluate these, a method for solving the wave equation has to be developed, which is done in Sect. 3.2. Section 4 uses the developed method to analyze the effect of variations in the wind and stability in the free atmosphere on wind-farm operation. Finally, the paper concludes with a summary and suggestions for further research in Sect. 5.

2 Gravity-wave interaction with wind farms

Wind farms form a blockage to the ABL flowing through and around them, thereby pushing the inversion layer, and the free atmosphere above, upwards. These displacements can trigger gravity waves, which may influence the ABL flow by inducing pressure gradients. The first part of this section discusses these induced pressures, while the second part gives an overview of how the TLM models ABL flow and how it incorporates gravity-wave effects.

2.1 Pressure feedback induced by gravity waves

The vertical displacement η_t of the capping inversion topping the ABL leads to two types of atmospheric gravity waves. The first type, called inversion gravity waves, propagate horizontally along the inversion layer, similar to surface water waves. Due to the strong stratification in the inversion layer, any vertical displacement η_t of the inversion layer will be counteracted by the buoyancy-generated changes in pressure $p_{\mathrm{1}}^{\prime }$ (Gill, 1982; Smith, 2010):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{p_{\mathrm{1}}^{\prime }}{\stackrel{\mathrm{‾}}{\mathit{\rho }}}}=g^{\prime }{\mathit{\eta }}_t,\end{array}$$

where $\stackrel{\mathrm{‾}}{\mathit{\rho }}$ is the unperturbed density, and $g^{\prime }=g\mathrm{\Delta }\mathit{\theta }/\mathit{\theta }$ is the reduced gravity, determined by the potential temperature θ and the inversion strength Δθ.

The second type of waves propagates vertically through the free atmosphere above the capping inversion if it is stably stratified, as is usually the case. The free atmosphere perceives η_t similarly to large-scale topographies, and internal gravity waves are generated (Smith, 2010). In return, the internal gravity waves also induce pressure gradients on the ABL. This pressure feedback is easily expressed by using double Fourier transforms in the horizontal directions (i.e.,  $(x,y)\to (k,l)$) of the displacement and the pressure. For each wavenumber (k,l), the perturbation in the free atmosphere is a plane wave, with the vertical velocity given by Nappo (2012):

$$\begin{array}{}\text{(2)}& w^{\prime }(x,y,z)=-\mathit{ı}\mathrm{\Omega }\widehat{{\mathit{\eta }}_t}\exp \left[\mathit{ı}\right(kx+ly+mz\left)\right],\end{array}$$

where $\mathrm{\Omega }=-\mathit{\kappa }\cdot \stackrel{\mathrm{‾}}{\mathit{u}}=-(\stackrel{\mathrm{‾}}{u}k+\stackrel{\mathrm{‾}}{v}l)$ is the intrinsic frequency of the gravity waves, with $\stackrel{\mathrm{‾}}{u}$ and $\stackrel{\mathrm{‾}}{v}$ the x and y components of the background velocity, and where m is the vertical wavenumber of the internal gravity waves. For each separate wavenumber (k,l), the pressure perturbation ${\widehat{p}}_{\mathrm{2}}^{\prime }$ of the plane wave is proportional to the ABL displacement ${\widehat{\mathit{\eta }}}_t$, with the relation determined by the stratification coefficients $\widehat{\mathrm{\Phi }}$ (Smith, 2010):

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{{\widehat{p}}_{\mathrm{2}}^{\prime }}{\stackrel{\mathrm{‾}}{\mathit{\rho }}}}=\widehat{\mathrm{\Phi }}{\widehat{\mathit{\eta }}}_t.\end{array}$$

For uniform free atmospheres, these coefficients are given by Smith (2010):

$$\begin{array}{}\text{(4)}& \widehat{\mathrm{\Phi }}={\displaystyle \frac{\mathit{ı}\left(N_g^{\mathrm{2}}-{\mathrm{\Omega }}^{\mathrm{2}}\right)}{m}},\end{array}$$

where $N_g=\sqrt{\frac{g}{\mathit{\theta }}\frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}z}}$ is the Brunt–Väisälä frequency, which is constant with height in the uniform case. Further, m can be found using the dispersion relation (Gill, 1982):

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

If ${\mathrm{\Omega }}^{\mathrm{2}}<N_g^{\mathrm{2}}$, then m² is positive, and the waves propagate vertically. If ${\mathrm{\Omega }}^{\mathrm{2}}>N_g^{\mathrm{2}}$, then m² is negative, and the waves become evanescent. To find Φ, the sign of m has to be known. For propagating waves (m²>0), the sign of m is determined by the radiation condition, which states that the energy flux of a wave, which is directed along $\frac{\partial \mathrm{\Omega }}{\partial m}$, should be upwards, away from the perturbation (Baines, 1998). It then follows that $\text{sign}\left(m\right)=-\text{sign}\left(\mathrm{\Omega }\right)$. For evanescent waves (m²<0), the positive root has to be chosen for the perturbation to die out with altitude (Smith, 1980).

2.2 Three-layer model

The TLM is based on earlier work by Smith (2010), who first analyzed the impact of gravity waves on wind-farm operation. It represents the ABL as being neutrally stable and capped by an inversion layer, with the flow above being stably stratified. Although stably stratified atmospheres have often been modeled this way (Klemp and Lilly, 1975; Durran, 1990; Vosper, 2004; Smith, 2007; Sachsperger et al., 2015), Smith (2010) was the first to apply such a model to wind-farm operation. His model assumes that the flow in the ABL does not vary with height so that the perturbed variables can be replaced by their height averages. As a result, horizontal flow divergence and convergence will lead to changes in the ABL height. Another assumption is that the ABL flow is assumed to be hydrostatic. The inversion layer is modeled as a zeroth-order jump in potential temperature Δθ at the top of the ABL, and the free atmosphere above is assumed to be uniformly stratified.

The TLM improves on the model by Smith in several ways, the most important of which is to divide the ABL in two separate layers. The resulting lower and an upper layer are denoted by subscripts 1 and 2, respectively. The two layers are separated by a pliant surface, similar to the interface between the ABL and the free atmosphere. The wind-farm forcing terms are added in the momentum equations for the lower layer while only affecting the upper layer through interaction through the pliant surface. For this reason, the lower layer is also called the wind-farm layer. The resulting approximation of the ABL is visualized in Fig. 1. The flow in the two layers is governed by the following two-dimensional depth-averaged linearized momentum and continuum equations (Allaerts and Meyers, 2019):

$$\begin{array}{}\text{(6)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{u}}_{j,\mathrm{1}}{\displaystyle \frac{\partial u_{i,\mathrm{1}}^{\prime }}{\partial x_j}}& +{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathit{\rho }}}}{\displaystyle \frac{\partial \left(g^{\prime }+\mathrm{\Phi }\right)*{\mathit{\eta }}_t}{\partial x_i}}=f_{\mathrm{c}}{\mathit{ϵ}}_{ij\mathrm{3}}{u}_{j,\mathrm{1}}^{\prime }\\ & +{\mathit{\nu }}_{t,\mathrm{1}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_{i,\mathrm{1}}^{\prime }}{\partial x_j\partial x_j}}+{\displaystyle \frac{D_{ij}}{H_{\mathrm{1}}}}{\mathrm{\Delta }}_{\mathrm{1}}^{\mathrm{2}}{u}_j^{\prime }-{\displaystyle \frac{C_{ij}}{H_{\mathrm{1}}}}{u}_{j,\mathrm{1}}^{\prime }+{\displaystyle \frac{F_i}{H_{\mathrm{1}}}}\end{array}\text{(7)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{u}}_{j,\mathrm{2}}{\displaystyle \frac{\partial u_{i,\mathrm{2}}^{\prime }}{\partial x_j}}& +{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathit{\rho }}}}{\displaystyle \frac{\partial \left(g^{\prime }+\mathrm{\Phi }\right)*{\mathit{\eta }}_t}{\partial x_i}}=f_{\mathrm{c}}{\mathit{ϵ}}_{ij\mathrm{3}}{u}_{j,\mathrm{2}}^{\prime }\\ & +{\mathit{\nu }}_{t,\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_{i,\mathrm{2}}^{\prime }}{\partial x_j\partial x_j}}-{\displaystyle \frac{D_{ij}}{H_{\mathrm{1}}}}{\mathrm{\Delta }}_{\mathrm{1}}^{\mathrm{2}}{u}_j^{\prime }\end{array}\text{(8)}& {\displaystyle }{\stackrel{\mathrm{‾}}{u}}_{j,\mathrm{1}}{\displaystyle \frac{\partial {\mathit{\eta }}_{\mathrm{1}}}{\partial x_j}}+H_{\mathrm{1}}{\displaystyle \frac{\partial u_{j,\mathrm{1}}^{\prime }}{\partial x_j}}=\mathrm{0}\text{(9)}& {\displaystyle }{\stackrel{\mathrm{‾}}{u}}_{j,\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\eta }}_{\mathrm{2}}}{\partial x_j}}+H_{\mathrm{2}}{\displaystyle \frac{\partial u_{j,\mathrm{2}}^{\prime }}{\partial x_j}}=\mathrm{0}.\end{array}$$

The atmospheric base state is governed by the mean depth-averaged wind speeds ${\stackrel{\mathrm{‾}}{u}}_{i,\mathrm{1}}$ and ${\stackrel{\mathrm{‾}}{u}}_{i,\mathrm{2}}$ (with $i=\mathrm{1},\mathrm{2}$) and the layer heights H₁ and H₂ of the wind farm and upper layer, respectively, and $u_{i,\mathrm{1}}^{\prime }$, $u_{i,\mathrm{2}}^{\prime }$, η₁, and η₂ represent the perturbations to this reference state. The total inversion layer displacement η_t is given by the sum of η₁ and η₂. Following Allaerts and Meyers (2019), H₁ is taken to be twice the hub height of the turbines throughout this work. Further, ${\mathrm{\Delta }}_{\mathrm{1}}^{\mathrm{2}}{u}_j^{\prime }=u_{j,\mathrm{2}}^{\prime }-u_{j,\mathrm{1}}^{\prime }$ is the perturbation of the velocity difference between the wind farm and the upper layer, linked to the perturbation of the friction at the interface by the matrix D_ij. Similarly, the matrix C_ij relates the perturbation of the friction at the ground to the velocity perturbations in the lower layer. Finally, f_c is the Coriolis parameter, and F_i is the wind-turbine-forcing term.

Figure 1 Schematic representation of the three-layer model. Figure from Allaerts and Meyers (2019), reproduced with permission.

In Eqs. (6)–(9), the pressure has been substituted using the equations discussed in Sect. 2.1. As the ABL is assumed to be hydrostatic, the pressure perturbation is equal to the pressure induced by the gravity waves generated by the changes in ABL height ${\mathit{\eta }}_t={\mathit{\eta }}_{\mathrm{1}}+{\mathit{\eta }}_{\mathrm{2}}$:

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{p^{\prime }}{\stackrel{\mathrm{‾}}{\mathit{\rho }}}}=\left(g^{\prime }+\mathrm{\Phi }\right)*{\mathit{\eta }}_t,\end{array}$$

where * is the convolution operator. In this way, the TLM incorporates the gravity-wave effects through a pressure boundary condition, which is determined by the stratification coefficients Φ. These coefficients do not depend on the values of η_t or p^′ and can thus be calculated separately from the set of Eqs. (6) to (9). To avoid the computationally expensive convolution, the TLM is solved using a spectral method with a Fourier–Galerkin discretization (Allaerts and Meyers, 2019). Finally, we note that the hydrostatic assumption in the boundary layer ($\frac{\partial p}{\partial z}=\mathrm{0}$) is only reasonable as long as pressure effects within the ABL are negligible compared to those of the gravity waves. In particular in cases where a capping inversion is absent, we have noticed that this assumption may not be valid, which can lead to unphysical perturbations. Therefore, in the current work, we will only consider atmospheric cases with a capping inversion. Extension of the boundary layer (Eqs. 6 and 7) to include hydrodynamic effects is a topic of further research.

The turbines are represented individually using an actuator disk model. To incorporate their interactions, the TLM is coupled with a wake model, which in this work is a Gaussian wake model coupled with linear superposition of velocity deficits (Bastankhah and Porté-Agel, 2014; Niayifar and Porté-Agel, 2016). The forces f_i,k for each individual turbine k are written as a first-order Taylor expansion around the background inflow velocity ${\stackrel{\mathrm{‾}}{u}}_{\mathrm{fs}}$ in order to incorporate the effect of the velocity perturbation (Allaerts and Meyers, 2019):

$$\begin{array}{}\text{(11)}& f_{i,k}\left(u_{\mathrm{fs}}\right)=f_{i,k}\left({\stackrel{\mathrm{‾}}{u}}_{\mathrm{fs}}\right)+J_i^{f,k}\left(u_{\mathrm{fs}}^{\prime }\right),\end{array}$$

where $J_i^{f,k}$ is the Jacobian of f_i,k. The background inflow velocity is taken to be the mean wind speed in the wind-farm layer, while the velocity perturbation is evaluated at a distance of 10 D upstream of the farm, where D is the diameter of the turbines. Finally, the turbine forces are filtered on the numerical grid with a Gaussian filter (Allaerts and Meyers, 2019):

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{F}_i& =\underset{\mathrm{0}}{\overset{L_x}{\int }}\underset{\mathrm{0}}{\overset{L_y}{\int }}G(x-x^{\prime },y-y^{\prime })\\ & \sum _k^{N_t}{f}_{i,k}\mathit{\delta }(x-x_k,y-y_k)\mathrm{d}{x}^{\prime }\mathrm{d}{y}^{\prime }\end{array},\end{array}$$

where L_x×L_y is the size of the domain, (x_k,y_k) denotes the positions of the turbines, and $G(x-x^{\prime },y-y^{\prime })$ is a 2D Gaussian kernel (Allaerts and Meyers, 2019):

$$\begin{array}{}\text{(13)}& G(x,y)={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }{L}^{\mathrm{2}}}}\exp \left(-{\displaystyle \frac{x^{\mathrm{2}}+y^{\mathrm{2}}}{L^{\mathrm{2}}}}\right).\end{array}$$

We use a filter length of L=1 km.

3 Extension to vertically non-uniform free atmospheres

In reality, the free atmosphere is not uniform, and the stratification strength and wind speed can strongly depend on altitude. This will of course impact internal gravity-wave propagation through the atmosphere and thus the pressure feedback of these waves in the ABL. Currently, the TLM does not incorporate this as the simplified version of the internal wave equation on which Eq. (4) is based is only valid for uniformly stratified free atmospheres with a constant wind velocity. This section derives expressions for the stratification coefficients for vertically non-uniform atmospheres, where both the stratification strength and wind speed can vary.

3.1 Gravity waves in vertically non-uniform flows

The internal gravity-wave equation in vertically non-uniform atmospheres with continuous background velocities is (Teixeira, 2014)

$$\begin{array}{}\text{(14)}& \begin{array}{rl}\left[{\displaystyle \frac{{\mathrm{D}}^{\mathrm{2}}}{\mathrm{D}{t}^{\mathrm{2}}}}\right({\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial y^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}})& -{\displaystyle \frac{\mathrm{D}}{\mathrm{D}t}}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}\left({\displaystyle \frac{\mathrm{D}}{\mathrm{D}t}}\right)\\ & +N_g^{\mathrm{2}}({\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial y^{\mathrm{2}}}})]w^{\prime }=\mathrm{0}\end{array},\end{array}$$

where $\frac{\mathrm{D}}{\mathrm{D}t}$ is the inertial derivative, which for steady systems simplifies to $\frac{\mathrm{D}}{\mathrm{D}t}={\stackrel{\mathrm{‾}}{u}}_i\frac{\partial }{\partial x_i}$, and where w^′ is the vertical velocity perturbation of the wave. For plane waves, the solution for this equation can be written as (Teixeira, 2014):

$$\begin{array}{}\text{(15)}& w^{\prime }(x,y,z)=W\left(z\right)\exp \left[\mathit{ı}\right(kx+ly\left)\right].\end{array}$$

Equation (14) then reduces to the Helmholtz equation (Gill, 1982):

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}{z}^{\mathrm{2}}}}W\left(z\right)+m^{\mathrm{2}}W\left(z\right)=\mathrm{0},\end{array}$$

where m² is given by

$$\begin{array}{}\text{(17)}& m^{\mathrm{2}}=(k^{\mathrm{2}}+l^{\mathrm{2}})({\displaystyle \frac{N_g^{\mathrm{2}}}{{\mathrm{\Omega }}^{\mathrm{2}}}}-\mathrm{1})-{\mathrm{\Omega }}^{-\mathrm{1}}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}\mathrm{\Omega }}{\mathrm{d}{z}^{\mathrm{2}}}}.\end{array}$$

It is important to note that m² is not a constant as it depends on the altitude through N_g(z) and Ω(z). The above derivation is only valid if the vertical background velocities are continuous. In points where the background velocities or their derivatives to altitude have a discontinuity, w^′ itself can be discontinuous, as is elaborated in Sect. 3.2.2.

The relation between w^′ and p^′ is given by

$$\begin{array}{}\text{(18)}& \left[{\displaystyle \frac{\mathrm{D}}{\mathrm{D}t}}{\displaystyle \frac{\partial }{\partial z}}-\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u}}{\partial z}}{\displaystyle \frac{\partial }{\partial x}}+{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{v}}{\partial z}}{\displaystyle \frac{\partial }{\partial y}}\right)\right]w^{\prime }=({\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial y^{\mathrm{2}}}}){\displaystyle \frac{p^{\prime }}{\stackrel{\mathrm{‾}}{\mathit{\rho }}}}.\end{array}$$

For plane waves, this relation is

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{{\widehat{p}}^{\prime }}{\stackrel{\mathrm{‾}}{\mathit{\rho }}}}={\displaystyle \frac{\mathit{ı}}{k^{\mathrm{2}}+l^{\mathrm{2}}}}\left(\mathrm{\Omega }{\displaystyle \frac{\mathrm{d}W}{\mathrm{d}z}}-{\displaystyle \frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}z}}W\right).\end{array}$$

Using $W\left(H\right)=-\mathit{ı}\mathrm{\Omega }{\widehat{\mathit{\eta }}}_t$ (cf. Eq. 2), the definition of the stratification coefficient then leads to

$$\begin{array}{}\text{(20)}& \widehat{\mathrm{\Phi }}={\left.\left[{\displaystyle \frac{\mathrm{\Omega }}{k^{\mathrm{2}}+l^{\mathrm{2}}}}\left({\displaystyle \frac{\mathrm{\Omega }}{W}}{\displaystyle \frac{\mathrm{d}W}{\mathrm{d}z}}-{\displaystyle \frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}z}}\right)\right]\right|}_{z=H}.\end{array}$$

It is easily verified that the above equation and the expression for m² (Eq. 17) simplify to Eqs. (4) and (5) if the free atmosphere is uniformly stratified.

3.2 Piecewise methods

To evaluate the expression for the stratification coefficients derived in the previous section, the Helmholtz equation for vertically non-uniform atmospheres has to be solved. This is no longer trivial as the vertical wavenumber now varies with altitude. In earlier studies, one of the main approaches to solving the Helmholtz equation has been the so-called piecewise or multilayer methods (Tolstoy, 1973; Gossard and Hooke, 1975; Baines, 1998; Smith et al., 2002; Teixeira, 2014; Pütz et al., 2019). These techniques are based on modeling the continuously varying atmosphere as a discrete set of layers, in which the flow's properties allow for easy computation of the internal wave field. This approach is also commonly used in acoustics, where it is called the fast field program (FFP) method (Salomons, 2001). To avoid confusion with the lower and upper layers of the TLM, the layers used by piecewise methods are called sublayers from now on.

3.2.1 General principle

The basic principle of piecewise methods is to represent the atmosphere as a discrete number of sublayers (Tolstoy, 1973; Gossard and Hooke, 1975; Baines, 1998). In these sublayers, the atmosphere's properties should vary in such a way that solutions to the Helmholtz equation are easily found while still approximating the actual profiles. By increasing the number of sublayers, most aspects of the real flow's behavior can be accurately simulated. A suitable approach is to use sublayers in which m² has simple profiles for which analytic solutions can be found. In this work, the atmosphere is approximated in a piecewise-uniform fashion so that m² is piecewise-constant and given in each sublayer by Eq. (5). The internal wave field in each sublayer is a superposition of exponential functions, corresponding to upwards- and downwards-traveling waves.

The main advantage of this method is that realistic wave patterns can be obtained with a relatively small number of sublayers. Within a sublayer, only 2 degrees of freedom have to be determined. Another advantage compared to other methods such as Wentzel–Kramers–Brillouin (WKB) theory is that piecewise methods can account for wave reflection, although they cannot incorporate weakly non-linear effects (Gill, 1982).

As the number of sublayers has to be limited for computational reasons, not all of the atmosphere can be approximated. An appropriate height H_n has to be chosen, above which the atmosphere is considered uniform so that W(z>H_n) is an exponential function. In general, the piecewise method can be summarized by the following equations (Tolstoy, 1973; Pütz et al., 2019):

$$\begin{array}{}\text{(21)}& {\displaystyle }{m}^{\mathrm{2}}\left(z\right)\approx {\stackrel{\mathrm{̃}}{m}}^{\mathrm{2}}\left(z\right)\mathit{\triangleq }\left\{\begin{array}{cc}{m}_{\mathrm{0}}^{\mathrm{2}},& H<z<H_{\mathrm{1}}\\ \mathrm{⋮}& \mathrm{⋮}\\ m_j^{\mathrm{2}},& H_j<z<H_{j+\mathrm{1}}\\ \mathrm{⋮}& \mathrm{⋮}\\ m_n^{\mathrm{2}},& H_n<z\end{array}\right.\text{(22)}& {\displaystyle }W\left(z\right)\approx \stackrel{\mathrm{̃}}{W}\left(z\right)\mathit{\triangleq }\left\{\begin{array}{cc}{W}_{\mathrm{0}}\left(z\right)=W_{\mathrm{0}+}\exp \left(\mathit{ı}{m}_{\mathrm{0}}z\right)+W_{\mathrm{0}-}\exp (-\mathit{ı}{m}_{\mathrm{0}}z),& H<z<H_{\mathrm{1}}\\ \mathrm{⋮}& \mathrm{⋮}\\ W_j\left(z\right)=W_{j+}\exp \left(\mathit{ı}{m}_{j}z\right)+W_{j-}\exp (-\mathit{ı}{m}_{j}z),& H_j<z<H_{j+\mathrm{1}}\\ \mathrm{⋮}& \mathrm{⋮}\\ W_n\left(z\right)=W_{n+}\exp \left(\mathit{ı}{m}_{n}z\right)+W_{n-}\exp (-\mathit{ı}{m}_{n}z),& H_n<z\end{array},\right.\end{array}$$

where H_j denotes the altitudes of the interfaces between sublayers, and n is the number of interfaces between the sublayers. The values for the 2(n+1) coefficients W_j±, which determine the internal gravity-wave field, depend on additional conditions at the boundaries and sublayer interfaces, which are discussed in Sect. 3.2.2. The atmospheric variables are evaluated in the middle of each sublayer so that ${\stackrel{\mathrm{̃}}{m}}_j^{\mathrm{2}}=\left[m\right(H_j+\mathrm{\Delta }{H}_j/\mathrm{2}){]}^{\mathrm{2}}$. This results in second-order convergence when approximating continuously varying atmospheres, as shown in Appendix A. Once the coefficients W_j± have been determined, Eq. (20) for the stratification coefficients can be evaluated as

$$\begin{array}{}\text{(23)}& \widehat{\mathrm{\Phi }}={\displaystyle \frac{{\mathrm{\Omega }}_{\mathrm{0}}}{k^{\mathrm{2}}+l^{\mathrm{2}}}}\left({\mathrm{\Omega }}_{\mathrm{0}}\mathit{ı}{m}_{\mathrm{0}}{\displaystyle \frac{W_{\mathrm{0}+}-W_{\mathrm{0}-}}{W_{\mathrm{0}+}+W_{\mathrm{0}-}}}-{\displaystyle \frac{\mathrm{d}{\mathrm{\Omega }}_{\mathrm{0}}}{\mathrm{d}z}}\right).\end{array}$$

3.2.2 Boundary and sublayer interface conditions

The values of the coefficients W_j± are determined by the displacement of the inversion layer, the boundary conditions imposed at the interfaces between the sublayers, and the radiation condition. The stratification coefficients are calculated in Fourier components by evaluating the response to a change in ABL height η_t. For each wavenumber k,l, the relation between $\mathit{\eta }\left(z\right)=\widehat{\mathit{\eta }}\left(z\right)\exp \mathit{ı}(kx+ly)$ and W(z) is given by Nappo (2012):

$$\begin{array}{}\text{(24)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{D}}{\mathrm{D}t}}\mathit{\eta }\left(z\right)& =-\mathit{ı}\mathrm{\Omega }\widehat{\mathit{\eta }}\left(z\right)\exp \mathit{ı}(kx+ly)\\ & =W\left(z\right)\exp \mathit{ı}(kx+ly)\end{array}.\end{array}$$

This leads to a boundary condition for W₀(H):

$$\begin{array}{}\text{(25)}& -\mathit{ı}{\mathrm{\Omega }}_{\mathrm{0}}{\widehat{\mathit{\eta }}}_t=W_{\mathrm{0}+}\exp \left(\mathit{ı}{m}_{\mathrm{0}}z\right)+W_{\mathrm{0}-}\exp (-\mathit{ı}{m}_{\mathrm{0}}z).\end{array}$$

This relation is already incorporated in Eq. (20). The displacement η and the pressure perturbation p^′ have to be continuous over the interfaces between sublayers (Gossard and Hooke, 1975; Klemp and Lilly, 1975; Gill, 1982; Baines, 1998; Smith et al., 2002; Vosper, 2004). For each interface at z=H_j between the sublayers j−1 and j, this leads to two conditions that Gossard and Hooke (1975) call the kinematic and dynamic conditions, respectively:

$$\begin{array}{}\text{(26)}& {\displaystyle }{\mathit{\eta }}_{j-\mathrm{1}}\left(H_j\right)={\mathit{\eta }}_j\left(H_j\right)\text{(27)}& {\displaystyle }{p}_{j-\mathrm{1}}^{\prime }\left(H_j\right)=p_j^{\prime }\left(H_j\right).\end{array}$$

Combining Eqs. (24) and (19) allows the kinematic and dynamic conditions to be written in terms of W(z):

$$\begin{array}{}\text{(28)}& {\displaystyle }{\left.{\displaystyle \frac{W_{j-\mathrm{1}}}{{\mathrm{\Omega }}_{j-\mathrm{1}}}}\right|}_{z=H_j}={\left.{\displaystyle \frac{W_j}{{\mathrm{\Omega }}_j}}\right|}_{z=H_j}\text{(29)}& {\displaystyle }\begin{array}{rl}& {\left.\left({\mathrm{\Omega }}_{j-\mathrm{1}}{\displaystyle \frac{\mathrm{d}{W}_{j-\mathrm{1}}}{\mathrm{d}z}}-{\displaystyle \frac{\mathrm{d}{\mathrm{\Omega }}_{j-\mathrm{1}}}{\mathrm{d}z}}{W}_{j-\mathrm{1}}\right)\right|}_{z=H_j}\\ & ={\left.\left({\mathrm{\Omega }}_j{\displaystyle \frac{\mathrm{d}{W}_j}{dz}}-{\displaystyle \frac{\mathrm{d}{\mathrm{\Omega }}_j}{\mathrm{d}z}}{W}_j\right)\right|}_{z=H_j}.\end{array}\end{array}$$

Pütz et al. (2019) argue that if the background velocities change continuously with altitude, the kinematic and dynamic interface conditions should be replaced by conditions ensuring that the vertical velocity is continuous:

$$\begin{array}{}\text{(30)}& {\displaystyle }{\left.W_{j-\mathrm{1}}\right|}_{z=H_j}={\left.W_j\right|}_{z=H_j}\text{(31)}& {\displaystyle }{\left.{\displaystyle \frac{\mathrm{d}{W}_{j-\mathrm{1}}}{\mathrm{d}z}}\right|}_{z=H_j}={\left.{\displaystyle \frac{\mathrm{d}{W}_j}{\mathrm{d}z}}\right|}_{z=H_j}\end{array}$$

since evaluating Ω and its derivative evaluated at the center of each sublayer in Eqs. (28) and (29), as the use of a constant wavenumber implies, results in a discontinuous profile for W. As Pütz et al. (2019) found, these discontinuities cause the piecewise method to converge to a different solution than when solving Eq. (16) with a simple finite-difference solver. However, the physical reasoning behind the kinematic and dynamic boundary conditions is generally valid, indicating that they should always be used. The difference between the setup of Pütz et al. (2019) and the older literature is solved if in Eqs. (28) and (29) not only W_j−1 and W_j are evaluated at z=H_j, but Ω_j−1 and Ω_j as well. While in the piecewise-constant method $m_j^{\mathrm{2}}$ should be taken at the center of each layer and held constant throughout, the background velocities and their derivatives should be evaluated at the interfaces when setting up the interface conditions. If changes in background velocity within layers is taken into account in this way, the kinematic and dynamic matching conditions will automatically simplify to those proposed by Pütz et al. (2019) when appropriate.

Above the highest sublayer, the atmosphere is assumed to be uniformly stratified. This results in the same situation as discussed in Sect. 2.1, with the upper boundary condition determining the sign of m in the propagating regime. For the propagating and evanescent regimes, respectively, the root is chosen so that the radiation condition is satisfied or that the perturbation dies out with altitude.

Combined with the radiation condition applied at height H_n, Eqs. (25), (28), and (29) provide 2(n+1) relations, which is enough to solve for W_j±. As the equations are linear, and the kinematic and dynamic equations at each interface only involve the wave fields in the adjacent sublayers, the system of equations for each wavenumber can be solved as a banded matrix of size 2(n+1) with a bandwidth of 5. Solving the system has a computational complexity of 𝒪(n). The total computation for a single profile of m²(z) with n=100, including the building of the matrix, takes roughly 0.4 s on a personal laptop with 16 GB of RAM and an Intel core i7 2.60 GHz, using scipy routines sped up with the numba package (Lam et al., 2015; Virtanen et al., 2020). To determine the pressure boundary condition in the TLM, the stratification coefficients for all the different wavenumbers have to be calculated, leading to a total complexity of 𝒪(N_xN_yn). Finally, we emphasize that all stratification coefficients can be precomputed as their values do not depend on the solution of the TLM. Once identified, the values of $\widehat{\mathrm{\Phi }}$ can be used in Eq. (10) to set the pressure boundary condition in the TLM.

3.2.3 Inversion and critical layers

If there are inversion layers in the free atmosphere, these can be modeled by discontinuities in $\stackrel{\mathrm{‾}}{\mathit{\theta }}$ corresponding to the inversion strengths $\mathrm{\Delta }\stackrel{\mathrm{‾}}{\mathit{\theta }}$. This can be incorporated in piecewise methods by adding $g^{\prime }=\frac{g\mathrm{\Delta }\stackrel{\mathrm{‾}}{\mathit{\theta }}}{\stackrel{\mathrm{‾}}{\mathit{\theta }}}$ to the right-hand side of Eq. (27) (Baines, 1998).

As the wind can vary, the situation can arise where ${\stackrel{\mathrm{‾}}{u}}_g=-{\stackrel{\mathrm{‾}}{v}}_{g}k/l$ so that the intrinsic frequency Ω=0. In a reference frame moving with the wind speed, the frequency of the waves is then zero, and the mean flow is no longer perturbed. It is clear from Eq. (17) that this corresponds to a singularity. The height at which this occurs is called a critical level, and exactly what happens is hard to predict with linear theory. In general, the wave disappears, and the energy it carried is absorbed by the mean flow (Gossard and Hooke, 1975; Gill, 1982; Baines, 1998). Therefore, critical levels are modeled as fully absorbing sublayers, as is commonly done in the literature (Smith et al., 2002; Wells and Vosper, 2010). This is implemented by applying the radiation condition at interface j when Ω changes sign across interface j+1 or becomes zero in sublayer j+1. Since there is no difference for the flow below the critical level between wave energy being absorbed aloft or just radiating outward indefinitely, this is equivalent for our purposes.

3.3 Verification

To verify the implementation of the piecewise-constant method, it was compared to a second-order finite-difference (FD) code with a central difference scheme on various continuously varying background velocities and buoyancy frequencies. The piecewise method consistently outperformed the FD code, achieving second-order convergence as expected through the proof in Appendix A and having small errors even at coarse grids.

We also reproduced results from Wells and Vosper (2010), who calculated the linear gravity response for an idealized atmosphere to a small 2D ridge, described by

$$\begin{array}{}\text{(32)}& h={\displaystyle \frac{h_r}{{\left(x/L_r\right)}^{\mathrm{2}}+\mathrm{1}}},\end{array}$$

where h_r=10 m and L_r=10 km. These values were chosen so that the effects of non-linearity would remain small (Wells and Vosper, 2010). The vertical velocity perturbation was computed on a 1000 km long domain with 2048 grid points; a 1D version of the grid that is used in Sect. 4; and up to a height of 20 km with a grid spacing of 125 m, which is comparable to what is used in Sect. 4.3. The idealized background velocity and Brunt–Väisäla frequency profiles, as well as the gravity-wave response, are shown in Fig. 2. The contour plots agree fairly well with those obtained by Wells and Vosper (2010), indicating that the piecewise method is correctly implemented and suited to this application. This test case is further discussed in Sect. 4.

Figure 2 Idealized upper-atmosphere velocity (a) and Brunt–Väisäla frequency (b) used by Wells and Vosper (2010) and the vertical velocity field as calculated with the piecewise-constant method (black lines) and as found by Wells and Vosper (2010) (gray lines) (c). The contour interval is 0.01 m s⁻¹, and the solid and dashed lines denote positive and negative values, respectively.

Wells and Vosper (2010) also analyzed the hydrostatic drag response to the same ridge for atmospheres with a two-layer buoyancy-frequency structure. This is a classic test case, and similar setups have been discussed by Gill (1982), Leutbecher (2001), and Teixeira and Argaín (2020), among others. Wells and Vosper (2010) considered a case where the background wind $\stackrel{\mathrm{‾}}{u}\left(z\right)$ is constant and one where it is given by

$$\begin{array}{}\text{(33)}& \stackrel{\mathrm{‾}}{u}\left(z\right)=u_{\mathrm{0}}+u_{\mathrm{1}}\sin \left({\displaystyle \frac{\mathit{\pi }z}{\mathrm{2}{z}_i}}\right),\end{array}$$

where z_i=10 km. In the constant wind case, the drag is computed for a range of u₀. In the case with vertical wind shear, it is computed for a range of u₁, with u₀ kept at 5 m s⁻¹. In reproducing their results, we used a step profile for the Brunt–Väisäla frequency:

$$\begin{array}{}\text{(34)}& N_g\left(z\right)=\left\{\begin{array}{ll}{N}_{\mathrm{1}},& z<z_i\\ N_{\mathrm{2}},& z>z_i\end{array},\right.\end{array}$$

where N₁=0.01, and $N_{\mathrm{2}}=\mathrm{0.02}{\mathrm{s}}^{-\mathrm{1}}$. Our results, obtained on the same grid as above with our method adapted to the hydrostatic regime, and those obtained by Wells and Vosper (2010) and Leutbecher (2001) for the constant wind case are shown in Fig. 3. Again, the good agreement indicates that the multilayer method performs well. An overview of all the upper-atmospheric profiles used for verification is given in Table 1.

Wells and Vosper (2010)

Table 1 An overview of the different upper-atmospheric flow profiles used for verification. The upper atmosphere set up by Wells and Vosper (2010) is also used in Sect. 4.1 and 4.2.

Figure 3 Mountain wave drag on a small ridge with a two-layer Brunt–Väisäla frequency profile and constant background wind (a) and vertical wind shear (b), normalized by the drag for constant background profiles. The black lines with squares show our results, in which the wave drag is normalized with the drag for constant background wind u₀ and stratification strength N₁. The blue lines with circles show the results of Wells and Vosper (2010), and the red line with triangles in the left figure shows the results of Leutbecher (2001).

4 Effects of vertically varying wind and stability

By using the piecewise-constant method developed in Sect. 3.2 to evaluate the stratification coefficients, the TLM can now take the variation in the stratification strength and wind speed with altitude into account. The impact of this non-uniformity on the interaction between gravity waves and wind farms is now analyzed in three ways. In Sect. 4.1, the impact of atmospheric profile variations on the stratification coefficients is analyzed by comparing $\widehat{\mathrm{\Phi }}$ for the idealized atmosphere used in Sect. 3.3 to the $\widehat{\mathrm{\Phi }}$ for a uniform upper atmosphere. The physical phenomenon causing the differences is identified as internal gravity-wave interference.

We further investigate how this influences the interaction between wind farms and gravity waves. To this end, Sect. 4.2 discusses an example case of ABL flow with a uniform upper atmosphere from Allaerts and Meyers (2019). By combining this case with the atmospheric profiles analyzed in previous sections, we investigate the effect of vertical variations in the atmospheric conditions on wind-farm–gravity-wave interaction. Finally, Sect. 4.3 presents a case study that estimates the overall impact of such variations by simulating 1 year of operation of the Belgian–Dutch offshore wind-farm cluster.

4.1 Stratification coefficients

To determine the effects of vertical non-uniformity on the pressure feedback of internal gravity waves, the stratification coefficients are calculated for the upper atmosphere used by Wells and Vosper (2010), shown in Fig. 2. The same grid as in Sect. 3.3 is used. The coefficients $\widehat{\mathrm{\Phi }}$ are calculated for both the original non-uniform and a uniform atmosphere. The latter was obtained by height-averaging the profiles of the velocity and the Brunt–Väisäla frequency, resulting in $U=\mathrm{20.1}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ and $N_g=\mathrm{0.0113}{\mathrm{s}}^{-\mathrm{1}}$. Figure 4 shows the results. Only $k<\mathrm{1}{\mathrm{km}}^{-\mathrm{1}}$ is shown so that the important details are clearly visible. Since this is a 2D case, the transversal wavenumber l is zero. Comparing the two profiles for $\mathrm{|}\mathrm{|}\widehat{\mathrm{\Phi }}\mathrm{|}\mathrm{|}$, it is clear that the vertical non-uniformity has a large impact on the stratification coefficients, with two peaks appearing at k≈0.36 and $k\approx \mathrm{0.72}{\mathrm{km}}^{-\mathrm{1}}$. We will now show that these changes in the profile of $\mathrm{|}\mathrm{|}\widehat{\mathrm{\Phi }}\mathrm{|}\mathrm{|}$ are caused by resonance in the free atmosphere.

Figure 4 The stratification coefficients for a uniform atmosphere with U=20.1 m s⁻¹ and N_g=0.0113 s⁻¹ (a), the stratification coefficients (b) and the gravity-wave resonance parameter (see Eq. 35) (c) for the vertical non-uniform atmosphere of Wells and Vosper (2010) used in Sect. 3.3. Only $k<\mathrm{1}{\mathrm{km}}^{-\mathrm{1}}$ is shown so that the important details are clearly visible. Since this is a 2D case, the transversal wavenumber l is zero. Above $k=\mathrm{0.408}{\mathrm{km}}^{-\mathrm{1}}$, the gravity waves become evanescent within some sublayers, leading to oscillations in A that do not correspond to resonant behavior.

4.1.1 Internal wave resonance

While in a uniform atmosphere the wind farm can only trigger waves with an upwards group velocity, changes in the stratification and wind speed can cause waves to reflect. This allows both up- and downgoing waves to propagate throughout the atmosphere. As up- and downgoing internal gravity waves pass through each other, they interfere, potentially causing resonance. The resulting large wave amplitudes can drastically affect the pressure feedback of the waves (Teixeira, 2014). The ratio A between the mean wave energy at a given wavenumber for atmospheric profile with and without vertical variations is an effective measure for gravity-wave resonance as long as the waves are in the propagating regime (Gill, 1982). It is easily evaluated numerically using the following expression:

$$\begin{array}{}\text{(35)}& A=\underset{H}{\overset{H_n}{\int }}{e}_{HD}\left(z\right)\mathrm{d}z\cdot {\left[{\displaystyle }\underset{H}{\overset{H_n}{\int }}{e}_U\left(z\right)\mathrm{d}z\right]}^{-\mathrm{1}},\end{array}$$

where e is the sum of the kinetic and potential energy densities (Gill, 1982):

$$\begin{array}{}\text{(36)}& e={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\stackrel{\mathrm{‾}}{\mathit{\rho }}\left(\stackrel{\mathrm{‾}}{{u^{\prime }}^{\mathrm{2}}+{v^{\prime }}^{\mathrm{2}}+{w^{\prime }}^{\mathrm{2}}}\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{g}^{\mathrm{2}}\stackrel{\mathrm{‾}}{{{\mathit{\rho }}^{\prime }}^{\mathrm{2}}}/\stackrel{\mathrm{‾}}{\mathit{\rho }}{N}_g^{\mathrm{2}}.\end{array}$$

For a sublayer in a piecewise-constant method, the integrals are straightforward to determine analytically, allowing A to be evaluated.

Figure 4 shows A for the atmosphere of Wells and Vosper (2010). It is clear that $\mathrm{|}\mathrm{|}\widehat{\mathrm{\Phi }}\mathrm{|}\mathrm{|}$ and A have a similar profile, with the same peaks occurring at $k\approx \mathrm{0.36}{\mathrm{km}}^{-\mathrm{1}}$ and $k\approx \mathrm{0.72}{\mathrm{km}}^{-\mathrm{1}}$. Above $k=\mathrm{0.408}{\mathrm{km}}^{-\mathrm{1}}$, the gravity waves become evanescent within some sublayers, leading to oscillations in A that do not correspond to resonant behavior. Despite this, it is clear from the figure that the profiles of $\mathrm{|}\mathrm{|}\widehat{\mathrm{\Phi }}\mathrm{|}\mathrm{|}$ and A follow the same pattern, showing that the pressure feedback is largely determined by constructive and destructive interference of the internal gravity waves.

4.2 Gravity-wave ABL interaction

We investigate how variations in the wind and stability change the effects wind-farm operation has on the ABL flow by revising an example case used by Allaerts and Meyers (2019). By combining this case with the non-uniform upper atmosphere used in the previous sections, we identify the mesoscale flow perturbations triggered by the wind farm. One of these changes caused by the variations in the atmospheric conditions – the increase in the appearance and strength of resonant lee waves – is further analyzed.

4.2.1 Example cases

To analyze how the changes in the stratification coefficients impact the flows around wind farms, a flow case discussed earlier by Allaerts and Meyers (2019) is analyzed. This example case is simulated once with uniform upper atmospheres and once with the upper atmosphere discussed in Sects. 3.3 and 4.1.

The example case is set up to have a Froude number of Fr=1.1, with the Froude number given by Allaerts and Meyers (2019),

$$\begin{array}{}\text{(37)}& \mathit{Fr}={\displaystyle \frac{{\stackrel{\mathrm{‾}}{u}}_B}{\sqrt{g^{\prime }H}}},\end{array}$$

where ${\stackrel{\mathrm{‾}}{u}}_B$ is a velocity scale for the ABL (Allaerts and Meyers, 2019),

$$\begin{array}{}\text{(38)}& {\stackrel{\mathrm{‾}}{u}}_B={\left({\displaystyle \frac{H_{\mathrm{1}}}{H}}{\displaystyle \frac{\mathrm{1}}{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{1}}^{\mathrm{2}}}}+{\displaystyle \frac{H_{\mathrm{2}}}{H}}{\displaystyle \frac{\mathrm{1}}{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{2}}^{\mathrm{2}}}}\right)}^{-\mathrm{1}/\mathrm{2}}.\end{array}$$

The Froude number represents the ratio of the advection terms in the momentum balance to the pressure gradient generated by the inversion waves (Durran, 1990). It is also a measure for the ratio of the wind speed to the velocity of inversion waves and therefore determines whether or not inversion waves can travel upstream from a stationary perturbation. This leads to the classification of flows as being either subcritical Fr<1, when inversion waves can affect the upstream flow, or supercritical Fr>1, when they cannot (Smith, 2010). Allaerts and Meyers (2019) analyzed two cases: one is subcritical, with Fr=0.9, while the other is supercritical, with Fr=1.1. Here, we focus on the supercritical case.

The example case is set up as follows. The mean flow in the two layers of the ABL is based on the analytical boundary layer model of Nieuwstadt (1983), with a cubic eddy viscosity profile ${\mathit{\nu }}_{\mathit{\tau }}=\mathit{\kappa }{u}_{*}z(\mathrm{1}-z/H{)}^{\mathrm{2}}$ and the Von Kármán constant κ=0.41 (Allaerts and Meyers, 2019). This model finds a velocity profile given a non-dimensional surface roughness length ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}}=z_{\mathrm{0}}/H$ and non-dimensional boundary layer height $h_{*}=H{f}_c/u_{*}$. Combined with the stratification parameters g^′ and N_g, this gives all the inputs required for the TLM. Allaerts and Meyers (2019) presume a conventionally neutral boundary layer, determining these stratification parameters with the non-dimensional groups $N_g/f_{\mathrm{c}}$ and $g^{\prime }H/A{u}_{*}^{\mathrm{2}}$, where A=500 is an empirical constant (Csanady, 1974; Tjernstrom and Smedman, 1993). Since we want to analyze the impact vertical variations in the upper-atmospheric profile have on the flow, the original inputs used by Allaerts and Meyers (2019) are altered so that the upper atmosphere corresponds to the one analyzed in Sect. 4.1. The correct geostrophic wind and Brunt–Väisäla frequency are obtained by modifying u_* and $N_g/f_{\mathrm{c}}$, respectively, and $g^{\prime }H/A{u}_{*}^{\mathrm{2}}$ is changed so that the final result still has a Froude number of 1.1. H is increased so that h_* and ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}}$ remain the same. The inputs are summarized in Table 2. Finally, the flow in the lower layer is aligned with the x direction.

Table 3 gives an overview of the wind-farm configuration used by Allaerts and Meyers (2019), which is chosen to be comparable to the Belgian–Dutch wind-farm cluster in area and installed capacity. The turbines are placed in a staggered pattern with respect to the x direction, and the relative spacing is equal in the x and y directions. The simulations were performed on a 1000 km by 400 km grid, with 2000 by 800 grid points, which is the same as the grid spacing used in Sect. 3.3. For the vertically non-uniform simulation, the same vertical grid as in Sect. 3.3 and 4.1 is used.

Allaerts and Meyers (2019)

Table 2 Flow parameters of the flow case based on the one used by Allaerts and Meyers (2019). The atmospheric state corresponds to a P_N of 1.1 and Fr of 1.1. For the vertically non-uniform case, the upper-atmospheric profiles of Wells and Vosper (2010), shown in Fig. 2, were used.

Table 3 Wind-farm configuration of the reference flow cases, as analyzed by Allaerts and Meyers (2019).

The results of the uniform and non-uniform simulations are shown in Fig. 5. When comparing the results, it is clear that the mesoscale disturbances are much larger in the non-uniform case, with a stronger reduction in the farm inflow velocity. The farm's effects also spread out in the transversal directions as V-shaped patterns appear over large distances. Finally, the wind farm seems to trigger strong lee waves in its wake. These waves also appear in the reference case, although they are weaker there.

Figure 5 Planform view of the inversion layer displacement η_t (a, d), pressure perturbation p^′ (b, e), and velocity reduction in the lower layer u^′ (c, f) in the reference uniform (a–c) and vertically varying (d–f) flow cases. The wind-farm region is indicated by the black rectangles.

4.2.2 Resonant lee waves in the ABL

We further analyze the lee waves that appear in Fig. 5 and show that they are associated in this case with internal wave reflection. We follow the ideas developed by Allaerts and Meyers (2019) by analyzing the equation for the total displacement η_t. When the wind in the ABL is aligned with the x axis $\left({\stackrel{\mathrm{‾}}{v}}_{\mathrm{1}}={\stackrel{\mathrm{‾}}{v}}_{\mathrm{2}}=\mathrm{0}\right)$, that equation corresponds to Allaerts and Meyers (2019):

$$\begin{array}{}\text{(39)}& \begin{array}{rl}& \left(-\mathrm{1}+F{r}^{-\mathrm{2}}+P_N^{-\mathrm{1}}{\displaystyle \frac{\mathrm{\Phi }}{G{N}_g}}\right)*{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\eta }}_t}{\partial x^{\mathrm{2}}}}+\left(F{r}^{-\mathrm{2}}+P_N^{-\mathrm{1}}{\displaystyle \frac{\mathrm{\Phi }}{G{N}_g}}\right)\\ & *{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\eta }}_t}{\partial y^{\mathrm{2}}}}={\displaystyle \frac{H_{\mathrm{1}}}{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{1}}^{\mathrm{2}}}}\mathrm{\nabla }\cdot {\mathrm{RHS}}_{\mathrm{1}}+{\displaystyle \frac{H_{\mathrm{2}}}{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{2}}^{\mathrm{2}}}}\mathrm{\nabla }\cdot {\mathrm{RHS}}_{\mathrm{2}},\end{array}\end{array}$$

where G is the geostrophic wind velocity. Furthermore, RHS₁ and RHS₂ are the right-hand sides of the Eqs. (6) and (7), respectively, and only depend on η_t through the turbulent viscosity, which is a relatively weak effect (Allaerts and Meyers, 2019). As discussed before, the Froude number Fr indicates the strength of the inversion waves. The second non-dimensional group P_N governs the pressure induced by internal gravity waves and is given by

$$\begin{array}{}\text{(40)}& P_N={\displaystyle \frac{{\stackrel{\mathrm{‾}}{u}}_B^{\mathrm{2}}}{G{N}_{g}H}}.\end{array}$$

It reflects the fact that the pressure induced by a vertical displacement scales linearly with GN_g in the hydrostatic regime. It is then clear that the left-hand side of Eq. (39) represents the forcing induced by flow advection $\left(-\frac{{\partial }^{\mathrm{2}}{\mathit{\eta }}_t}{\partial x^{\mathrm{2}}}\right)$ and the corresponding gravity waves $\left(\left({\mathit{Fr}}^{-\mathrm{2}}+P_N^{-\mathrm{1}}\mathrm{\Phi }/G{N}_g*\right){\mathrm{\nabla }}_H^{\mathrm{2}}{\mathit{\eta }}_t\right)$. These are balanced by the right-hand side of the equation, which represents the other terms in the momentum equations. If the advection terms and the pressure contributions from the gravity waves balance each other, the left-hand side of Eq. (39) becomes zero. This corresponds to a resonant state as η_t can be non-zero without external forcing, which explains the lee waves that appear in Fig. 5.

The left-hand side of Eq. (39) is easily expressed in Fourier components, leading to the definition of the two-dimensional lee-wave resonance parameter R:

$$\begin{array}{}\text{(41)}& R={\left(\cos \mathit{\lambda }\right)}^{\mathrm{2}}-{\mathit{Fr}}^{-\mathrm{2}}-{\displaystyle \frac{H}{{\stackrel{\mathrm{‾}}{u}}_B^{\mathrm{2}}}}\widehat{\mathrm{\Phi }},\end{array}$$

where λ is the angle the horizontal wave vector makes with the k axis so that ${\left(\cos \mathit{\lambda }\right)}^{\mathrm{2}}=k^{\mathrm{2}}/(k^{\mathrm{2}}+l^{\mathrm{2}})$. The parameter R is a non-dimensional parameter indicating the flow's resistance to the occurrence of two-dimensional lee-wave resonance. If R=0, the pressure perturbations induced by η_t≠0 and the accompanying gravity waves balance the advection terms, and resonant lee waves appear. Equation (41), in combination with Eq. (4), shows that for uniform upper atmospheres, this type of resonance can only take place if the internal gravity waves are evanescent as $\widehat{\mathrm{\Phi }}$ has to be real for R to be zero. Physically, this corresponds to propagating waves not being able to trap the perturbation energy below the capping inversion (Vosper, 2004). In contrast, the wave reflection in vertically non-uniform atmospheres can cause part of the wave energy to be trapped by being reflected back. This corresponds to the stratification coefficients as given by Eq. (23) potentially having a real component, even in the propagating regime. The additional constraint found by Allaerts and Meyers (2019), that the flow has to be subcritical, is also not necessary when the waves can be reflected.

In order to apply this theory to the case from Sect. 4.2.1, Fig. 6 shows the wavenumber spectra for η_t and R. Figure 6 shows that low values of ||R|| correspond to high values of $\mathrm{|}\mathrm{|}{\widehat{\mathit{\eta }}}_t\mathrm{|}\mathrm{|}$, indicating resonant lee waves in the ABL. The variations in wind and stability in the upper atmosphere can cause ||R|| to decrease several orders of magnitude in the propagating wave regime. This happens when the internal gravity-wave interference is destructive, which lowers the vertical energy flux of the waves considerably. This keeps the perturbation energy contained in the lower atmosphere, leading to low ||R||, which for uniform atmospheres only happens in the evanescent wave regime (Vosper, 2004; Allaerts and Meyers, 2019). In contrast, constructive interference can lead to very large values for $\mathrm{|}\mathrm{|}\widehat{\mathrm{\Phi }}\mathrm{|}\mathrm{|}$, and thus for ||R||. Furthermore, when ||R||≫0, a large imbalance exists between the advection and pressure forces in the ABL. Changes in ABL height of the wavelengths at which this occurs cannot exist in an equilibrium system without external forcing. This is clearly visible in Fig. 6 as high values of ||R|| correspond to low values of $\mathrm{|}\mathrm{|}{\widehat{\mathit{\eta }}}_t\mathrm{|}\mathrm{|}$. Internal gravity-wave resonance thus prohibits large displacements of the inversion layer at the wavelengths at which it occurs by creating large pressure gradients that cannot be counteracted by the flow acceleration.

Figure 6 $\mathrm{|}\mathrm{|}{\widehat{\mathit{\eta }}}_t\mathrm{|}\mathrm{|}/H$ (a, c) and ||R|| (b, d) for the reference flow case with a uniform (a, b) and a non-uniform free atmosphere (c, d). Only the wavenumbers $\mathrm{|}\mathrm{|}k\mathrm{|}\mathrm{|}\le \mathrm{1.26}{\mathrm{km}}^{-\mathrm{1}}$ and $\mathrm{|}\mathrm{|}l\mathrm{|}\mathrm{|}\le \mathrm{1.26}{\mathrm{km}}^{-\mathrm{1}}$ are shown so that the important details are clearly visible. The left column shows parts of the Fourier transforms of the left column of Fig. 5.

While the above analysis explains how vertical variations in the atmospheric profile change the interaction between internal gravity waves and the ABL flow, it does not offer insight on how to predict the resulting impact on wind-farm performance. It is not clear what parameters could describe this. Extensive research has been done on internal gravity-wave resonance, with most of it focusing on flow around topographies. However, this has to be used with caution as the overall flow is then analyzed in the context of mountain wave drag, where the height displacement is given by the shape of the terrain (Teixeira, 2014). In contrast, the wind-farm forcing leads to this displacement through its interaction with the ABL. Additionally, it is itself influenced by the mesoscale effects it triggers, leading to an additional feedback loop. Parameters successfully describing internal wave resonance may therefore not be able to predict how vertical non-uniformity will impact the interaction with wind farms.

4.3 Overall impact

To determine the impact of varying wind speeds and stability on wind-farm energy production, we follow the approach of Allaerts et al. (2018) by simulating 1 year of wind-farm operation of the Belgian–Dutch wind-farm cluster. This analysis is performed with both uniform and non-uniform upper atmospheres, and the two are compared. The TLM input is based on ERA5 reanalysis data of the year 2016 from Hersbach et al. (2018), which are available at hourly frequency, resulting in 8784 flow cases. Like Allaerts et al. (2018), we use data from the grid point nearest to the wind-farm cluster, at 51.6^∘ N, 3.0^∘ E, and the same approach in determining the TLM input from the atmospheric data. We further assume that all turbines are DTU 10 MW reference turbines, a commonly analyzed model for which the power curve is readily available (Bortolotti et al., 2019), and place them within the same trapezoidal area as Allaerts et al. (2018). To ensure that we obtain the same total installed capacity of 3.8 GW, we scale the number of turbines and the turbine spacing from 475 and 7.36 to 380 and 8.23, respectively. The thrust coefficient of the turbines is set based on the height-averaged background velocity in the wind-farm layer of the TLM.

The simulations were performed on a 2000 by 2000 grid, on a 1000 km by 1000 km domain, the same grid density as in the previous sections. Variations in the atmospheric profiles were taken into account up to the tropopause. Within the troposphere, the potential temperature and velocity profiles were modeled as first- and third-order splines, respectively, through the ERA5 data points, around which the sublayers were spaced as well, resulting in 30 to 50 sublayers for each case. The tropopause altitude and the stratification in the stratosphere were determined with a two-line, piecewise-linear regression fit on the temperature profiles between the ABL height and 15 km. The velocity in the stratosphere was then determined by height-averaging the profile from the tropopause up to 15 km. As an example, Fig. 7 shows the profiles of θ, u, and v for the upper atmosphere at 00:00 GMT on 1 May 2016.

Figure 7 The atmospheric profiles of potential temperature (a), velocity in the x (b) and y direction (c). The blue crosses are the ERA5 data, and the black lines are the inputs to the TLM, showing the piecewise-constant approximation made in the discretization.

Figure 8 The maximum capping inversion displacement (a) and the inflow velocity perturbation (b) for all analyzed cases with both uniform (x axis) and non-uniform (y axis) upper atmospheres.

From the 8746 cases, we only use those where the atmosphere is statically stable at every altitude in the free atmosphere. Additionally, the cases without capping inversion, (cf. earlier discussion following Eq. 10 in Sect. 2.2) or with a capping inversion situated lower than twice the turbine hub height were left out, leaving 3890 cases. This filtering was necessary as for the removed cases the assumptions made in the derivation of the model are not valid, as discussed in Sect. 2.2. The results of the simulations are shown in Fig. 8 and summarized in Table 4. On average, the difference between the results with non-uniform and uniform upper atmospheres is small. Despite this however, the impact on individual cases is often significant. In 16.6 % of the analyzed cases, the difference between the uniform and the non-uniform inflow perturbation was more than 30 % of the uniform perturbation case. We therefore conclude that vertical variations in the upper-atmospheric profiles are important to take into account when analyzing individual flow cases.

Table 4 Average perturbations over all the analyzed flow cases for both uniform and non-uniform upper atmospheres.

Figure 9 The difference between the maximum capping inversion displacements in the non-uniform and uniform simulations as normalized by the ABL height (a) and the displacement in the uniform simulation (b), plotted against the parameters Fr and P_N. The lack of a clear trend in both figures indicates that Fr and P_N are not related to the impact of the vertical non-uniformity.

Figure 9 shows that the differences between the simulations seem to be independent from parameters that were good predictors of the TLM's behavior in previous studies, such as Fr and P_N. We also investigated parameters that were found to correlate well with mountain wave drag in vertically non-uniform atmospheres, such as c₁ as found by Klemp and Lilly (1975) (Eq. 16), but could not identify a meaningful correlation in our case. Directional shear effects, such as those investigated by Teixeira et al. (2008), might explain the discrepancy.

5 Conclusions

The goal of this study was to extend the applicability of a wind-farm gravity-wave model to vertically non-uniform free atmospheres. This was done by changing the expressions for the stratification coefficients $\widehat{\mathrm{\Phi }}$ to results derived from the internal wave equation for general stratified flows. By applying the well-known piecewise method with large numbers of sublayers, general stratification and velocity profiles can be incorporated into the model.

The effects of the variations in background wind and stability were studied by analyzing how free-atmospheric wave reflection influences the wave pressure feedback, the ABL flow, and overall wind-farm performance. Firstly, the stratification coefficients for the idealized atmosphere used by Wells and Vosper (2010) were compared to those for a uniform atmosphere based on it. The differences were found to be caused by constructive and destructive interference of the internal gravity waves in the free atmosphere. In a second step, this vertically non-uniform atmosphere was combined with a flow case used by Allaerts and Meyers (2019). An analysis of the appearance of resonant lee waves led to a qualitative understanding of the vertical non-uniformity's effects on the interaction between gravity waves and the ABL flow. Due to destructive internal wave interference, resonant lee waves can appear. On the other hand, internal gravity-wave resonance dampens inversion layer displacement for the wavelengths at which it occurs.

Finally, the extended TLM was used to simulate 1 year of operation of the Belgian–Dutch wind-farm cluster, repeating a similar analysis by Allaerts et al. (2018), in order to determine overall impact of variations in the atmospheric profiles on the interaction between the ABL flow and the gravity waves to be determined. While this impact was found to be small when averaged out over all flow cases, for individual flow cases it is often of the same order of magnitude as the total flow perturbation. In 16.6 % of the analyzed flows, the variations in the atmospheric profiles caused a relative difference in the upstream velocity reduction of more than 30 %. It is unclear how the effect of this non-uniformity for individual flow cases can be predicted since it does not simply scale with Fr or P_N.

The results of this study show that vertical atmospheric non-uniformity could play a major role in the interaction between wind farms and gravity waves. This suggests that variations with altitude of the free atmosphere's wind and stability should be taken into account when simulating wind-farm operation in specific atmospheric conditions and may be important for the optimization of turbine control in the future (see, e.g., Lanzilao and Meyers, 2021b). In the future, we foresee further improvements of the TLM, among others, including hydrodynamic effects in the boundary layer and upgrading the wake model to include the improved wake merging model by Lanzilao and Meyers (2021a). Next to that, we plan further validation against detailed large-eddy simulations (similar to Allaerts and Meyers, 2019) and data from operational wind farms.

Appendix A: Convergence analysis

Adding more sublayers leads to a better approximation of the actual profiles of atmospheric variables. This would then also result in a better approximation of W(z). To estimate the rate of convergence, z_* is defined as the altitude where $\frac{\partial \stackrel{\mathrm{̃}}{W}}{\partial z}\left(z_{*}\right)=\frac{\partial W}{\partial z}\left(z_{*}\right)$. If n is sufficiently large, such an altitude exists in each sublayer. The error $E=W-\stackrel{\mathrm{̃}}{W}$ can then be written as a Taylor expansion around z_*:

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}E\left(z\right)& =W\left(z\right)-\stackrel{\mathrm{̃}}{W}\left(z\right)\\ & =W\left(z_{*}\right)+{\displaystyle \frac{\partial W}{\partial z}}\left(z_{*}\right)\mathrm{\Delta }z+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{d^{\mathrm{2}}W}{\mathrm{d}{z}^{\mathrm{2}}}}\left(z_{*}\right)\mathrm{\Delta }{z}^{\mathrm{2}}\\ & -\stackrel{\mathrm{̃}}{W}\left(z_{*}\right)-{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{̃}}{W}}{\mathrm{d}z}}\left(z_{*}\right)\mathrm{\Delta }z-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{d^{\mathrm{2}}\stackrel{\mathrm{̃}}{W}}{\mathrm{d}{z}^{\mathrm{2}}}}\left(z_{*}\right)\mathrm{\Delta }{z}^{\mathrm{2}}+\mathrm{\dots }\\ & =W\left(z_{*}\right)\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{m}^{\mathrm{2}}\left(z_{*}\right)\mathrm{\Delta }{z}^{\mathrm{2}}\right)\\ & -\stackrel{\mathrm{̃}}{W}\left(z_{*}\right)\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\stackrel{\mathrm{̃}}{m}}^{\mathrm{2}}\left(z_{*}\right)\mathrm{\Delta }{z}^{\mathrm{2}}\right)+\mathrm{\dots }\end{array}\end{array}$$

Using the above result, the change in error over the jth sublayer Δ_jE can be written as

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{\mathrm{\Delta }}_{j}E& =E\left(H_{j+\mathrm{1}}\right)-E\left(H_j\right)\\ & =W\left(z_{*}\right)\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{m}^{\mathrm{2}}\left(z_{*}\right)\mathrm{\Delta }{z}_{j+\mathrm{1}}^{\mathrm{2}}\right)\\ & -\stackrel{\mathrm{̃}}{W}\left(z_{*}\right)\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\stackrel{\mathrm{̃}}{m}}^{\mathrm{2}}\left(z_{*}\right)\mathrm{\Delta }{z}_{j+\mathrm{1}}^{\mathrm{2}}\right)\\ & -W\left(z_{*}\right)\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{m}^{\mathrm{2}}\left(z_{*}\right)\mathrm{\Delta }{z}_j^{\mathrm{2}}\right)\\ & +\stackrel{\mathrm{̃}}{W}\left(z_{*}\right)\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\stackrel{\mathrm{̃}}{m}}^{\mathrm{2}}\left(z_{*}\right)\mathrm{\Delta }{z}_j^{\mathrm{2}}\right)+\mathrm{\dots }\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{m}^{\mathrm{2}}\left(z_{*}\right)\left(\mathrm{\Delta }{z}_j^{\mathrm{2}}-\mathrm{\Delta }{z}_{j+\mathrm{1}}^{\mathrm{2}}\right)\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\stackrel{\mathrm{̃}}{m}}^{\mathrm{2}}\left(z_{*}\right)\left(\mathrm{\Delta }{z}_j^{\mathrm{2}}-\mathrm{\Delta }{z}_{j+\mathrm{1}}^{\mathrm{2}}\right)+\mathrm{\dots }\end{array}\end{array}$$

Because $\mathrm{\Delta }{z}_{j+\mathrm{1}}+\mathrm{\Delta }{z}_j={\mathrm{\Delta }}_{j}H$, with Δ_jH the thickness of the jth sublayer, a maximum value for |ΔE_j| is given by

$$\begin{array}{}\text{(A3)}& \mathrm{|}{\mathrm{\Delta }}_{j}E\mathrm{|}\le {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left|m^{\mathrm{2}}\left(z_{*}\right)-{\stackrel{\mathrm{̃}}{m}}^{\mathrm{2}}\left(z_{*}\right)\right|{\mathrm{\Delta }}_j{H}^{\mathrm{2}}+\mathrm{\dots }\end{array}$$

If $m^{\mathrm{2}}={\stackrel{\mathrm{̃}}{m}}^{\mathrm{2}}$ somewhere in the sublayer, substituting a Taylor expansion for m² then directly leads to

$$\begin{array}{}\text{(A4)}& \mathrm{|}{\mathrm{\Delta }}_{j}E\mathrm{|}\le {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left|{\displaystyle \frac{\mathrm{d}{m}^{\mathrm{2}}}{\mathrm{d}z}}\left(z_{*}\right)\right|{\mathrm{\Delta }}_j{H}^{\mathrm{3}}+\mathrm{\dots }\end{array}$$

From this analysis, it is clear that the piecewise-constant method has the best result when the distance to z_* is minimized for all z. Therefore, when approximating a general continuously varying atmosphere, sublayers should be evenly spaced, and the atmospheric state should be evaluated halfway through each sublayer when calculating ${\stackrel{\mathrm{̃}}{m}}_j^{\mathrm{2}}$. In that case, the maximum error can be expected to scale with n|Δ_jE| and ${\mathrm{\Delta }}_{j}H\sim n^{-\mathrm{1}}$, resulting in a second-order rate of convergence, as also found with a different derivation by Pütz et al. (2019). This is confirmed by comparison with a finite-difference solver.

It is notable that the above derivation is not limited to piecewise-constant methods but is valid for general piecewise methods. Therefore, as long as $\stackrel{\mathrm{̃}}{W}\left(z\right)$ is not set up so that $\stackrel{\mathrm{̃}}{W}\left(z_{*}\right)=W\left(z_{*}\right)$, the convergence rate will remain second-order. For example, we also developed a piecewise-linear methods using Airy functions instead of exponential functions that did not outperform the piecewise-constant one. This is because when linearly approximating a general function m²(z), $\stackrel{\mathrm{̃}}{W}\left(z_{*}\right)$ and W(z_*) will still not coincide. As a result, there is no improvement for general atmospheric profiles.