Consistent modelling of wind turbine noise propagation from source to receiver

The unsteady nature of wind turbine noise is a major reason for annoyance. The variation of far-ﬁeld sound pressure levels is not only caused by the continuous change in wind turbine noise source levels but also by the unsteady ﬂow ﬁeld and the ground characteristics between the turbine and receiver. To take these phenomena into account, a consistent numerical technique that models the sound propagation from the source to receiver is developed. Large eddy simulation with an actuator line technique is employed for the ﬂow modelling and the corresponding ﬂow ﬁelds are used to simulate sound generation and propagation. The local blade relative velocity, angle of attack, and turbulence characteristics are input to the sound generation model. Time-dependent blade locations and the velocity between the noise source and receiver are considered within a quasi-3D propagation model. Long-range noise propagation of a 5 MW wind turbine is investigated. Sound pressure level time series evaluated at the source time are studied for varying wind speeds, surface roughness, and ground impedances within a 2000 m radius from the turbine.


I. INTRODUCTION
As a result of increasing demand for renewable energy, fewer suitable land-based sites are available for wind farms. Given that noise is a primary obstacle for gaining broad public acceptance, the accurate noise assessment of wind turbines is a necessity. Accurate predictions of far-field wind turbine noise require the knowledge of the source levels and a realistic representation of the medium between the turbines and the receivers in which the sound propagation takes place. This is a complex task as both phenomena depend on a wide range of parameters.
The classical approach for far-field noise predictions assumes an overall source power level for a wind turbine dependent on its rotor diameter (RD) and uses a propagation relationship based on hemispherical spreading. Even though this approach neglects many physical processes, it has been the standard for some years. 1 A more advanced method is the Nord2000 that uses a semianalytical ray tracing model that models refraction effects using a linear approximation for the sound speed profile. 2 There are other corrections applied for undulating terrain and ground impedance. Even though this model was demonstrably more accurate than many other models, 3 there are certain shortcomings. For example, the source model is a monopole at hub-height used to represent the wind turbine irrespective of RD and the source strength is independent of the inflow conditions. The ray tracing method has been used for predicting noise from single wind turbines as well as wind farms in Refs. 4 and 5. Effects of various parameters, such as wind direction, ground impedance, and turbulence, were studied. Moreover, fullscale meteorological experiments and micro-scale models were used as input for the ray-tracing model and the effects of wakes were investigated. 6 Apart from the ray-tracing models, the parabolic equation (PE) method has also been used to address far-field wind turbine noise. Results obtained from the PE method were compared with various engineering models in Ref. 7 and suggestions for further research, such as realistic source representation and inclusion of terrain, were proposed. Flow fields obtained from Reynoldsaveraged Navier-Stokes simulations were incorporated with the PE method in Ref. 8 and the variation of the wake effects on far-field noise for different incoming shears were studied. A similar methodology was employed in Refs. 9 and 10 using large eddy simulation (LES) flow fields with unsteady source representation. While high-fidelity models used in these studies modelled the propagation physics more accurately, they did not consider the source characteristics in detail.
A commonly used method to model the wind turbine noise source was described in Refs. 11 and 12. The method divides the wind turbine blades into airfoil segments and sums the contribution of each segment's noise levels calculated using semiempirical relationships. This source model along with simple propagation calculations were used in Ref. 13 and the results were compared to near-field experiments. A similar source model was employed and coupled with the ray-tracing method in Ref. 14. Even though these models represented the source more accurately, the propagation effects were not well studied. a) Electronic mail: wjzhu@yzu.edu.cn The aforementioned models tend to focus on one side of the phenomenon, i.e., either a detailed source representation is used but propagation physics are disregarded or the opposite. Another major shortcoming of most modelling techniques is their steady nature. Listening room experiments [15][16][17] and dose-response relationship studies 18,19 showed a correlation between the source unsteadiness and the annoyance, particularly for wind turbines. Additionally, field experiments emphasized the considerable modification of far-field noise levels and characteristics under various atmospheric conditions. 20,21 A recent study showed qualitative results relating the wind direction and the number of complaints received from various wind farms. 22 They highlighted that the downwind receivers, who characterized the noise as thumping, were the most annoyed for distances greater than 2.5 km.
Although previous field experiments have underlined the effect of atmospheric conditions, the existing prediction tools do not consider the complicated flow around wind turbines and its effect on sound generation and propagation. There is a need for a model that includes the interaction between the incoming turbulent flow and the wind turbine, and the propagation physics consistently. The modular methodology proposed in this article attempts to fill the aforementioned gap. Here, the name "modular" refers to the fact that the models are loosely coupled such that the flow, source, and propagation models can be replaced with alternatives in the future. The technique uses a sound-generation model based on semiempirical modelling of airfoil noise within an aeroelastic tool and a propagation model based on the PE method. LES with an actuator line (AL) technique is employed for flow modelling. Using the LES flow fields as inputs to the aeroelastically coupled aeroacoustic simulations provides the ability to model wind turbine noise generation including its structural dynamics and interaction with the incoming turbulent flow. Although the aeroacoustic source models were obtained with the assumption of rigid bodies, the angle of attack change due to the flexibility of the blade is considered in the model. The constant change in the source levels and the spectral characteristics due to interactions between rotating blades, unsteady inflow, and turbine operational conditions can be taken into account. After each time step of the source calculation, each blade is represented with a single point source within the PE domains and independent simulations are carried out for each receiver. This approach ensures that the directivity, source strength, and realistic flow field between the source and receiver are taken into account. Using the timedependent flow fields and source locations and strengths, successive PE simulations are conducted. Subsequently, the nearand far-field noise levels of a 5 MW wind turbine for different incoming shear values and wind speeds are investigated.
The article is organized as follows. Section II describes the numerical methods and gives details about their implementation. Section III provides the modular methodology. Results are presented in Sec. IV, which is followed by the conclusions in Sec. V.

