# Wind turbine sound propagation: Comparison of a linearized Euler equations model with parabolic equation methods

Jules Colas, Ariane Emmanuelli, Didier Dragna, Philippe Blanc-Benon, Benjamin Cotté, Richard J. A. M. Stevens

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Jules Colas, Ariane Emmanuelli, Didier Dragna, Philippe Blanc-Benon, Benjamin Cotté, et al.. Wind turbine sound propagation: Comparison of a linearized Euler equations model with parabolic equation methods. Journal of the Acoustical Society of America, 2023, 154 (3), pp.1413-1426. 10.1121/10.0020834 . hal-04209955v2

HAL Id: hal-04209955

https://hal.science/hal-04209955v2

Submitted on 14 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Wind turbine sound propagation: Comparison of a linearized Euler equations model with parabolic equation methods

Jules Colas,<fnref n="1" /> Ariane Emmanuelli,<fnref n="1" /> Didier Dragna,<fnref n="1" /> Phillippe Blanc-Benon,<fnref n="2" /> Benjamin Cotté,<fnref n="3" /> and Richard Stevens<fnref n="4" />

Noise generated by wind turbines is significantly impacted by its propagation in the atmosphere. Hence, for annoyance issues an accurate prediction of sound propagation is critical to determine noise levels around wind turbines. This study presents a method to predict wind turbine sound propagation based on linearized Euler equations. We compare this approach to the parabolic equation method, which is widely used since it captures the influence of atmospheric refraction, ground reflection, and sound scattering at a low computational cost. Using the linearized Euler equations is more computationally demanding but can reproduce more physical effects as fewer assumptions are made. An additional benefit of the linearized Euler equations is that they provide a time-domain solution. To compare both approaches, we simulate sound propagation in two distinct scenarios. In the first scenario, a wind turbine is situated on flat terrain; in the second, a turbine is situated on a hilltop. The results show that both methods provide similar noise predictions in the two scenarios. We find that while some differences in the propagation results are observed in the second case, the final predictions for a broadband extended source are similar between the two methods.

©2023 Acoustical Society of America. [http://dx.doi.org/10.1121/10.0020834]

[Editor: Vladimir E. Ostashev]

Pages: 1–15

## I. INTRODUCTION

The noise produced by wind turbines is one of the main constraints for the installation of new wind farms. It can also entail curtailment plans when annoyance issues emerge after the construction of the farm leading to energy production loss up to 70 % during the night (Du-mortier et al., 2015). Therefore, accurate prediction tools are needed to assess the wind farm noise during the development and operational phases. The modeling of the aerodynamic noise sources and the propagation of sound in the atmosphere are two critical fields of study for the prediction of wind turbine noise. Several models for the noise emitted by a wind turbine have been developed in recent years (Barlas et al., 2017a; Cotté, 2019). These models aim at capturing both the mean sound pressure level (SPL) emitted by the source, and the amplitude modulation (AM) induced by the rotation of the blades. The unsteady nature of the source is considered to be one of the main annoyance causes as the global SPL is usually quite low (Hansen et al., 2019). To predict the SPL field around the turbine, outdoor sound propagation models need to consider ground effects, atmospheric

absorption, and refraction induced by the variability of wind and temperature within the atmospheric boundary layer (ABL). The effect of the ABL flow on wind turbine sound propagation has been extensively studied (Barlas et al., 2018; Heimann and Englberger, 2018). The evolution of the temperature and wind speed gradients during the day modifies the sound propagation and, hence, the SPL around the wind turbine. Furthermore, the presence of the wind turbine creates a wake that acts as a waveguide and tends to increase the SPL at a specific location downwind of the turbine (Barlas et al., 2017b). Finally, topography was also shown to play a significant role in wind turbine sound propagation (Heimann et al., 2018; Sessarego and Shen, 2020; Shen et al., 2019).

Various numerical methods have been used to calculate wind turbine sound propagation. Engineering tools based on simple empirical models are more suited for operational purposes. Methods based on ray acoustics can consider atmospheric refraction (Heimann and Englberger, 2018; Prospathopoulos and Voutsinas, 2007). However, they are usually less precise in terms of SPL and rely on a high-frequency approximation. The parabolic equation (PE) methods have been used extensively both

<footnote n="1">¹Univ Lyon, Ecole Centrale de Lyon, CNRS, Univ Claude Bernard Lyon 1, INSA Lyon, LMFA, UMR5509, 69134 Ecully CEDEX, France</footnote>

<footnote n="2">²Univ Lyon, CNRS, Ecole Centrale de Lyon, INSA Lyon, Univ Claude Bernard Lyon 1, LMFA, UMR5509, 69134 Ecully CEDEX, France</footnote>

<footnote n="3">³Institute of Mechanical Sciences and Industrial Applications (IMSIA), ENSTA Paris, CNRS, CEA, EDF, Institut Polytechnique de Paris</footnote>

<footnote n="4">⁴Physics of Fluids Group, Max Planck Center Twente for Complex Fluid Dynamics, J. M. Burgers Center for Fluid Dynamics, University of Twente, P. O. Box 217, 7500 AE Enschede, The Netherlands</footnote>
for their good accuracy at long range and for their low computational cost (Gilbert and White, 1989). They consist in solving a one-way-wave equation in the frequency domain usually in a two-dimensional (2D) geometry, although the formulation holds in three dimensions (3D). PE methods have thus become the state of the art for wind turbine noise propagation (Barlas et al., 2017b; Kayser et al., 2020). Nevertheless, they suffer from several limitations. First, the solution is only valid for propagation angles close to the main propagation direction. Second, the one-way wave equation neglects back-scattering. This could lead to inaccuracies in the presence of topography with steep slopes. Finally, approximations are usually made in order to consider propagation in an inhomogeneous moving medium. The most common approach is to consider an atmosphere at rest with an effective sound speed, that includes the effect of both temperature and of the wind velocity gradient. This approach can be inaccurate if the wind speed is too high or if its vertical component is not negligible (Dallois et al., 2001). Improved PE methods consider the effect of the ABL on sound propagation by including the mean flow terms when deriving the one-way wave equation (Dallois et al., 2001; Ostashev et al., 2020). The range of application of PE methods for wind turbine noise has already been investigated in several studies, by comparing PE results to analytical solutions or measurements (Kayser et al., 2023; Lee et al., 2016; Nyborg et al., 2022). However, validation against analytical solutions is restricted to simple cases. In addition, outdoor experiments are not perfectly controlled, as, for instance, the wind field is only known partially or the wind turbine noise sources have to be modeled. Comparison against measurements remains thus overall limited.

A more advanced sound propagation method is to solve the linearized Euler equations (LEE) directly. This method is also widely used for outdoor propagation (Blumrich and Heimann, 2002; Dragna and Blanc-Benon, 2014; Salomons et al., 2002; Van Renterghem, 2014). It takes the effect of the mean flow on sound propagation accurately into account. In particular, it considers the vertical component of the wind speed, which is neglected in the effective sound speed approach. The LEE method overcomes the PE limitations, which include restricted angular validity and restrictions on backscattering. In addition, as a time-domain method, it provides a broadband solution and can be used to consider unsteady effects. The main drawback is the higher computational cost of the LEE compared to the PE method, which explains why it has not yet been considered for wind turbine noise propagation.

The objectives of this paper are to introduce a LEE model for wind turbine noise propagation and to assess the advantages of using a LEE model with respect to state-of-the-art PE methods. The LEE model can account for topography, ground impedance, and inhomogeneous mean flow. For comparison, two different PE implementations are considered: a vector PE able to handle strong wind variation but limited to flat terrain, and

a PE formulation able to consider topography and wind gradients through an effective sound speed approach. We compare the results of the LEE and PE approaches for two realistic cases: first a wind turbine on flat ground and second a wind turbine located on top of a steep hill. For both cases, the mean flow in which the sound propagates is obtained from previously computed large eddy simulations (LES) by Liu and Stevens (2020) and an extended moving source model based on Amiet's theory (Cotté, 2019) is used.

The paper is organized as follows. In Sec. II the complete methodology for the wind turbine noise prediction is described including the LEE and PE methods. Then Sec. III details the two cases studied. In Sec. IV the SPL and AM obtained with both methods are compared for the two cases. Finally, concluding remarks are given in Sec. V.

## II. METHOD

### A. General methodology

The general framework used to compute wind turbine noise at a receiver location is summarized in Fig. 1. It is based on the coupling of three different models: a LES code used to obtain the mean wind velocity in the ABL, a source model based on Amiet's theory, and a propagation model. The propagation models studied in this work are then detailed in Sec. II B.

<figure><img src="image_3.png" /><figcaption>FIG. 1. Diagram of the complete prediction methodology</figcaption></figure>

First, the LES code is used to compute a realistic ABL. The interaction between the flow and the wind turbine is modeled with an actuator disk method. This code has been extensively tested both for the computation of realistic atmospheric boundary layers (Gadde et al., 2021) and for the simulation of wind farm flow, with or without topography (Gadde and Stevens, 2021; Liu and Stevens, 2020). Although the simulations are unsteady, only the mean velocity fields are used in the
following. By doing so the turbulence scattering is not considered in this study, although it is known to have an impact on the propagation of wind turbine noise (Barlas et al., 2017b). The flow data are then fed into an extended source model and a propagation model based on either the LEE or on the PE.

The source model derives from Amiet’s strip theory: Each blade is divided into several segments, considered as uncorrelated sources. For each segment the turbulent inflow and trailing edge noise are computed using the model developed by Tian and Cotté (2016). The resulting SPL in free field, denoted SPL<sub>ff</sub>, depends on the wind turbine geometry (chord length and type of the blade segment) and on the incoming velocity profile and turbulent spectrum. Additionally, to determine the propagation effects, such as ground reflection or atmospheric refraction, the SPL relative to the free field, denoted ΔL, must be computed for each source-receiver pair.

Finally, source and propagation effects are considered by computing the SPL at a receiver produced by one segment i at a given angular position of the blades β such that:

$$SPL^i(\mathbf{x}, \omega, \beta) = SPL_{ff}^i(\mathbf{x}, \omega, \beta) + \Delta L^i(\mathbf{x}, \omega, \beta) - \alpha(\omega)R, \quad (1)$$

where x is the receiver coordinates, ω the angular frequency, α the atmospheric absorption coefficient, and R the distance from source to receiver. The contributions of each segment are combined to obtain the total SPL at the receiver:

$$SPL(\mathbf{x}, \omega, \beta) = 10 \log_{10} \left( \sum_{i=1}^{N_s} 10^{SPL^i(\mathbf{x}, \omega, \beta)/10} \right), \quad (2)$$

where N<sub>s</sub> is the number of segments used to discretize the three blades. The overall SPL (OASPL) are then computed by summing the SPL over ω for a given angular position of the blades. In addition, amplitude modulation can be quantified by measuring the difference between the minimum and maximum OASPL values during one rotation of the blades. This method can become very expensive as the number of propagation simulations is proportional to the number of receivers, angular positions, and blade segments.

To limit the computational cost, the method proposed by Cotté (2019) considers a number of fictive source heights distributed over a vertical line in the rotor plane, as depicted in Fig. 2.a. The propagation simulations are performed only for these fictive sources. Subsequently, the SPL<sup>i</sup> in Eq. (1) is computed using the value of ΔL corresponding to the closest fictive source, see Fig. 2b. Note that the number of propagation simulations is reduced by considering a small number of source positions and by positioning all fictive sources in the same propagation plane. Hence, the same ΔL results are utilized to compute SPL for all receivers in this plane. This method is based on the assumption that the dimension of the wind turbine rotor is relatively small in relation to the propagation distances. In order to use this method

it is also important to ensure that the number of fictive source heights is sufficient to obtain convergence of the results i.e. that the average SPL and AM during a rotation do not depend on the number of sources used. In Cotté (2019), it was shown that a minimum of 7 source heights is necessary to obtain this convergence in the downwind direction.

<figure><img src="image_4.png" /><figcaption>FIG. 2. (a) Calculation of the SPL at a receiver due to one blade segment. (b) Total SPL from all segments is determined using the closest fictive source to compute the ΔL.</figcaption></figure>

## B. Acoustic propagation simulation

In this section the models used to predict the sound propagation in the ABL are introduced. First, the LEE model developed for this study is presented and then the PE implementations are described.

### 1. LEE

a. Acoustic model. This propagation model is based on the finite-difference solution of the LEE in a curvilinear mesh (Dragna and Blanc-Benon, 2014). Thus, the acoustic field is solved separately from the mean flow field and does not affect it. The LEE are derived from the conservation of mass, momentum, and energy for atmospheric propagation by neglecting terms
of order $(|\mathbf{V}_0|/c_0)^2$ (Ostashev et al., 2005),

$$ \begin{aligned} & \frac{\partial p}{\partial t} + (\mathbf{V}_0 \cdot \nabla)p + \rho_0 c_0^2 \nabla \cdot \mathbf{v} = 0, \\ & \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{V}_0 \cdot \nabla)\mathbf{v} + (\mathbf{v} \cdot \nabla)\mathbf{V}_0 + \frac{\nabla p}{\rho_0} = 0, \end{aligned} \quad (3) $$

where $t$ is the time; $p$ and $\mathbf{v} = (u, w)$ are the acoustic pressure and velocity; and $\rho_0$, $c_0$, and $\mathbf{V}_0 = (u_0, w_0)$ are the mean density, sound speed, and velocity. The 2D LEE in the $(x, z)$ plane can be written in conservative form as

$$ \mathbf{U}_t + \mathbf{E}_x + \mathbf{F}_z + \mathbf{H} = 0, \quad (4) $$

where $\mathbf{U} = [p, \rho_0 u, \rho_0 w]^T$ is the vector of unknowns; the partial derivatives are denoted $i_j = \partial i / \partial j$; and $\mathbf{E}$, $\mathbf{F}$, and $\mathbf{H}$ are the Eulerian fluxes defined as

$$ \begin{aligned} & \mathbf{E} = [u_0 p + \rho_0 c_0^2 u, u_0 \rho_0 u + p, u_0 \rho_0 w]^T, \\ & \mathbf{F} = [w_0 p + \rho_0 c_0^2 w, w_0 \rho_0 u, w_0 \rho_0 w + p]^T, \\ & \mathbf{H} = [-p(\nabla \cdot \mathbf{V}_0), \rho_0(\mathbf{v} \cdot \nabla)u_0, \rho_0(\mathbf{v} \cdot \nabla)v_0]^T. \end{aligned} \quad (5) $$

The acoustic source is a pulse with a Gaussian spatial distribution, introduced via the initial conditions of Eq. (4):

$$ \mathbf{U}(x, z, t = 0) = \left[ S_0 \exp \left( -\frac{R^2}{B^2} \right), 0, 0 \right]^T, \quad (6) $$

where $R = \sqrt{x^2 + (z - z_S)^2}$, $B^2 = (3\Delta x)^2 / \log(2)$, $z_S$ is the height of the center of the pulse, $S_0 = 1$ Pa is the amplitude of the source, and $\Delta x$ is the streamwise grid spacing. This source aims to represent a broadband monopole, with frequency content up to $f = 0.6c_0/B$; see appendix Eq. (A.22).

**b. Curvilinear transformation.** To take into account a ground profile $h$, a transformation of the coordinate system is applied from Cartesian coordinates $(x, z)$ to curvilinear coordinates $(\xi, \eta)$ such that

$$ \begin{aligned} x(\xi, \eta) &= \xi, \\ z(\xi, \eta) &= h(\xi) + \frac{\eta}{z_{\max}} [z_{\max} - h(\xi)], \end{aligned} \quad (7) $$

where $h(x)$ is the terrain elevation and $z_{\max}$ is the maximum height of the domain. Hence, Eq. (4) becomes:

$$ \left(\frac{\mathbf{U}}{J}\right)_t + \left(\frac{\xi_x \mathbf{E} + \xi_z \mathbf{F}}{J}\right)_\xi + \left(\frac{\eta_x \mathbf{E} + \eta_z \mathbf{F}}{J}\right)_\eta + \frac{\mathbf{H}}{J} = 0, \quad (8) $$