II. NUMERICAL METHODS
In this study, the National Renewable Energy Laboratory (NREL) 5 MW reference wind turbine is used. The wind turbine has a rotor diameter of 126 m and a hub-height of 90 m. The structural and aerodynamic characteristics of the wind turbine are detailed in Ref. 23 and the necessary inputs are released with FAST v8 simulation tool. 24 The controller is not used because the rotational speed of the wind turbine is kept constant at 7.8 rpm and 10.17 rpm for a 6 m/s and 9 m/s hub-height wind speed, respectively. The next three subsections detail the flow, source, and propagation models.

A. Flow model
The Technical University of Denmark's (DTU) pseudospectral incompressible Navier-Stokes solver is used to model the flow around the wind turbine. 25 The code solves the filtered continuity and momentum conservation equations in three-dimensional (3D) form with the Smagorinsky subgrid scale model. 26 The flow is driven by a pressure gradient wherein the turbulence is maintained solely because of shear stresses rather than buoyancy forces. This will also impact the sound propagation model as the temperature distribution, and thus the speed of sound is assumed to be constant. To capture the desired boundary layer characteristics, a wall model is applied at the bottom boundary. With this setup, the simulations are run until the boundary layer is established for a given hub-height wind speed (U h ) and a surface roughness value (z 0 ). After this and with the same initial conditions, two simulations are conducted with and without the turbine for each case listed in Table I.
The wind turbine rotor is modelled via the AL technique, 27 in which a body force distributed along the blade "lines" are used to represent the rotor blades. The AL technique is an efficient method to mimic the effect of wind turbines on the flow. It utilizes the tabulated airfoil lift and drag coefficients and the instantaneous velocity fields to determine the blade forces. Comparisons with field and wind tunnel experiments showed that the AL simulations can capture the wind turbine blade loading as well as the flow around it within an acceptable accuracy. 28,29 To avoid the downwind flow affecting the flow upwind of the wind turbine due to periodic boundary conditions in the horizontal direction, a buffer zone is applied to smoothly adjust the flow from the very-far-wake downwind to that of the precursor simulation. Figure 1 shows a snapshot of the streamwise velocity in case 1 interpolated onto slices at five azimuthal angles between the turbine and the selected receiver locations.

B. Aerodynamic noise source model
For the wind turbine noise source model, we use an aeroacoustics module based on NREL AirFoil Noise (NAFNOISE) code 30 within the FAST v8 modular framework. In this article,  we call this integrated code FAST v8þAA (aeroacoustics), which allows the modelling of aerodynamic noise generated by the blades in consideration of wind turbine structural dynamics and its interaction with the incoming turbulent flow. Using the blade element theory, each blade is divided into a number of two-dimensional (2D) airfoil elements. The total noise level at a given receiver location is predicted as the sum of the contributions from all the blade elements. This prediction method implicitly assumes that the noise generation mechanisms for all blade elements and all models are uncorrelated. Limitations of this assumption were investigated in Ref. 31, in which a correction method for 3D correlation was proposed. Considering that the span-over-chord ratio for each element is around 1.5 in the outer part of the blade, the error caused by this assumption does not exceed 2 dB for 20 Hz and decreases with frequency (negligibly small for frequencies above 100 Hz). Because it is rather complicated to calculate sound pressure levels (SPL) with 3D correlated formulae, we use the uncorrelated version in this study. Furthermore, it is worthwhile to mention that the convective amplification and Doppler effects are not taken into account in the present model. Figure 2 shows a snapshot of the flow field that the turbine is exposed to, two selected receivers, and the symbolic rays traced from the blade elements. For each element and time step, the local angle of attack and relative velocity are output from the aeroelastic solver. The changing angle of attack and the turbulence intensity are the main contributors of the noise modulation. The time-averaged angle of attack in the rotor area is shown in Fig. 3. The vertical wind shear is the dominant factor causing the angle of attack variation and case 1 has a more drastic change over one revolution than case 2. In addition, there is an asymmetry along the y axis because of a relatively small horizontal wind shear and the effect of the tower on the incoming flow.
In the present study, only two types of aerodynamic noise have been included: turbulent boundary layer trailingedge and turbulent inflow noise. The models and methods used for obtaining the necessary inputs are explained below.