where $J = |\xi_x \eta_z - \xi_z \eta_x|$ is the Jacobian of the transformation. This transformation was proposed by Gal-Chen and Somerville (1975) to recover a flat top boundary.

**c. Numerical scheme.** The formulation in Eq. (8) can be written as

$$ \mathbf{U}_t = \mathcal{F}(\mathbf{U}), \quad (9) $$

where $\mathcal{F}$ is a function of $\mathbf{U}$ and its spatial derivatives. A fourth-order six-step Runge-Kutta (RK) algorithm (Berland et al., 2006) is used to integrate the solution

<figure><img src="image_5.png" /><figcaption>FIG. 3. Snapshots of the moving frame at three different instants in time as the wave propagates inside the domain.</figcaption></figure>

from $\mathbf{U}(t_n)$ to $\mathbf{U}(t_{n+1})$ with $t_n$ the discrete time. The spatial derivatives needed at each step of the RK algorithm are computed using a fourth order 11-point stencil finite-difference centered scheme (Bogey and Bailly, 2004) whose coefficients are optimized to minimize the dispersion error over a large range of wavenumber. For the points near the boundary, non-centered 11-points stencil schemes are used (Berland et al., 2007). Selective filters are also applied to remove grid-to-grid oscillations.

**d. Moving frame.** A moving frame method (Dragna and Blanc-Benon, 2014) is used to keep an affordable computational cost for long range propagation. In this method, the distance between the wavefront and the right boundary of the domain is computed at each iteration and the acoustic variables are shifted to maintain this distance constant throughout the simulation (see Fig. 3). This approach allows for a large reduction of the computational domain size. Thus, in the simulations, a 300 m × 300 m moving frame is used to calculate the sound propagation in a 3000 m × 300 m domain, which reduces the computational cost by one order of magnitude.

**e. Boundary conditions.** A broadband impedance condition is used at the bottom boundary of the domain to model realistic ground absorption. The impedance condition is developed in the time domain (Rienstra, 2006) such that it writes as a convolution. By taking the surface impedance as a rational function of the frequency, it is possible to substitute this convolution by additional differential equations that can be solved along the RK scheme (Troian et al., 2017). This saves computational resources, while preserving high-order accuracy. In this work, the method was not implemented on the impedance but on the reflection coefficient for stability reasons.

At the top of the domain a perfectly matched layer (PML) is implemented to simulate unbounded atmosphere. Inside this layer, the derivative in $z$ in Eq. (4) is modified such that

$$ \frac{\partial}{\partial z} \rightarrow \frac{1}{\kappa + \sigma/i\omega} \frac{\partial}{\partial z}, \quad (10) $$
with

$$ \begin{aligned} \sigma(z) &= \sigma_0 [(z - z_{\text{PML}})/L]^3, \\ \kappa(z) &= 1 + (\kappa_0 - 1) [(z - z_{\text{PML}})/L]^3, \end{aligned} \quad (11) $$

where $z_{\text{PML}}$ is the coordinate at the bottom of the layer and $L$ is its thickness. As for the impedance condition, this transformation can be written in the time domain as a convolution that is integrated along the RK scheme with additional differential equations (Cosnefroy, 2019; Komatitsch and Martin, 2007).

Finally, on the left side of the domain, the waves must also leave the computational domain without producing any reflections. To achieve this, the sound speed is gradually reduced in a thin layer to ensure that potentially reflected waves move slower than the domain, which is moving at the speed of sound. Hence, they can not re-enter the domain (Cosnefroy, 2019).

f. Post-processing.. During the computation, pressure time signals $p(t)$ are recorded at several receiver locations. Then, the Fourier transform $\hat{p}(\omega)$ of each signal is computed, and the relative sound pressure is retrieved by dividing $\hat{p}(\omega)$ by the solution for a Gaussian pulse in free field,

$$ \Delta L(\mathbf{x}, \omega) = 10 \log_{10} \left( \frac{|\hat{p}(\mathbf{x}, \omega)|^2}{|\hat{p}_{\text{ff}}(\mathbf{x}, \omega)|^2} \right), \quad (12) $$

where $\hat{p}_{\text{ff}}$ is given by:

$$ \hat{p}_{\text{ff}}(\mathbf{x}, \omega) = \frac{\pi k_0 B^2 S_0}{4c_0} \exp\left(-\frac{k_{\text{eff}}^2 B^2}{4}\right) H_0^{(1)}(k_0 R), \quad (13) $$

$H_0^{(1)}$ is the zeroth order Hankel function of the first kind, $k_0$ is the wave number, and $k_{\text{eff}} = k_0/(1+M)$ with $M = u_0/c_0$ the Mach number at the wind turbine hub's height. The term $k_{\text{eff}}$ accounts for the effect of the flow on the source which induces a shift in frequency. The derivation of this correction term can be found in the Appendix.

Once $\Delta L$ is computed, it can be used in the SPL prediction methodology explained in Sec. II A.

## 2. PE method

a. PE formulations.. Two 2D PE formulations are implemented. The first one is a curvilinear formulation of the wide-angle PE with the effective sound speed approximation (WAPE$_{topo}$). It is derived from the methodology described by Sack and West (1995) for the generalized terrain parabolic equation (GT-PE), but with the coordinate transformation shown in Eq. (7). The WAPE$_{topo}$ allows us to consider topography with slope that does not exceed around 30° (Salomons, 2001). The moving atmosphere is modeled by defining an effective sound speed profile such that

$$ c_{\text{eff}}(x, z) = c_0(T_0(x, z)) + u_0(x, z), \quad (14) $$

where $c_0$ is the sound speed that depends on the mean temperature $T_0$. It is worth mentioning that the WAPE$_{topo}$ is equivalent to the classical effective sound

speed WAPE formulation for flat ground. The accuracy of the WAPE$_{topo}$ strongly depends on quantities computed at the ground (ground profile and its derivatives, ground surface impedance...), where small numerical errors can appear and accumulate (Sack and West, 1995). Errors can also arise from the effective sound speed approach with phase errors accumulating over the distance.

The second PE is a wide-angle vector PE formulation (WAPE$_{vec}$), that incorporates a vector wind field and is valid for arbitrary high Mach numbers. It is presented in Sec. VI. A of Ostashev et al. (2020). Here the velocity field is considered with fewer assumptions. However, it was not derived in curvilinear coordinates and hence can not consider topography. While WAPE$_{vec}$ can consider both the vertical and horizontal components, only the horizontal component of the wind velocity $u_0$ is considered in our implementation.

The two methods employ second-order finite-difference schemes in the $z$-direction and a Crank-Nicolson algorithm to advance the solution to $x + \Delta x$ from the solution at $x$. Numerically, this involves inverting a tridiagonal matrix at each step of the resolution, which is done efficiently using the Thomas algorithm. The starting field is the second-order starter presented in Salomons (2001), which represents a monopole source.

b. Boundary condition. At the top of the domain a PML is implemented (Collino, 1997) by modifying the partial derivative along the vertical coordinate as for the LEE method [see Eq. (10)]. At the bottom of the domain an impedance boundary condition is used. The implementation of the boundary conditions leads to modifications of the matrix coefficients, while conserving the tridiagonal shape. The boundary conditions are similar for the two PE formulations but, because of the differences in the equations, the resulting modifications of the matrix coefficients are not identical.

## III. CASES STUDIED

The study examines two scenarios: one with a wind turbine on a flat surface and one with a wind turbine positioned on a hilltop. The latter scenario is particularly relevant, as wind turbines are often installed on hills to maximize energy output. Therefore, it is crucial to accurately model sound propagation in such complex settings. For both cases, LES results from Liu and Stevens (2020) for a truly neutral ABL are used. The mean flow fields from these simulations are normalized by the friction velocity $u_*$ and by the diameter of the wind turbine. In this work, the wind turbine hub height and diameter are both set to 100 m. The friction velocity is set to $u_* = 0.512 \text{ m s}^{-1}$ to get a wind velocity equal to 10 m s$^{-1}$ at the hub height. For simplicity, the temperature and the sound speed are assumed constant, which is a valid assumption when considering a neutral atmosphere but would not be for a stable or unstable one. For the second case, the hill is defined such that:

$$ h(x) = h_{\text{max}} \cos^2 \left( \frac{\pi x}{2l} \right) \quad \text{for } -l < x < l, \quad (15) $$
<figure><img src="image_6.png" /><figcaption>FIG. 4. Mean axial velocity field $u_0$ for the two cases of interest: (a) a wind turbine over flat ground, (b) on top of a hill. Mean vertical velocity field $w_0$ for (c) a wind turbine over flat ground, (d) on top of a hill.</figcaption></figure>

where $h(x)$ is the terrain elevation, $h_{max} = 100$ m, and $l = 260$ m is the half-width of the hill. The mean axial velocity component is plotted for both cases in Fig. 4(a) and 4(b). A velocity deficit can be observed just after the turbine in both cases. In the first case, a ducting effect is expected due to the wind turbine wake (Barlas et al., 2017b). The acoustic waves are trapped inside the wake and then redirected towards the ground when the unperturbed ABL is recovered. This generates a focusing zone at the ground with high SPL. Ducting in the presence of the hill is expected to be stronger as the wake is more pronounced and more directed towards the ground. The vertical component of the wind speed, plotted Fig. 4c, is small in the flat case (less than 4 % of $u_0$), so that the vector PE and the LEE are expected to produce similar results. In the case with topography, the hill induces a vertical component of the wind speed just before the wind turbine of a few meters per second, shown in Fig. 4d, which is not taken into account by the WAPE$_{topo}$. This difference of a few meters per second upstream of the hill could introduce some discrepancies between the LEE and WAPE$_{topo}$ methods in the final OASPL. In addition, errors could arise for both cases due to the one-way approximation and the angular validity of the two PE methods.

The variable porosity model (Attenborough et al., 2011) is used with a flow resistivity of 50 kN s m⁻⁴ and a porosity change rate of 100 m⁻¹ to model a grassy ground. The PE simulations are performed for a set of frequencies used to compute third octave band spectrum between 50 Hz and 1 kHz. Values of ΔL for the same frequencies are extracted from the broadband results of

the LEE simulation for comparison. These frequencies are gathered in Table I. For all simulations, the numer-

<table><caption>TABLE I. Frequencies *f* used to compute the third octave band spectrum with the PE methods. *f*<sub>c</sub> is the central frequency of each band.</caption><tr><td>f<sub>c</sub> (Hz)</td><td>50</td><td>63</td><td>80</td><td>100</td><td>125</td><td>160</td><td>200</td><td>250</td><td>315</td><td>400</td><td>500</td><td>630</td><td>800</td><td>1000</td></tr><tr><td>f (Hz)</td><td>50</td><td>63</td><td>80</td><td>100</td><td>125</td><td>160</td><td>192</td><td>241</td><td>297</td><td>373</td><td>467</td><td>588</td><td>741</td><td>926</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>208</td><td>260</td><td>315</td><td>391</td><td>489</td><td>616</td><td>770</td><td>962</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>334</td><td>409</td><td>512</td><td>645</td><td>800</td><td>1000</td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>429</td><td>536</td><td>675</td><td>831</td><td>1039</td><td></td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>864</td><td>1080</td><td></td></tr></table>

ical parameters are set to obtain accurate results up to 1 kHz and 3 km downwind of the turbine. For the LEE, this implies setting a grid size of 0.05 m and a CFL of 0.5, which corresponds to 127000 time iterations to complete the simulation with a computational moving domain of 36 × 10⁶ points. The computational time for one simulation is approximately 1200 central processing units (CPU) hours, which is equivalent to 3 days on a 16-core machine. For WAPE$_{vec}$ and WAPE$_{topo}$, the grid size depends on the computed frequency. For the flat case convergence was reached for Δx = Δz = λ/10, which requires 1 CPU hour to compute the results at all frequencies. The hill case requires a higher resolution to capture sharp gradients and convergence is reached for Δx = Δz = λ/50, which leads to an increased computational time of 30 CPU hours.

The source model used for the wind turbine is the same as described in Tian and Cotté (2016). The only difference is that the wind turbine is scaled up to be of 100 m in diameter instead of 93 m. The blade is decomposed into eight segments, and 36 angular positions are considered for one rotation of the blades. Finally, it was found by Cotté (2019) that seven fictive source heights are sufficient to obtain a convergence on SPL and AM for the flat case at ground level using the approach described in Sec. II A. However, for the case with a hill, 30 source heights are required for convergence. This significant increase in the number of sources needed to achieve convergence can be explained by the strong dependence between the focusing pattern and the source height in the presence of the hill. Thus, changing the source height by a few meters can affect the position of the focusing zone at the ground of more than 100 m. Therefore, it is necessary to reduce the step between the simulated source heights.

# IV. COMPARISON

## A. Over flat ground

### 1. Relative SPL

The value of ΔL represents the effect of the ground absorption and the mean flow on the propagation with respect to the solution in the free field. The ΔL fields are
<figure><img src="image_7.png" /><figcaption>FIG. 5. $\Delta L$ computed with LEE (left) and WAPE$_{vec}$ (right) at 100 Hz for several source heights (circles): 142 m (top), 100 m (middle), and 58 m (bottom).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 6. $\Delta L$ computed for a source at 142 m and a receiver at 2 m for three different frequencies.</figcaption></figure>

depicted in Fig. 5 for LEE and WAPE$_{vec}$ at 100 Hz for three source heights corresponding to the highest (142 m) and lowest (58 m) fictive sources and to the fictive source at hub height (100 m). The results for the WAPE$_{topo}$ are not plotted as they are almost identical with those of the WAPE$_{vec}$. In this case, the two main effects are the occurrence of constructive and destructive interference patterns due to ground reflection and the ducting of acoustic waves by the wake of the wind turbine. In particular, note the large increase in $\Delta L$ at $x = 1.7$ km at the ground for the source at 142 m height. It is clear that the same interference patterns and the same ducting are obtained in the WAPE$_{vec}$ and LEE results. The influence of the wind speed gradient and of the wake length and intensity on the propagation for different source heights is well captured with both methods as the $\Delta L$ levels are similar. Still, a small difference between the two methods can be observed at the very beginning of the domain ($x < 100$ m). For the LEE method, the moving frame does not allow the reflected wave to reach the top of the domain close to the source, leading to a zone without any interferences and a value of $\Delta L$ close to 0. In the case

of PE, the classic cone due to the angle of validity of the method can be seen at the very beginning of the domain. Hence, the two methods do not produce the same results in the near field.

The relative sound pressure level is plotted for a line of receivers 2 m above the ground in Fig. 6 for three different frequencies. At all frequencies, a strong peak is visible at 1.7 km which corresponds to the distance where the focused wave hits the ground. The values of $\Delta L$ obtained with WAPE$_{topo}$ and WAPE$_{vec}$ methods are compared with the LEE results. The curves are almost identical for the three methods. The difference between the $\Delta L$ obtained with the three methods does not exceed 0.5 dB for $x > 1$ km. At 512 Hz and 1000 Hz a small shift is found between the position of the peaks in the WAPE$_{topo}$ and WAPE$_{vec}$ results. The latter seems, as expected, to better account for the wind field and the results are closer to those of the LEE. For other source or receiver heights (not plotted here), it was observed that the difference between WAPE$_{topo}$ and LEE results tends to increase with distance. It was also seen that as the
<figure><img src="image_9.png" /><figcaption>FIG. 7. $\Delta L$ difference computed between the LEE and (a,b) WAPE$_{vec}$ and (c,d) WAPE$_{topo}$ for (a,c) an atmosphere at rest and (b,d) for an ABL profile with 10 m s$^{-1}$ at hub height. The results are averaged over $N_f = 35$ frequencies, see Table I, and are computed for flat terrain.</figcaption></figure>

interference pattern gets more complex, the difference is more visible.

The average difference of $\Delta L$ is plotted in Fig. 7 to quantify the error made with PE formulations. This difference is defined with:

$$ \epsilon = \frac{1}{N_f} \sum_{i=1}^{N_f} |\Delta L_{\text{LEE}}(f_i) - \Delta L_{\text{PE}}(f_i)|, \quad (16) $$