Turbulent boundary layer trailing-edge noise (TBLTE):
The total SPL of noise generated from the interaction of the turbulent boundary layer with the trailing edge is calculated through the summation of the different contributions of noise on the pressure side, the suction side, and angle-dependent noise. These three noise mechanisms can be modelled semiempirically using scaling laws [see Brooks, Pope, and Marcolini (BPM) 32 ]. Different from the classical BPM noise generation model, the boundary layer characteristics in this study are obtained from Q 3 UIC (DTU's viscousinviscid interactive boundary layer flow solver), which is more accurate than XFOIL (a viscous-inviscid interactive code for subsonic isolated airfoils). 33 For a range of Reynolds numbers and angles of attack, the boundary layer thickness values are calculated and used as inputs to FAST v8. At each time step, each blade element's boundary layer thickness is interpolated for the corresponding airfoil, local angle of attack, and Reynolds number. Additionally, the angle of attack values used for the separation flag (a conditional switch used in the classical BPM model) is modified according to the lift coefficient curves of the blade airfoils used.

Turbulent Inflow Noise:
The turbulent inflow noise model developed in Refs. 34 and 35 is used for inflow noise estimation, including a correction for airfoil thickness proposed in Ref. 36. The turbulence intensity (TI) is calculated by taking into account the inflow characteristics obtained from the flow simulations (Sec. II A) and instantaneous airfoil locations. First, TI is calculated in the 2D plane at 1 RD upstream of the turbine for each 10-min data set and stored. This calculation is used as input to FAST v8 and then at each time instant the TI value at each airfoil location is interpolated and assigned. The integral turbulent length scale (L) is calculated using the relationship in Ref. 37 that can be expressed as a function of a nondimensional roughness parameter, a r , and height, z: L ¼ 25z 0:35 a À0:063 r . A more upto-date expression for length-scale distribution with height can be found in Ref. 38. However, in this study we consider only the given expression with the a r values 0.01 and 0.2 to represent open sea flat terrain and open country, respectively.
As explained in Sec. II A, the flow solver models only the shear-generated turbulence. This approach results in higher turbulence levels for higher shear cases. Figure 4 shows the averaged streamwise velocity and TI for cases 1 and 2. The corresponding SPL spectra for a ground-level receiver 2 RD downstream of the turbine are depicted in Fig.  5. It is observed that the low-frequency content (below 130 Hz) is dominated by the turbulent inflow noise. Therefore, case 1 levels are higher than those in case 2 within this frequency range. The other dominant noise source is from the TBLTE suction side. It dominates the midrange frequencies (between 200 and 800 Hz). The difference between the two cases in this frequency range is negligibly small. This observation is in line with Ref. 39. While the incoming shear is not that effective on sound generation, it has a significant effect on sound propagation, thus the far-field levels. 40 This effect will be investigated in detail in Sec. IV.

C. Sound propagation model
For the sound propagation modelling, we use DTU's WINDSTAR-PRO (wind turbine simulation tool for aerodynamic noise propagation). The tool implements a variety of 2D PE models in a FORTRAN environment with a message passing interface parallelization strategy in frequency and realization/time step. The code has been validated and used for various propagation calculations in Ref. 9. In this study, the two-dimensional, wide-angle, Crank-Nicholson, parabolic equation is used with the starter function and the implementation details given in Ref. 41. The method uses the effective speed of sound approach wherein the moving atmosphere is replaced by a hypothetical motionless medium with an effective sound speed where v x is the wind velocity component along the direction of propagation between the source and receiver. 42 In this study, the speed of sound is kept constant [c(x, z) ¼ 340 m/s] in the whole domain, and the propagation phenomena (i.e., refraction, diffraction, scattering) is solely due to the wind speed and its fluctuations around the wind turbine. This approach is representative of a neutrally stratified atmosphere on a day with high wind or thick cloud layers. 43 Figure 6 shows the time-averaged effective sound speed profiles at multiple locations in the midvertical plane obtained from the simulations listed in Table I. It is observed that each case has a distinct shear even though the hub-height wind speeds are the same for the first two and last two. Additionally, for the same surface roughness increasing the hub-height wind speed from 6 to 9 m/s results in increasing shear. It is clear that cases 3 and 4 have higher shear values than cases 1 and 2, respectively. While the upstream profiles are the main parameters for the upwind propagation characteristics, the wake evolution is the determining factor for the downwind propagation (see Sec. IV).
For the Crank-Nicholson parabolic equation calculations, the spatial resolution in both directions is set to oneeighth of the wavelength (Dx ¼ Dz ¼ k/8, where k is the wavelength of the considered frequency). Only flat terrain is considered and the ground impedance is characterized using the four-parameter model proposed in Ref. 44  respectively. These two ground covers were selected to define land-based and offshore conditions. The latter is important because the noise from offshore wind farms can be audible in long distances, as a result of hard ground and various atmospheric phenomena. 45 The other parameters of the impedance model are kept constant: pore shape factor (s p ¼ 0.75), Prandtl number (N Pr ¼ 0.72), grain shape factor (n 0 ¼ 0:5), porosity (X ¼ 0.3), ratio of specific heats (c ¼ 1.4), and density (q ¼ 1.19 kg/m 3 ).
All simulations are carried out for 1/3-octave band centre frequencies from 20 to 800 Hz, because the frequencies above have a negligible contribution (less than 0.1 dB for case 1 calculations) to the overall SPL due to atmospheric absorption as well as the dominant part of the spectra shown in Fig. 5. The corresponding sound pressure levels in each band are summed logarithmically to obtain the overall SPL: Here, N is the number of frequencies used and L p (f i ) is the sound pressure level defined as where the first two terms on the right-hand side (source power level and geometrical spreading) are obtained from the source model explained in Sec. II B for each receiver location. The third term represents the atmospheric absorption wherein the absorption coefficient is calculated according to International Standards Organization 9613-1 for air at 20 C with 80% relative humidity. The last term is the relative SPL that represents the deviation of a source from the free field as a result of ground and atmospheric effects. This last term is calculated using the PE method. Figure 7 shows a schematic of nine 2D PE domains from three blades to three receiver locations at one time instant. For each frequency and receiver, each blade is modelled as a monopole source, depicted with the red spheres in the figure.

III. MODULAR METHODOLOGY
This section is devoted to the explanation of the modular methodology that allows for a loose coupling of the three models explained in Sec. II.
(1) First, the flow field is simulated using LES. A 2D slice in spanwise and vertical directions (y-z slice) 1 RD upstream of the wind turbine is stored at each time step (0.02 s) to be used as input for the source simulations. The entire 3D flow field is also stored with a sampling frequency of 10 Hz to be used for propagation simulations. (2) The flow field sampled upstream of the turbine is used as input to FAST v8, forcing a fully aeroelastic turbine to a realistic atmospheric flow. While the time step used in the aeroelastic simulations is 0.02 s, the integrated aeroacoustic module is called every 0.1 s. The smaller time step for the aeroelastic simulations is due to the FAST v8 stability requirements. Additionally, we observe that a 10 Hz resolution for the aeroacoustic source simulations is sufficient because higher frequency content is calculated semiempirically.  Table I and for two ground covers (grassland and hard ground).
For each case, approximately 23 Â 10 6 independent propagation simulations are executed. The computational time required for sound propagation is 5000 processor hours. The flow solver requires 6000 processor hours after establishing the desired boundary layer. The stored flow fields are used for propagation simulations with various ground covers.
It is worthwhile noting that the SPL obtained from the successive PE simulations are evaluated at the source time. Even though the receiver time is the most common for noise measurements, it is disregarded in this paper. The term "time-dependent SPL" is used to refer to the PE output. Additionally, the simulations that follow the six-step procedure are referred to as "coupled simulations" from now on.   Figure 8 shows a snapshot of the overall SPL obtained from the coupled simulations for case 1. The contour plot is obtained via linear interpolation of the results on an equidistant grid with 20 m radial and 2 azimuthal spacing. It is observed that the noise directivity of the wind turbine is well captured, considering the crosswind levels are significantly lower than the upwind and downwind levels. This point will be elaborated on in Sec. IV. Additionally, it is also observed that even though the downwind and upwind levels have initially similar values, the atmospheric effects take over after a certain distance at which the upwind levels become much lower than the downwind ones.

IV. RESULTS
Using the methodology explained in Sec. III, three sets of calculations are carried out as described below. Note that, unless stated otherwise, the results show overall sound pressure levels summed from 20 to 800 Hz.

A. Time-averaged results
First, the output obtained from the SPL S is investigated. Figure 9 shows the time-averaged SPL for case 1. It is observed that the levels upwind and downwind of the turbine are considerably higher than the crosswind levels (approximately 25 dB). This pattern is caused by the wind turbine noise directivity and is observed for all cases listed in Table  I. However, this result differs from field experiments in Refs. 46 and 47, which show that the difference between the crosswind and downwind levels do not exceed 7 dB. The mismatch between the model output and the experiments is due to many reasons, such as background noise, binning of the wind direction (i.e., 610 ) to represent the crosswind results, and inaccuracy of the directivity functions used in calculations. Nevertheless, because the main focus of our work is the propagation effects on wind turbine noise, this large difference due to directivity is not further investigated.
Note that SPL S includes only the atmospheric absorption and geometrical spreading for the propagation. As a result, the sound pressure levels decrease with distance equally upwind and downwind of the turbine. This pattern changes when the atmospheric propagation effects are considered. Figure 10 shows the time-averaged output obtained from the coupled simulations (SPL SþP G ). Similar to Fig. 9, the SPL in the crosswind direction is less than that in the downwind or upwind directions. Different than Fig. 9, the upwind levels are lower than the downwind levels after a certain distance.
A more detailed comparison of all cases is depicted in Fig. 11, which shows the average upwind and downwind SPL SþP G for the grassland. It follows from the figure that the SPL at the first receiver location for cases 1 and 2 are approximately 9 dB lower than cases 3 and 4 in either direction. This is an expected result, because the rotational speed of the wind turbine increases with the increasing wind speed within this operation range. Subsequently, the tip speed increases from 51 to 66 m/s, which results in increased sound pressure levels. Also note that SPL SþP HG has a similar averaged footprint, though not shown here. Furthermore, for the same hub-height wind speed, it is observed that an increase in shear causes a slight increase (approximately 1 dB) in source levels at the closest receiver location, 2 RD away from the turbine. This result was discussed and attributed to the increased turbulent inflow noise in the low-frequency content in Sec. II B. Figure 11 shows the A-weighted SPL, which decreases the SPL considerably, because the low-frequency content dominates the overall levels (see Fig. 5). However, the trends of the weighted and unweighted SPL distributions with    Table I. From top to bottom: overall levels and frequency-dependent levels. In each subplot, the upper half domain is with the grass-covered land: D SPLG ¼ SPL SþPG À SPL S . The lower half domain is the hard ground: D SPLHG ¼ SPL SþPHG À SPL S . Wind direction is from left to right and the turbine is located at the center of the domain. distance are very similar. Because the focus of the study is on the propagation effects only, the unweighted results are reported further.
The small SPL differences between different shears in the near field vary significantly with distance depending on the atmospheric conditions. To isolate the propagation effects, the differences in SPL values are depicted in Fig. 12. D SPL G ¼ SPL SþP G À SPL S for the grassland case and D SPL HG ¼ SPL SþP HG À SPL S for the hard ground case are shown. In all the subplots, D SPL G is the top and D SPL HG is the bottom half domain.
We start with the common features of all the flow cases. First, it is observed that D SPL HG values are consistently higher than D SPL G . This is due to the more absorbing character of grassland than the hard ground. The SPL difference between two ground covers for the same flow case increases with distance. The difference values reach approximately 6 dB for the overall SPL at a 2 km distance.
Second, it is observed that the higher the shear, the earlier the shadow zones start upwind of the turbine [see  Fig. 12. Even though these zones are insonified because of diffraction and scattering (depending on the turbulence scales and acoustic wavelength 48 ), the timeaveraged levels of the coupled simulations are consistently lower than the source-only simulations.
A counterintuitive observation that can be deduced from Figs. 11 and 12 is that the sound pressure levels upwind of the wind turbine are higher than the downwind up to a certain distance. This distance varies with the considered case. For example, the upwind levels are higher than the downwind levels up to 1200 m for case 1, whereas this value is 950 m for case 4. This difference is due to the combined effects of refraction and ground reflection/absorption. The spectra at multiple locations are investigated to gain insight for the upwind and downwind SPL difference using one more set of simulations wherein the flow upstream of the turbine is used to simulate the downwind propagation. This means that the effects of the wake are neglected downwind of the turbine, and referred to as "no-wake simulations." Figure 13 shows the spectra for case 4. It is observed that for the upwind propagation, the ground caused dips in the near field are much less pronounced than both of the downwind ones. This is the main reason for the high upwind levels. As we get farther away from the turbine, the upwind levels start to decrease and eventually become lower than the downwind ones. Downwind of the turbine the levels obtained from the no-wake simulations are initially higher than the ones including the realistic wake flow. However, after a certain distance the opposite is true. This observation is in agreement with some of the results obtained in Refs. 6, 8, and 9. The main reason for this result is the wake-induced refraction behind a wind turbine, which plays a role similar to a SOFAR channel in underwater acoustics, 52 wherein the sound waves are ducted within the minima of the speed of sound. After a certain distance, waves are refracted downwards as the wake recovery takes place. This distance depends mainly on the considered flow case. For example, for case 1 the wake is not persistent enough because of the high incoming turbulence. Hence, a relatively narrow SPL amplification region (1300-1600 m) and low amplification levels are observed (approximately 1 dB). However, for case 4, in which the turbulence is low, the wake-induced SPL amplification levels reach up to 3 dB and the SPL amplification region is smeared to a larger area (1000-2000 m). The wake effects will be investigated with the time-dependent results in Sec. IV B. In Fig. 12(a), if we compare the upper and lower half, in the crosswind direction, the increased levels reach up to 8 and 4 dB for the hard ground and grassland, respectively. It is also observed that the regions where D SPL values are greater than or equal to zero, are longer in the crosswind direction than in the downwind direction for both ground covers (note that this is not valid in the downwind regions with wake-induced SPL amplifications explained previously). For the same flow and ground impedance case, the main difference between the D SPL values in the crosswind and downwind directions is caused by the flow field. Although sound waves propagate through the wake-induced flow field downwind of the turbine, the crosswind sound waves are affected only by the turbulent perturbations because the mean crosswind velocity is close to zero and we assume a neutral atmosphere in which there is no temperature gradient. This means that the refraction-related propagation effects play a larger role in SPL attenuation in the downwind case than in the crosswind case at certain distances. This outcome does not necessarily mean that the neutral atmosphere causes less attenuation. On the contrary, it highlights the complexity of the propagation phenomena. Nevertheless, because the wind turbine noise is dominated by a dipole emission pattern the overall levels are significantly lower in the crosswind direction than in the downwind direction.

B. Time-dependent results
In this study, the quantification of the wind turbine noise amplitude modulation (AM) is done using the method proposed by the UK Institute of Acoustics Noise Working Group on Wind Turbine Noise Amplitude Modulation. The details of the method can be found in Ref. 53. Only a brief description of the calculation procedure is given here. The method is based on transforming a SPL time series of 10 s blocks into the frequency domain to detect the blade passage frequency and its next two harmonics. Afterward, some threshold checks are applied (i.e., a harmonic is kept only if its reconstructed time signal has a peak-to-trough ratio bigger than 1.5 dB). The harmonics that pass the threshold checks are used for the conversion from frequency to time domain. Figure 14(a) shows a sample 10 s block and Fig.  14(b) shows the power spectrum that is calculated by taking the square of the absolute discrete Fourier transform (DFT), F(x), of the time signal and dividing it by the square of number of points. Figure 14(c) shows the real part of the DFT for all frequencies (white bars) as well as the ones that are included for reconstructing a new signal after the aforementioned threshold checks (red bars). Figure 14(d) shows the original and the reconstructed signal. The modulation depth is then determined by subtracting the fifth percentile (L 95 ) of the reconstructed signal from the 95th percentile (L 5 ). Dashed lines in Fig. 14(d) show these percentiles. This method is applied to the three cases described earlier and the AM levels are referred to as AM S , AM SþP G , and AM SþP HG . Different from the Institute of Acoustics method, the detected AM levels in each 10 s block are averaged over a 10-min simulation. Figure 15 shows the AM levels of the source-only and coupled simulations with grassland for case 3. The sourceonly simulation captures the near-field AM in the crosswind direction. This is relatively well understood 54 to be caused by the directivity and the blade rotation. The coupled simulation, on the other hand, shows a different trend. The results indicate a significant increase in AM levels when the propagation effects are taken into account in the far-field downwind or upwind of the turbine.
If we look closely at some of the time series (see Fig.  15) certain observations can be deduced. The time series obtained from receiver number 1 (downwind 1260 m) shows that the relatively small modulations observed from the source-only simulation are enhanced with the propagation effects. This increase in the downwind AM levels can be explained with the random phase and amplitude fluctuations caused by the turbulence that is only considered with the coupled simulation. It is worth noting that the turbulence scales smaller than the LES grid size are not resolved, and thus the full scattering phenomenon might not be captured, especially at the highest frequencies considered. Apart from the atmospheric and wake-induced small-scale turbulence, another reason could be the overall unsteady behaviour of the wind turbine wake (also known as wake meandering 55 ). This may cause the sound waves refract through full or partial wake, which results in increased SPL fluctuations at the selected downwind locations. For a better understanding of this phenomenon, the top view of the instantaneous streamwise velocity in case 3 at two different time instants are shown in Fig. 16. It is clear that the wake has a dynamic nature. Note that, the time scale of the large-scale wake movement is much longer than the blade rotation time scale. Therefore, at each instant sound waves emitted from each source propagate through different refractive regions. It is hard to distinguish which effect is the main contributor of the enhanced downwind AM. However, a comparison of the AM levels obtained from the simulations with and without the wake (see Fig. 17) emphasizes the importance of the wind turbine wake.
The time series of receiver number 2 (upwind 1260 m) in Fig. 15 shows an overall SPL attenuation for the coupled simulation. However, the modulation depths are considerably higher than in the source-only simulation. One reason for this could be the constant change of the upwind extent of shadow zones caused by the blade rotation. For the sake of argument, let us consider the idealized upwind propagation without turbulence. It is well known that for a low-elevation source the shadow zones start within a shorter distance than an elevated source. 49 This means that under idealized conditions, the start of the shadow zone is determined only by the highest source. In this study, each blade is represented with one source within the propagation model and the source is often located around the tip of each blade (see Sec. II B). Two time instants are selected during the turbine rotation: A and B. At instant A, one source goes through the top tip height (153 m), whereas the other two sources are at a 58 m height. This means the shadow region start is dominated by the source at 153 m. At instant B, one source goes through the bottom tip height (28 m) and the other two sources are located at a 117 m height. This is the instant when the highest source among the three is at the lowest height. Therefore, the shadow zone will start closer to the turbine than in instant A. Figures 18(a) and 18(b) show the contour plots of the logarithmically summed SPL for 800 Hz obtained from three independent 2D PE simulations for this idealized case. If we measure the noise around the start of a shadow zone, the amplitude fluctuations would be significantly larger (i.e.,  at instant A we hear something but at instant B we do not hear anything). Figure 19 shows the time signals at multiple locations. The AM levels calculated around the shadow zone boundaries [Figs. 19(b) and 19(c)] are much higher than the ones before [ Fig. 19(a)] or after [ Fig. 19(d)] the shadow zone. The idealized situation changes when the turbulence is accounted for, because it is not possible to determine a distinct boundary of the shadow zone. Nevertheless, with the rotation of the blades and because of the turbulent flow field between the sources and receivers, these regions constantly change shape (see Mm. 1) resulting in increased upwind AM levels. An important caveat is that these simulations do not model the background noise. It is very likely that the attenuated levels at the upwind would be lower than the background noise in field measurements. Therefore, it may not be always possible to measure these SPL fluctuations.
Mm. 1. Time evolution of overall SPL at a 2m receiver height for flow case 1 obtained from the combined generation and propagation simulations at source time. This is a file of type "mov" (31.1 MB).
A more detailed study is shown in Fig. 20. Similar to the Sec. IV A, to isolate the propagation effects on AM, the difference contour plots are depicted (D AM G ¼ AM SþP G À AM S for grassland and D AM HG ¼ AM SþP HG À AM S for hard ground). It is observed that in case 3 the region with high upwind AM is the widest, namely, from 500 to 1900 m. The same thing is valid for the downwind AM. Additionally, the case 3 levels are higher than those in all other cases. The comparison of case 1 and 2 shows that case 1 has narrower AM regions with higher levels than case 2. Figure 20 shows that the shear, wake flow evolution, and ground characteristics play a role in AM levels both upwind and downwind of the turbine.

V. CONCLUSION
In this article, a consistent and modular modelling technique for wind turbine noise propagation was developed. A flow solver based on LES/AL, a wind turbine noise generation model based on FAST v8/NAFNoise, and a propagation model based on the PE method were loosely coupled. The unsteady flow fields were used for both noise generation and propagation calculations. The focus of the present study was on the propagation effects for varying wind shear, wind speed, and ground covers. Various phenomena, such as shadow zone existence, absorbing character of grass versus hard ground were captured when the propagation effects were included in the noise calculations. A clear increase in the source level was observed with increasing wind speed (9 dB difference at 2 RD away from the turbine, from U h ¼ 6 m/s to U h ¼ 9 m/s). It was shown that for the same hub-height wind speed, increased shear resulted in a small increase (approximately 1 dB) in the source levels. However, small SPL differences between various cases in the near field were enhanced in the far field as a result of the propagation effects. This observation was valid both for the steady and unsteady investigations. The SPL time signals evaluated at the source time showed enhanced far-field AM levels both upwind and downwind of the turbine.
The developed modular technique for the coupled model in this article allows for substitution of each part of the model with an alternative. This approach can be a good way to overcome the limitations, such as inaccurate airfoil noise calculations, simple wind turbine representations within the flow, or propagation models. Although the submodels are easily replaceable with their higher fidelity versions, more accurate modelling of the flow around wind turbines in complex terrain is one of the major challenges that lies ahead for accurate noise prediction.