with $N_f$ the total number of frequencies computed. In order to assess the effect of the mean flow on the results, an additional simulation for an atmosphere at rest is performed for all three methods. The difference for the atmosphere at rest is shown in Fig. 7a and 7c. Note that the error is almost identical for the WAPE$_{vec}$ and the WAPE$_{topo}$. The error is very large close to source due to the angle of validity of the PE and can even go above 3 dB. Further away, the propagation angle lies within the angle of validity of the WAPE$_{topo}$, and the error reduces to less than 0.5 dB. The introduction of a mean flow does not significantly modify the error for the WAPE$_{vec}$ (Fig. 7b). The effect of the flow is very well accounted for in this method, as expected. On the contrary, for WAPE$_{topo}$ the addition of a mean flow (Fig. 7d) greatly increases the error. The difference is not only significant in the near field but in the entire domain with differences going above 3 dB. It can be observed that close to the ground the error is still relatively small, which corresponds to the case presented in Fig. 6. It is worth noting that, because we are considering a flat ground, the error in the WAPE$_{topo}$ method is only due to the effective sound speed approach.

<figure><img src="image_10.png" /><figcaption>FIG. 8. OASPL at $\beta = 0^\circ$ for the flat case computed using (a) LEE and (b) WAPE$_{topo}$.</figcaption></figure>

## 2. SPL

In this section, the SPL obtained by combining $\Delta L$ with the source model described in Sec. II A are compared for the LEE and the two PE methods. The OASPL field for one angular position of the blades obtained from LEE and WAPE$_{topo}$ methods is presented in Fig. 8. It corresponds to the contribution of all blade segments for one angular position (here $\beta = 0^\circ$, i.e. for one of the three blades pointing upwards) summed over the frequency bands between 50 Hz and 1 kHz. The results from WAPE$_{vec}$ are omitted as they are again almost identical with those of WAPE$_{topo}$. Several zones of large OASPL produced by the different sources distributed along the blades can be observed in both cases. The OASPL obtained from the WAPE$_{topo}$ and the LEE methods for several positions of the blades is available in video Mm. 1. As the position of the blades changes, different $\Delta L$ calculations are activated, leading to distinct focusing zones on the ground. This mechanism is responsible for the amplitude modulation in the far field (Barlas et al., 2018). Hence, it is not expected that strong discrepancies appear as the error observed in Fig. 7 would tend to average out with an extended broadband source.

Mm. 1. Evolution of the OASPL with the angle $\beta$, in the flat case for the WAPE$_{topo}$ and the LEE methods.

The mean OASPL over one rotation for receivers located 2 m above the ground is shown in Fig. 9a. The minimum and maximum values reached during the rotation are also delimited by color patches. The mean OASPL obtained with the three methods is very similar. The peak observed 1.7 km from the source is well captured by both PE methods. The small phase shift in $\Delta L$ observed for the WAPE$_{topo}$ is still present in the OASPL prediction but is less pronounced. The two amplitude modulation zones from 1.25 km to 1.9 km and from 1.9 km to 3 km are similar with all methods. It corresponds to the areas where the SPL varies due to the motion of the blades.

The narrow band spectrum at three downstream locations is shown in Fig. 9b for the LEE and the two PE methods.
<figure><img src="image_11.png" /><figcaption>FIG. 9. (a) Averaged OASPL over one rotation computed with the three propagation methods, minimum and maximum delimited with color patches. (b) Averaged spectrum obtained over a full rotation of the turbine blades for receivers at a height of 2 m, positioned at distances of *x* = 500 m, *x* = 1.5 km, and *x* = 2.9 km downstream of the turbine (listed from top to bottom).</figcaption></figure>

formulations. Spectra are nearly identical over the considered frequency range. The dips induced by ground absorption are equally captured. The error does not increase with frequency or distance, which shows the good convergence of the three methods for this frequency range and domain length.

## B. On top of a hill

### 1. Relative SPL

The presence of the hill influences sound propagation by its geometry and by inducing a more complex mean flow, which changes the refraction of acoustic waves in the ABL. This refraction pattern is shown in Fig. 10 for the LEE and the WAPE<sub>topo</sub> methods and for three source heights. As for the flat case, different focusing zones can be identified depending on the source height. The presence of a caustic at the bottom of the hill can also be observed. The interference dips look more pronounced for the WAPE<sub>topo</sub> and the levels are slightly higher for the LEE. Furthermore, important differences can be seen in the results, especially at 100 Hz, for a line of receivers at 2 m height and for a source at 242 m (Fig. 11). These

large discrepancies are due to a combined effect of the hill and of the mean flow. In fact, the steep slope of the hill is the main cause of the error in this case as discrepancies were already visible without any mean flow (not shown here). The maximum slope of the hill (reached at *x* = 130 m) is equal to 31°, which is at the limit of the WAPE<sub>topo</sub> validity range (Salomons, 2001). It is worth noting that the regularity of the ground profile or the precision with which its derivatives are calculated are not an issue here as they are analytical. The mean error in the whole domain for all frequencies is between 3 dB and 4 dB for all heights, which is 2 dB higher than for the flat case. Still it can be observed that the effect of the hill and the mean flow is the same for the LEE and WAPE<sub>topo</sub>. Thus, the shape of the peak at 1.8 km is similar for both methods at 512 Hz and 1000 Hz. Even if the WAPE<sub>topo</sub> introduces errors in this case, the effects of the ABL and the hill remain well captured both in terms of amplitude and position of the focusing.

### 2. SPL

As for the flat terrain case, the OASPL map depicted in Fig. 12 corresponds to the superposition of SPL from the different sources distributed along the blades. The caustic observed at the bottom of the hill in Fig. 10 is still present, as well as several sound focusing zones. A notable difference from the flat terrain case is that the focusing is significantly stronger. Additionally, as discussed in Sec. III, the presence of the hill results in greater amplitude modulation due to blade movement. This can be observed in Mm. 2, and it is clear that both methods capture this phenomenon.

Mm. 2. Evolution in the hill case of the OASPL with the angle $\beta$ for the WAPE<sub>topo</sub> and the LEE method.

To further investigate the discrepancies between the LEE and WAPE<sub>topo</sub> methods, the mean OASPL over one rotation is plotted in Fig. 13a along with the minimum and maximum OASPL for a line of receivers at 2 m height. The burst in mean OASPL at 1.8 km is similar for the two methods as well as the maximum of amplitude modulation with a strong peak at 2.5 km. Hence, it is clear that the same propagation effects are captured with both methods. However, the OASPL computed with LEE are slightly higher than those computed with the WAPE<sub>topo</sub> (less than 1 dB). Finally SPL spectra at three receiver positions are presented in Fig. 13b. The figure demonstrates that WAPE<sub>topo</sub> and LEE exhibit excellent agreement over the entire frequency range at *x* = 1.5 km and *x* = 2.5 km downstream of the turbine. However, for the receiver situated close to the turbine (at 500 meters downstream), WAPE<sub>topo</sub> slightly underpredicts the SPL compared to LEE.
<figure><img src="image_12.png" /><figcaption>FIG. 10. $\Delta L$ computed with LEE (left) and WAPE$_{topo}$ (right) at 100 Hz for several source heights (circles): 242 m (top), 200 m (middle), and 158 m (bottom).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 11. $\Delta L$ computed for a source at 242 m and a receiver at 2 m for three different frequencies.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 12. OASPL at $\beta = 0^{\circ}$ for the hill case computed using (a) LEE and (b) WAPE$_{topo}$.</figcaption></figure>

V. CONCLUSION

The use of a new method based on the linearized Euler equations (LEE) for predicting wind turbine noise propagation was investigated. The method includes a flow model based on LES and an extended source model

based on Amiet's theory. Comparison with the state-of-the-art PE methods (WAPE$_{vec}$ and WAPE$_{topo}$) was performed for two cases: a baseline case with a wind turbine over flat ground and a more complex case where the wind turbine is positioned at the top of a hill. For the baseline case, the sound pressure level relative to the free field and the OASPL are almost identical when using LEE and the WAPE$_{vec}$. It was noticed that, as expected, the WAPE$_{topo}$ introduces phase errors due to the effective sound speed approach. However, this has a minor impact on the prediction of the sound pressure levels at the ground and of OASPL. It is worth noting that the WAPE$_{vec}$ is still preferable as it is derived with less assumptions and that it could be improved further by considering the vertical component of the wind speed. In the case with topography, the WAPE$_{topo}$ generates a noticeable error in the relative sound pressure level computed for a given frequency (of the order of 3 dB). This is due to the steep slope of the hill. Nevertheless, the OASPL and AM obtained with the WAPE$_{topo}$ remain very close to those obtained with the LEE method. The main effects of the flow on the propagation are still well simulated
<figure><img src="image_15.png" /><figcaption>FIG. 13. (a) Averaged OASPL over one rotation computed with the LEE and WAPE<sub>topo</sub> methods, minimum and maximum delimited with color patches. (b) Averaged spectrum obtained over a full rotation of the turbine blades for receivers at a height of 2 m, positioned at distances of *x* = 500 m, *x* = 1.5 km, and *x* = 2.5 km downstream of the turbine (listed from top to bottom).</figcaption></figure>

by the WAPE<sub>topo</sub>, and that the shifts in the interference pattern observed at each frequency tend to average out for the prediction of overall levels.

Hence, we find that for situations with topography PE methods in general and WAPE<sub>topo</sub> in particular provide a suitable first approach to determine sound propagation from wind turbines.. It is still worth noting the advantages of using LEE. First, a time-domain solution is obtained which allows one to compute a broadband SPL spectrum. The flow is taken into account with fewer assumptions and higher wind speed can be considered without introducing errors due to the effective sound speed approach. The main drawback of this method is its computation cost, especially when numerous source heights must be considered.

## ACKNOWLEDGMENTS

The authors thank Luoqin Liu for providing access to the LES data of Liu and Stevens (2020). This work was performed within the framework of the LABEX CeLyA (ANR-10-LABX-0060) of Université de Lyon, within the program “Investissements d’Avenir” (ANR-16-IDEX-0005) operated by the French National Research Agency

(ANR). The authors were granted access to the HPC resources of PMCS2I (Pôle de Modélisation et de Calcul en Sciences de l'Ingénieur et de l'Information) of Ecole Centrale de Lyon, PSMN (Pôle Scientifique de Modélisation Numérique) of ENS de Lyon and P2CHPD (Pôle de Calcul Hautes Performances Dédiés) of Université Lyon I, members of FLMSN (Fédération Lyonnaise de Modélisation et Sciences Numériques), partner of EQUIPEX EQUIP@MESO. This work was supported by the Franco-Dutch Hubert Curien partnership (Van Gogh Programme No. 49310UM). For the purpose of Open Access, a CC-BY public copyright license has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission. The authors declare no conflicts of interest. Data from numerical simulations are available from the authors on reasonable request.

## APPENDIX: CALCULATION OF THE NOISE LEVEL RELATIVE TO THE FREE FIELD FOR THE LEE

The appendix details the derivation of the formula in Eq. (12) for calculating the noise levels relative to the free-field from the LEE solution.

### 1. Analytical derivation

We first derive the analytical solution in the frequency domain and in far-field for an impulsive source in a 2-D homogeneous uniformly moving medium and in free-field, as sketched in Fig. 14. The derivation closely follows that presented in Ostashev et al. (2005) for the case of a monochromatic point source.

<figure><img src="image_16.png" /><figcaption>FIG. 14. Sound propagation from a spatially distributed source in a 2-D homogeneous uniformly moving medium.</figcaption></figure>

### Solution as a convolution

A Cartesian system of coordinates **x** = (x, z) is used. A moving homogeneous atmosphere with a constant mean flow **V**_0 = (u_0, 0) is considered. The initial conditions are p(**x**, t = 0) = S(**x**) and **v**(**x**, t = 0) = 0 and the source spatial distribution S(**x**) is centered at the origin. Incorporating the initial conditions as source
terms, the LEE in Eq. (3) are written as:

$$ \frac{\partial p}{\partial t} + u_0 \frac{\partial p}{\partial x} + \rho_0 c_0^2 \nabla \cdot \mathbf{v} = S(\mathbf{x})\delta(t), \quad (A.1) $$

$$ \frac{\partial \mathbf{v}}{\partial t} + u_0 \frac{\partial \mathbf{v}}{\partial x} + \frac{1}{\rho_0} \nabla p = 0, \quad (A.2) $$

with $\delta$ the Dirac delta function. Combining Eqs. (A.1) and (A.2) leads to:

$$ \left( \frac{\partial}{\partial t} + u_0 \frac{\partial}{\partial x} \right)^2 p - c_0^2 \Delta p = \left( \frac{\partial}{\partial t} + u_0 \frac{\partial}{\partial x} \right) S(\mathbf{x}) \delta(t). \quad (A.3) $$

To translate the problem into the frequency domain, we employ the Fourier transform:

$$ \hat{p}(\mathbf{x}, \omega) = \int_{-\infty}^{\infty} p(\mathbf{x}, t) e^{i\omega t} dt. \quad (A.4) $$

Taking the Fourier transform of Eq. (A.3) and dividing by $c_0^2$ gives:

$$ \Delta \hat{p} - \left(-i k_0 + M \frac{\partial}{\partial x}\right)^2 \hat{p} = -\frac{1}{c_0} \left(-i k_0 + M \frac{\partial}{\partial x}\right) S(\mathbf{x}) \quad (A.5) $$

with $k_0 = \omega/c_0$ and $M = u_0/c_0$. Substituting in Eq. (A.5) the function $\hat{\phi}$ defined by:

$$ \hat{p} = -\frac{1}{c_0} \left(-i k_0 + M \frac{\partial}{\partial x}\right) \hat{\phi} \quad (A.6) $$

leads to:

$$ \Delta \hat{\phi} - \left(-i k_0 + M \frac{\partial}{\partial x}\right)^2 \hat{\phi} = S(\mathbf{x}). \quad (A.7) $$

For simplifying the previous equation, we use the transformation of coordinates $\mathbf{X} = (X, Z)$, with $x = X/\gamma$, $z = Z$ and $\gamma = 1/\sqrt{1-M^2}$. Finally, introducing the function $\hat{\psi}$ defined by $\hat{\phi} = \exp(-iK_0MX)\hat{\psi}$ with $k_0 = K_0/\gamma$ in Eq. (A.7) yields:

$$ \left( \frac{\partial^2}{\partial X^2} + \frac{\partial^2}{\partial Z^2} + K_0^2 \right) \hat{\psi} = Q(\mathbf{X}). \quad (A.8) $$

with $Q(\mathbf{X}) = \exp(iK_0MX)S(\mathbf{X})$. Eq. (A.8) corresponds to the 2D inhomogeneous Helmholtz equation. Its solution is written as a convolution of the source term $Q(\mathbf{X})$ and the Green's function, which gives:

$$ \hat{\psi}(\mathbf{X}) = -\frac{i}{4} \int H_0^{(1)}(K_0|\mathbf{X} - \mathbf{X}'|) Q(\mathbf{X}') d\mathbf{X}', \quad (A.9) $$

with $H_0^{(1)}$ the zeroth order Hankel function of the first kind.

## Far-field approximation

In order to evaluate the convolution in Eq. (A.9), a far-field approximation, known as the Fraunhofer approximation, is performed. It is assumed that the source-receiver distance is large compared to the characteristic size of the source B, i.e. $|\mathbf{X}| \gg B$. As a consequence, the source-receiver distance is approximated at

the zeroth order in the amplitude of the integrand, i.e. $|\mathbf{X} - \mathbf{X}'| = |\mathbf{X}|$, and at the first order in its phase, i.e. $|\mathbf{X} - \mathbf{X}'| = |\mathbf{X}| - (\mathbf{X} \cdot \mathbf{X}')/|\mathbf{X}|$. Under the Fraunhofer approximation, one has:

$$ H_0^{(1)}(K_0|\mathbf{X} - \mathbf{X}'|) = H_0^{(1)}(K_0|\mathbf{X}|) \exp(-iK_0|\mathbf{X}|) \exp\left[iK_0\left(|\mathbf{X}| - \frac{\mathbf{X} \cdot \mathbf{X}'}{|\mathbf{X}|}\right)\right], \quad (A.10) $$

which allows us to express the function $\hat{\psi}$ in Eq. (A.9) as:

$$ \hat{\psi}(\mathbf{X}) = \frac{-i}{4} H_0^{(1)}(K_0|\mathbf{X}|) \int S(X', Z') \exp[-iK_0(\cos \Theta - M)X'] \exp(-iK_0 \sin \Theta Z') dX' dZ' \quad (A.11) $$

Returning in the physical space, we get:

$$ \hat{\psi}(\mathbf{x}) = \frac{-i\gamma}{4} H_0^{(1)}(k_0\gamma\sqrt{\gamma^2 x^2 + z^2}) \int S(x', z') \exp[-i k_0 \gamma^2 (\cos \Theta - M) x'] \exp(-i k_0 \gamma \sin \Theta z') dx' dz', \quad (A.12) $$

where $\cos \Theta = X/|\mathbf{X}|$ and $\sin \Theta = Z/|\mathbf{X}|$ are related to $\cos \theta = x/|\mathbf{x}|$ and $\sin \theta = z/|\mathbf{x}|$ by:

$$ \cos \Theta = \frac{\cos \theta}{\sqrt{1 - M^2 \sin^2 \theta}}, \quad \sin \Theta = \frac{\sin \theta}{\gamma \sqrt{1 - M^2 \sin^2 \theta}}. \quad (A.13) $$

Using Eq. (A.12), we can write an expression for the function $\hat{\phi}$ in compact form as:

$$ \hat{\phi}(\mathbf{x}) = \frac{-i\gamma}{4} e^{ik_0 \gamma^2 M x} H_0^{(1)}(k_0 \gamma \sqrt{\gamma^2 x^2 + z^2}) \tilde{S}(k_0 \gamma^2 (\cos \Theta - M), k_0 \gamma \sin \Theta) \quad (A.14) $$

with $\tilde{S}$ the spatial Fourier transform of $S(\mathbf{x})$:

$$ \tilde{S}(k_x, k_z) = \int S(x, z) e^{-ik_x x - i k_z z} dx dz. \quad (A.15) $$

Finally, the pressure in the frequency domain is obtained from Eqs. (A.6) and (A.14):

$$ \hat{p}(\mathbf{x}) = A(\omega, M, \theta) \frac{-i\gamma^3}{4} e^{ik_0 \gamma^2 M x} \left( H_0^{(1)}(k_0 \gamma \sqrt{\gamma^2 x^2 + z^2}) - \frac{i\gamma x M}{\sqrt{\gamma^2 x^2 + z^2}} H_1^{(1)}(k_0 \gamma \sqrt{\gamma^2 x^2 + z^2}) \right) \quad (A.16) $$

where $A(\omega, M, \theta)$ represents the equivalent amplitude of the impulsive source:

$$ A(\omega, M, \theta) = \frac{i k_0}{c_0} \tilde{S}(k_0 \gamma^2 (\cos \Theta - M), k_0 \gamma \sin \Theta). \quad (A.17) $$
## 2. Source correction

At long range, we can assume that the acoustic waves have been launched at small angles, corresponding to $\theta \approx 0$. Denoting $A_M(\omega) \equiv A(\omega, M, \theta \approx 0)$, the source amplitude becomes:

$$A_M(\omega) = \frac{ik_0}{c_0} \tilde{S} \left( \frac{k_0}{1+M}, 0 \right). \quad (A.18)$$

### Case of a Gaussian spatial distribution

In this work, the impulsive source has a Gaussian spatial distribution

$$S(x, z) = S_0 \exp \left( -\frac{x^2 + z^2}{B^2} \right), \quad (A.19)$$

whose Fourier transform is:

$$\tilde{S}(k_x, k_z) = \pi B^2 S_0 \exp \left( -\frac{(k_x^2 + k_z^2)B^2}{4} \right). \quad (A.20)$$

The source amplitude in Eq. (A.18) is then given by:

$$A_M(\omega) = \frac{ik_0}{c_0} \pi B^2 S_0 \exp \left( -\frac{k_0^2 B^2}{4(1+M)^2} \right). \quad (A.21)$$

Note that for $M = 0$, the amplitude for the Gaussian impulse source in Eq. (A.21) leads to the same expression as given in Eq. (23) in Dragna et al. (2011) for the case of a homogeneous atmosphere at rest. From this expression we can also assess the frequency content of the Gaussian pulse. The maximum is reached for $f = c_0/(\sqrt{2\pi B})$ and a cut-off frequency $f_{10}$ can be defined such that the amplitude is equal to 10% of this maximum:

$$f_{10} \approx 0.6c_0/B. \quad (A.22)$$

The amplitude of the Gaussian impulsive source thus depends on the surrounding mean velocity. As an illustration, Figure 15 shows the source amplitude $A_M$ for $M = 0$ and for $M = 0.03$, which corresponds to the Mach number at the height of the wind turbine hub. The source amplitude tends to be shifted towards higher frequencies. The difference in dB is negligible at low frequencies. However, it becomes noticeable at high frequencies, even for this small Mach number. Thus, the difference reaches around 1.5 dB at 1 kHz.

### Evaluation of noise levels relative to the free field

The calculation of the noise levels relative to the free field from the LEE has to account for the modification of the source amplitude by the presence of the mean flow. For that, we estimate the free-field solution as:

$$\hat{p}_{ff}(\mathbf{x}, \omega) = -\frac{i}{4} A_M(\omega) H_0^{(1)}(k_0 R) \quad (A.23)$$

where the amplitude of the source is determined with Eq. (A.18) considering the Mach number at the source

<figure><img src="image_17.png" /><figcaption>FIG. 15. (a) Amplitude of the Gaussian impulsive source without mean flow and with a mean flow of Mach number equal to 0.03 and (b) difference in dB.</figcaption></figure>

height instead of $M = 0$. Eq. (A.23) is equivalent to Eq. (13).

Attenborough, K., Bashir, I., and Taherzadeh, S. (2011). "Outdoor ground impedance models," The Journal of the Acoustical Society of America 129(5), 2806-2819, doi: 10.1121/1.3569740.

Barlas, E., Wu, K. L., Zhu, W. J., Porté-Agel, F., and Shen, W. Z. (2018). "Variability of wind turbine noise over a diurnal cycle," Renewable Energy 126, 791-800, doi: 10.1016/j.renene.2018.03.086.

Barlas, E., Zhu, W. J., Shen, W. Z., Dag, K. O., and Moriarty, P. (2017a). "Consistent modelling of wind turbine noise propagation from source to receiver," The Journal of the Acoustical Society of America 142(5), 3297-3310, http://asa.scitation.org/doi/10.1121/1.5012747, doi: 10.1121/1.5012747.

Barlas, E., Zhu, W. J., Shen, W. Z., Kelly, M., and Andersen, S. J. (2017b). "Effects of wind turbine wake on atmospheric sound propagation," Applied Acoustics 122, 51-61, doi: 10.1016/j.apacoust.2017.02.010.

Berland, J., Bogey, C., and Bailly, C. (2006). "Low-dissipation and low-dispersion fourth-order Runge-Kutta algorithm," Computers & Fluids 35(10), 1459-1463, doi: 10.1016/j.compfluid.2005.04.003.

Berland, J., Bogey, C., Marsden, O., and Bailly, C. (2007). "High-order, low dispersive and low dissipative explicit schemes for multiple-scale and boundary problems," Journal of Computational Physics 224(2), 637-662, doi: 10.1016/j.jcp.2006.10.017.

Blumrich, R., and Heimann, D. (2002). "A linearized Eulerian sound propagation model for studies of complex meteorological effects," The Journal of the Acoustical Society of America 112(2), 446-455, doi: 10.1121/1.1485971.
Bogey, C., and Bailly, C. (2004). "A family of low dispersive and low dissipative explicit schemes for flow and noise computations," Journal of Computational Physics 194(1), 194-214, doi: 10.1016/j.jcp.2003.09.003.

Collino, F. (1997). "Perfectly matched absorbing layers for the paraxial equations," Journal of Computational Physics 131(1), 164-180, doi: 10.1006/jcph.1996.5594.

Cosnefroy, M. (2019). "Simulation numérique de la propagation dans l'atmosphère de sons impulsionnels et confrontations expérimentales. Numerical simulation of atmospheric propagation of impulse sound and experimental comparisons," Ph.D. thesis, Ecole Centrale de Lyon, 2019LYSEC014.

Cotté, B. (2019). "Extended source models for wind turbine noise propagation," The Journal of the Acoustical Society of America 145(3), 1363-1371, doi: 10.1121/1.5093307.

Dallois, L., Blanc-Benon, P., and Juvé, D. (2001). "A wide-angle parabolic equation for acoustic waves in inhomogeneous moving media: applications to atmospheric sound propagation," Journal of Computational Acoustics 9(02), 477-494, doi: 10.1142/S0218396X01000772.

Dragna, D., and Blanc-Benon, P. (2014). "Towards realistic simulations of sound radiation by moving sources in outdoor environments," International Journal of Aeroacoustics 13(5-6), 405-426, doi: 10.1260/1475-472X.13.5-6.405.

Dragna, D., Cotté, B., Blanc-Benon, P., and Poisson, F. (2011). "Time-Domain Simulations of Outdoor Sound Propagation with Suitable Impedance Boundary Conditions," AIAA Journal 49(7), 1420-1428, doi: 10.2514/1.J050636.

Dumortier, B., Vincent, E., and Deaconu, M. (2015). "Acoustic control of wind farms," The European Wind Energy Association Conference, Paris, France, November 17-20, 1-8, doi: 10.1121/1.47233730.

Gadde, S. N., and Stevens, R. J. A. M. (2021). "Interaction between low-level jets and wind farms in a stable atmospheric boundary layer," Physical Review Fluids 6(1), 014603, doi: 10.1103/PhysRevFluids.6.014603.

Gadde, S. N., Stieren, A., and Stevens, R. J. A. M. (2021). "Large-eddy simulations of stratified atmospheric boundary layers: comparison of different subgrid models," Boundary-Layer Meteorology 178(3), 363-382, doi: 10.1007/s10546-020-00570-5.

Gal-Chen, T., and Somerville, R. C. (1975). "On the use of a coordinate transformation for the solution of the Navier-Stokes equations," Journal of Computational Physics 17(2), 209-228, doi: 10.1016/0021-9991(75)90037-6.

Gilbert, K. E., and White, M. J. (1989). "Application of the parabolic equation to sound propagation in a refracting atmosphere," The Journal of the Acoustical Society of America 85(2), 630-637, doi: 10.1121/1.397587.

Hansen, K. L., Nguyen, P., Zajamšek, B., Catcheside, P., and Hansen, C. H. (2019). "Prevalence of wind farm amplitude modulation at long-range residential locations," Journal of Sound and Vibration 455, 136-149, doi: 10.1016/j.jsv.2019.05.008.

Heimann, D., and Englberger, A. (2018). "3D-simulation of sound propagation through the wake of a wind turbine: Impact of the diurnal variability," Applied Acoustics 141, 393-402, doi: 10.1016/j.apacoust.2018.06.005.

Heimann, D., Englberger, A., and Schady, A. (2018). "Sound propagation through the wake flow of a hilltop wind turbine-A numerical study," Wind Energy 21(8), 650-662, doi: 10.1002/we.2185.

Kayser, B., Cotté, B., Ecotière, D., and Gauvreau, B. (2020). "Environmental parameters sensitivity analysis for the modeling of wind turbine noise in downwind conditions," The Journal of the Acoustical Society of America 148(6), 3623-3632, doi: 10.1121/10.0002872.

Kayser, B., Mascarenhas, D., Cotté, B., Ecotière, D., and Gauvreau, B. (2023). "Validity of the effective sound speed approximation in parabolic equation models for wind turbine noise propagation," The Journal of the Acoustical Society of America 153(3), 1846-1854, doi: 10.1121/10.0017653.

Komatitsch, D., and Martin, R. (2007). "An unsplit convolutional Perfectly Matched Layer improved at grazing incidence for the seismic wave equation," Geophysical Journal International 72(1), 333-344, doi: 10.1190/1.2757586.

Lee, S., Lee, D., and Honhoff, S. (2016). "Prediction of far-field wind turbine noise propagation with parabolic equation," The Journal of the Acoustical Society of America 140(2), 767-778, doi: 10.1121/1.4958996.

Liu, L., and Stevens, R. J. A. M. (2020). "Effects of two-dimensional steep hills on the performance of wind turbines and wind farms," Boundary-Layer Meteorology 176(2), 251-269, doi: 10.1007/s10546-020-00522-z.

Nyborg, C. M., Fischer, A., Thysell, E., Feng, J., Søndergaard, L. S., Hansen, T. R., Hansen, K. S., and Bertagnolio, F. (2022). "Propagation of wind turbine noise: measurements and model evaluation," Journal of Physics 14, doi: 10.1088/1742-6596/2265/3/032041.

Ostashev, V. E., Wilson, D. K., Liu, L., Aldridge, D. F., Symons, N. P., and Marlin, D. (2005). "Equations for finite-difference, time-domain simulation of sound propagation in moving inhomogeneous media and numerical implementation," The Journal of the Acoustical Society of America 117(2), 503-517, doi: 10.1121/1.1841531.

Ostashev, V. E., Wilson, D. K., and Muhlestein, M. B. (2020). "Wave and extra-wide-angle parabolic equations for sound propagation in a moving atmosphere," The Journal of the Acoustical Society of America 147(6), 3969-3984, doi: 10.1121/10.0001397.

Propathopoulos, J. M., and Voutsinas, S. G. (2007). "Application of a ray theory model to the prediction of noise emissions from isolated wind turbines and wind parks," Wind Energy 10(2), 103-119, doi: 10.1002/we.211.

Rienstra, S. (2006). "Impedance models in time domain, including the extended Helmholtz resonator model," 12th AIAA/CEAS Aeroacoustics Conference (27th AIAA Aeroacoustics Conference), American Institute of Aeronautics and Astronautics, doi: 10.2514/6.2006-2686.

Sack, R. A., and West, M. (1995). "A parabolic equation for sound propagation in two dimensions over any smooth terrain profile: The generalised terrain parabolic equation (GT-PE)," Applied Acoustics 45(2), 113-129, doi: 10.1016/0003-682X(94)00039-X.

Salomons, E. M. (2001). Computational atmospheric acoustics (Springer Netherlands, Dordrecht).

Salomons, E. M., Blumrich, R., and Heimann, D. (2002). "Eulerian time-domain model for sound propagation over a finite-impedance ground surface. comparison with frequency-domain models," Acta Acustica 88(4), 483-492.

Sessarego, M., and Shen, W. Z. (2020). "Noise propagation calculation of a wind turbine in complex terrain," Journal of Physics: Conference Series 1452(1), 012063, doi: 10.1088/1742-6596/1452/1/012063.

Shen, W. Z., Zhu, W. J., Barlas, E., and Li, Y. (2019). "Advanced flow and noise simulation method for wind farm assessment in complex terrain," Renewable Energy 143, 1812-1825, doi: 10.1016/j.renene.2019.05.140.

Tian, Y., and Cotté, B. (2016). "Wind turbine noise modeling based on amiet's theory: effects of wind shear and atmospheric turbulence," Acta Acustica united with Acustica 102(4), 626-639, doi: 10.3813/AAA.918979.

Troian, R., Dragna, D., Bailly, C., and Galland, M.-A. (2017). "Broadband liner impedance eduction for multimodal acoustic propagation in the presence of a mean flow," Journal of Sound and Vibration 392, 200-216, doi: 10.1016/j.jsv.2016.10.014.

Van Renterghem, T. (2014). "Efficient outdoor sound propagation modeling with the finite-difference time-domain fdtd method: a review," International Journal of Aeroacoustics 13(5-6), 385-404, doi: 10.1260/1475-472X.13.5-6.385.