An inter-model comparison of parabolic equation methods for sound propagation from wind turbines

: The modeling of sound propagation for land-based wind turbines is a complex task that takes various parameters into account. Not only do the wind speed and wind direction affect the noise received at a certain position by changing the refraction of the sound, but also the terrain complexity, ground impedance, and receiver position relative to the source and ground all affect propagation. These effects are seen by the reﬂections of the sound at the ground surface causing interference of sound waves, or by the receiver being positioned in and out of noise shadow zones in the upwind far ﬁeld position, or in steep terrain irregularities. Several sound propagation models with different levels of ﬁdelity have been developed through time to account for these effects. This paper will focus on two different parabolic equation models, the Beilis-Tappert Parabolic Equation and the Generalized Terrain Parabolic Equation, through theoretical studies of varying terrain complexity, ground impedance, and sound speed proﬁles (upwind, downwind, and no wind). In addition, the propagation models are validated through spectral comparisons to noise measurements from two different campaigns considering loudspeaker noise and wind turbine noise, respectively. V C


I. INTRODUCTION
With the increasing demand for onshore wind energy, the annoyance caused by noise emitted from wind turbines is of high concern. Therefore, reliable methods for predicting noise levels at residences neighbouring onshore wind farms are needed. Parabolic Equation (PE) models (Salomons, 2001) have been widely used for determining sound propagation in underwater as well as atmospheric media. Applying a PE model rather than a simple model such as the ISO 9613-2 (ISO, 1997) standard model, yields the possibility of applying detailed sound speed and terrain information to the problem of sound propagation (Dallois et al., 2001;Ostashev et al., 2020). This makes the PE models highly suitable for investigations of the propagation of wind turbine noise.
The Beilis-Tappert Parabolic Equation (BTPE) model (Beilis and Tappert, 1979) shares features with the Greens Function Parabolic Equation (GFPE) method (Salomons, 1998), making it a computationally more efficient propagation model than the Generalized Terrain Parabolic Equation (GTPE) model (Sack and West, 1995). In both models, it is noticed that the grid spacing depends on the wavelength of the considered frequency, k. Thus, for higher frequencies, the models become more computationally expensive. Since the BTPE allows for larger horizontal grid spacings than the GTPE (Jensen et al., 2011;Parakkal et al., 2012), the computational time of the GTPE is substantially longer than the BTPE. A comparison between sound propagation models in varying terrains and refracting atmospheres is thus considered warranted and can be accomplished by including downward and upward refraction as well as different classes of ground flow resistivity. Thus, by varying these parameters in a controlled numerical setup, better knowledge about the modeling differences can be obtained.
Previously, several studies focusing on sound propagation from point sources closer to the ground and on wind turbine sound have been presented. The effect of the ground flow resistivity on sound propagation has previously been studied through measurements over various ground types by € Ohlund and Larsson (2015) and Conrady et al. (2018) as well as numerically by Kayser et al. (2019) supported by laboratory measurements. In addition, Dragna et al. (2015) compared the validity of different ground impedance models and Van Renterghem et al. (2014) modeled the offshore acoustic propagation at different sea states by varying the impedance.
The atmospheric conditions as well have a major impact on sound propagation and the effects were studied in relation to sound propagation from wind turbines by Barlas et al. (2018). In an upward refracting atmosphere, shadow zones are likely to appear at a certain distance from the source depending on the frequency. Evans and Cooper (2012) presented measurements of the propagation loss difference between upwind and downwind propagation, while Bolin et al. (2020) provide numerical studies of the shadow zone occurrence and influence of turbulence in upward refracting conditions. It is commonly known that in a downward refracting atmosphere, the sound rays will be bent towards the ground. This can, on one hand, cause a higher perception of noise for a receiver located in the downwind position of a wind turbine. On the other hand, it can result in more reflections on the ground leading to larger negative interference effects. Like for the upward refraction atmosphere, shadow zones can appear in the case of a downward refracting atmosphere as well depending on the frequency and altitude of the source along with the atmospheric conditions.
The effect on the GTPE model from the wind turbine wake and diurnal effects has been investigated through thorough studies by Barlas et al. (2018Barlas et al. ( , 2017b. These studies include high fidelity modeling with large eddy simulations and an unsteady representation of the noise sources at each blade tip. Studies of moving point sources distributed along the rotor have also been presented by Cott e (2018, 2019). Lee et al. (2016) similarly analysed the wake effects of a wind turbine on sound propagation by the use of an Actuator Disc source model and the Crank-Nicholson PE model (Gilbert and White, 1989) as well as compared the PE results in the upwind and downwind position to field measurements. Furthermore, the implementation of a turbulent velocity field, causing scattering of sound into the shadow zones, has been presented and analysed by Ostashev (1994) and Wilson et al. (1996). Since the use of higher fidelity flow modeling or turbulence in the PE models increases the computational time and complexity significantly, these approaches are left out of this study as they are considered too complex for inter-model comparisons. Therefore, only instantaneous computations with a steady point source and no turbulence are considered.
By the use of a moving source approach, Kayser et al. (2020) studied the sensitivity of the wide-angle parabolic equation to different parameters and found that the angle between the direction to the receiver and the wind direction as well as the ground flow resistivity have significant influences on the modeling results. Moreover, the number of earlier comparisons of atmospheric PE modeling results for complex terrain are limited and are often only presented for one type of sound speed profile.
Additional numerical studies of PE models have been conducted by Almquist et al. (2014) and Parakkal et al. (2012), which will also build the foundation for the studies in this paper. These studies, as well as the studies provided by Attenborough et al. (1995), are limited to a source height of 5-10 m, which is naturally not applicable to wind turbines.
In this study, the two well-known PE models, the GTPE model (Barlas et al., 2017a;Salomons, 2001) and the BTPE model (Jensen et al., 2011;Parakkal et al., 2012), are applied to the study of sound propagation from an elevated source at a height representative of the hub of a modern day wind turbine. This study investigates the effects of the terrain complexity and impedance on sound propagation in both no wind, upward and downward refracting atmospheres. The Jeltsch energy-conserving PE (JEPE) model (Karasalo and Sundstr€ om, 1996) is used as a third algorithm through the numerical cases. The present study is further supported by comparisons to two experimental campaigns. Previously, comparisons between the GFPE model and flat terrain (over water) measurements have been done (Bolin et al., 2014), while the GTPE model has been evaluated and compared to noise propagation measurements by Cao et al. (2022) and Nyborg et al. (2022). The measurement campaigns included in this paper consider the measured propagation of noise from a loudspeaker attached to the nacelle of a wind turbine in relatively flat terrain , and the measured noise of a wind farm in undulating terrain .
The PE models are thereby compared in various terrain and meteorological conditions providing a study including both numerical and experimental cases. The purpose of this study is thus to give a better understanding of the relative accuracy of using one PE model over the other and to examine the sensitivity of each model to sound speed profile and terrain.
This paper is organized as follows: Sec. II briefly introduces the implemented sound propagation models in question and highlights their main differences. Section III presents the results from the numerical cases by varying the refraction, ground impedance, and ground complexity in the different PE models. Section IV presents the two measurement campaigns: The loudspeaker campaign A, in which the results are evaluated through a statistical approach using 2 h of measurement data, and the wind farm campaign C in which the 10-min equivalent sound pressure level is used for two different weather scenarios. Last, an overall discussion and the main conclusions are given in Sec. V.

II. SOUND PROPAGATION MODELS
This section describes the different sound propagation models analysed in this study which are all originating from the PE method (West et al., 1992). The different models share similarities in solving a system of parabolic equations with the use of a square-root operator approximation. Each PE method introduces a different approximation to the square-root operator and uses different numerical methods to solve the system of equations. Thus, the PEs used in this study results in differences both theoretically and numerically.

A. GTPE
The GTPE is based on the wide angle parabolic equation (WAPE) method, which uses the Crank Nicholson approach to solve the parabolic equations system (Gilbert and White, 1989;Sack and West, 1995;Salomons, 2001). The GTPE models sound propagation over a complex terrain by introducing a coordinate shift at changes in the terrain height in the Helmholtz equation prior to solving the system of equations. Thus, the method "flattens" the terrain in the computations such that the vertical z coordinate is given relative to the complex terrain.
The GTPE method has been implemented in FORTAN using a message passing interface (MPI) (Barlas, 2017;Barlas et al., 2017a;Shen et al., 2019). A limiting elevation angle of 30 is obtained by the suggested 2nd order Gaussian starter function (Salomons, 2001). The grid resolution in both the range-and vertical direction is determined through a grid convergence study, and set to Dz ¼ Dr ¼ k=8 for all frequencies. The GTPE method assumes a motionless atmosphere with an effective speed of sound, expressed as c eff ¼ c 0 þ u, where c 0 is the adiabatic speed of sound and u is the wind speed field projected into the plane of propagation. For the upper boundary, an artificial absorbing layer with a thickness of 50k is assumed to avoid reflections of the sound back into the computational domain. The height of the computational domain to the absorbing layer has been set to 500 m for all computations in this paper. The effect of choosing larger computational heights of 1 km and 2 km was examined as part of the analysis, but did not show any effects on the results.

B. BTPE
Faster computational times compared to the GTPE can be expected for Fourier-based parabolic equation solvers as longer-range steps can be chosen and the Fourier-based algorithms are not computationally expensive. The BTPE model (Beilis and Tappert, 1979) has many similarities to the method of the GFPE (Gilbert and Di, 1993;Salomons, 1998Salomons, , 2001 and implements a split-step Fourier technique (Jensen et al., 2011), which has been adjusted for atmospheric sound propagation treating irregular terrain by a "phase screen" (Parakkal et al., 2012). The BTPE has the advantage over the GFPE that it handles irregular terrain (i.e., range-dependent ground altitude). In this paper the BTPE, as described in Parakkal et al. (2012), has been implemented in the software MATLAB (ver 2019b) and validated to the results presented in Parakkal et al. (2012). The range steps, Dr, have been chosen to 2:5k for frequencies at and below 50 Hz, 5k between 63 and 100 Hz, and 10k for frequencies of 125 Hz and above; the vertical discretization is Dz ¼ k=10 throughout and results at positions have been interpolated from these grids.
Where the GTPE uses the 2nd order Gaussian starter function, the BTPE applies a 4th degree polynomial starter function as suggested by Salomons (2001). The computational domain heights are set to 500 m for ranges up to 2000 m and 1000 m for ranges up to 5000 m. A virtual absorption layer at the upper boundary similar to the one suggested for the GFPE by Salomons (2001) has been implemented.

C. JEPE
The JEPE (Karasalo and Sundstr€ om, 1996) is based on formulating the PE as an ordinary Differential-Algebraic (DA) system of equations by discretizing the vertical coordinate and solving the DA-system by a fourth-order strongly damped scheme (Jeltsch, 1977). The JEPE was originally developed for underwater acoustic applications and later adapted to atmospheric sound propagation by introducing a local artificial absorption layer (Salomons, 2001) as an approximation of the non-reflecting condition at the upper boundary. Any irregular terrain is handled by using orthogonal boundary-adapted curvilinear coordinates (Abrahamsson, 1991), which simplifies the formulation of stable boundary conditions at the upper and lower boundaries. An adaptive range-step controlled by a local error estimate is used for the marching in range.

III. NUMERICAL CASES
This section briefly describes the parameters that are varied through the numerical case studies. The parameters include the atmospheric refraction, terrain, and ground impedance. In the last part of the section, the computational results of the GTPE and BTPE in the different conditions are presented and discussed. For these theoretical cases, results using the JEPE model (Karasalo and Sundstr€ om, 1996) were available and are therefore included as well.

A. Atmospheric boundary layer
The BTPE and the GTPE both assume a hypothetical motionless atmosphere with an effective speed of sound, consisting of the sound speed in a non-moving air, c 0 , and the wind speed profile in the sound propagation direction. Thus, the vector properties of the wind velocity are omitted and instead approximated as a scalar at each grid point. The effective speed of sound is determined by the logarithmic profile as suggested by Salomons (2001), with c 0 ¼ 340 m/s, the roughness length z 0 ¼ 0:1 m, and the refraction parameter b ¼ 61 for a downward and upward refracting atmosphere, respectively, and b ¼ 0 for a no-wind atmosphere. A logarithmic wind speed profile is expected to be representative of the considered scenario and is commonly used in boundary layer meteorology (Stull, 1988). Figure 1 shows the considered sound speed profiles. In the cases of irregular terrain, the sound speed profiles are following the ground elevation. Therefore, at a receiver height of 2 m above the ground surface, the speed of sound will be constant with the range, r. This is a simplified behaviour of the flow field compared to reality where the wind interacts with the ground causing instabilities in the flow field depending on the stability of the atmosphere (Finnigan et al., 2020). The altered flow field resulting from these interactions can, for example, be evaluated by large eddy simulations (LES) (Liu and Stevens, 2020). However, to keep the studies simple, this is considered out of scope of the presented work. Turbulence and temperature variations in the atmosphere directly affect the speed of sound profile (Barlas et al., 2018;Barlas et al., 2017b;Salomons, 2001), and both the GTPE and BTPE have the possibility of taking turbulence into account. Statistical approaches to include turbulence have previously been presented by Ostashev (1994) and Wilson et al. (1996), while Cott e and Blanc-Benon (2007) and Ostashev and Wilson (2000) investigated the contributions of the turbulence spectrum to acoustic scattering. By considering turbulence in a sound propagation model significant scattering of the sound field can appear. This can especially be of importance for upwind propagation where higher sound levels can occur in the shadow zone. When increasing the acoustic frequency, the shadow zone caused by refraction appears at a distance closer to the sound source. However, when computing sound propagation from a wind turbine, the high-frequency sound decays faster due to the large contribution from atmospheric absorption. When the sound source characteristics are further considered, the energy of the source is highest at lower frequencies. The contribution of the high-frequency noise to the overall sound level at receivers in the far field of a wind turbine is therefore of lower importance. On the other hand, when considering complex terrain, significant shadow zones can arise at positions behind steep hills. Considering turbulence in the computations over complex terrain would likewise result in scattering of sound into these shadow zones. This can be of higher importance when modeling the noise at receivers positioned in the near field of a wind turbine. Moreover, temperature fluctuations in the computational domain affect the adiabatic sound speed used in the determination of the effective sound speed. However, for simplicity and better reproducibility in the presented studies, the turbulence and temperature variations are not included.

B. Terrain complexity
To explore the effects of terrain complexity on the PE models, this study compares propagation over flat terrain, a Gaussian hill of increasing steepness (Parakkal et al., 2012), and an undulating terrain with four continuous hills (Almquist et al., 2014). The terrains have previously been studied by Barlas (2017), where both the GFPE and WAPE were compared. However, while the study by Barlas (2017) focused on a single low frequency and low source heights, this paper extends the study to include full range analysis for 50 and 500 Hz as well as spectral analysis ranging from a frequency of 20 Hz to 2 kHz at a few selected distances. This paper further considers the analysis of the effects from different ground impedance classes (Sec. III C) and three different sound speed profiles (Sec. III A). Furthermore, the observer angle for a receiver positioned 2 m above the ground in the different terrain scenarios is checked. The observer angle is defined as the angle between the ground and the straight line from the receiver position to the sound source. As the focus of this paper is kept at elevated sources for comparability with the hub height of a modern wind turbine, only results with the source at a height of 100 m are included. The elevated source height of H ¼ 100 m leads to larger elevation angles and since the BTPE and GTPE models both have limiting elevation angles, ground distances shorter than 200 m are considered invalid.
Two examples of the terrain following speed of sound in a downward refracting atmosphere are extracted from the GTPE and BTPE computations and shown in Fig. 2. From these figures, it is evident that the resolution of the complex terrain varies from model to model. This is a result of the large range step in the BTPE, which causes the undulating terrain in the BTPE (The right panel of Fig. 2) to be generally coarser resolved than in the GTPE (The left panel of Fig. 2).

C. Ground impedance
Three different ground impedance classes, listed in Table I, are chosen to investigate the sensitivity of altering the ground type at the lower boundary of the computational domain. The ground impedance classes represent a very soft, snow-like ground, an uncompacted ground, and a compacted dense ground. An acoustically hard surface like the compacted dense ground with higher sound reflecting properties can also be related to the impedance of water. The ground flow resistivities of the three ground categories are listed in the table that follows. To compute the frequencydependent ground impedances from the ground flow resistivities, the four parameter impedance model by Attenborough (1985) is used. The remaining parameters in the model are kept constant and are given by: the pore shape factor, s f ¼ 0:75, the Prandtl number, N Pr ¼ 0.7, the grain shape factor, n 0 ¼ 0:5, the porosity, X ¼ 0:3, the specific heat rate of air, c ¼ 1:4, and the ground density q ¼ 1:19 kg/m 3 . It should be noted that the ground acoustic impedance has been computed separately and used as an input to the different PE models to ensure that a similar lower boundary is used in the numerical cases. In the comparisons to loudspeaker and wind turbine noise measurements done later in this paper, the ground acoustic impedance is obtained in the PE models from the specified classes of ground flow resistivity.

D. Numerical results
The frequency depending sound pressure level, L p ð f i Þ, at a receiver is normally partly evaluated by the sound power level of the source and partly by the propagation of sound at a distance, d, from the sound source (Salomons, 2001;Wagner et al., 1996), In the definition of L p the frequency-dependent sound power level of the sound source is given by L W ð f i Þ, the atmospheric absorption is denoted by að f i Þd, the geometrical spreading by 10 log 10 ð4pd 2 Þ and the sound pressure level relative to the free field sound pressure level by DL p ð f i Þ. The numerical cases consider the comparison of the transmission loss (TL) from the different sound propagation models including the geometrical spreading and DL p , but excluding the atmospheric absorption. For the comparisons with measurements from the loudspeaker and wind turbine campaigns, atmospheric absorption is considered. For the theoretical cases, the sound power level at each frequency band, L W ð f i Þ, is set to 0 dB and the atmospheric attenuation (ISO, 1993) is neglected. Thus, the following analysis considers the TL as the effects caused by the geometrical spreading of the sound and the relative sound pressure level caused by meteorological and ground effects.

Range-dependent TL
Figures 3 and 4 present the range-dependent TL over a flat terrain obtained from the BTPE and GTPE models. The computations are given for frequencies of 50 and 500 Hz, and the different sound speed profiles and classes of ground flow resistivity. The increasing ground flow resistivity only has a minor effect on sound propagation for f ¼ 50 Hz. Moreover, due to the increasing reflecting properties of the ground, the TL is lower for r ¼ 250 kPa s/m 2 and lowest for r ¼ 2000 kPa s/m 2 yielding higher perceived noise levels in these cases. For f ¼ 50 Hz, the models have a slight offset of 1 dB but generally agree well. It is seen that the TL offset obtained between the GTPE and the BTPE and JEPE at no wind and downwind propagation becomes slightly more obvious for the higher ground flow resistivities. The PE models start to deviate in the upward refracting conditions at f ¼ 50 Hz at longer ranges which could be caused by numerical noise at very low dB values. One should however keep in mind, that the TL values in these regions are unrealistically low due to the omission of turbulence. The higher frequency of f ¼ 500 Hz introduces significant deviations in   the models by the interferences observed in the near field in Fig. 4. Comparing the BTPE and GTPE in downward refraction, these interferences between the direct and reflected sound have a phase shift that increases with ground flow resistivity. The differences in the interferences could be a result of the higher dispersion error of the PE models at large elevation angles. Thus, combined with a high ground flow resistivity, the interference patterns are amplified causing them to appear at further distances from the source. Furthermore, it should be noted that the different PE models are considered far field models, which makes the results in shorter distances to the source questionable. It is observed how the JEPE does not capture these interferences, further emphasizing that they are caused by the starting conditions of the BTPE and GTPE rather than realistic ground reflections. In addition, the shadow zone caused by upward refraction appears earlier for f ¼ 500 Hz, especially causing large deviations between the JEPE and the remaining models early on. The effects of increasing the ground impedance are emphasized when introducing terrain complexity. Surface plots of the BTPE and GTPE models are presented for a ground flow resistivity of r ¼ 250 kPa s/m 2 at 50 Hz in a downward refracting atmosphere in Fig. 5 and an upward refracting atmosphere in Fig. 6. It is seen how both the upward and downward refracting cases introduce interference patterns reaching beyond the height of the source position in the computational domain. The interferences consistently start at the front of a hill and are repeated for each hill along the range of the domain in the downward refracting case, whereas for the upward refracting case, a shadow zone caused by the refraction is appearing after the second hill, making the interference less prominent in the further range. It should be noted, that the use of stationary computations that omit the presence of atmospheric turbulence introduces the interference patterns observed in Figs. 5 and 6. Thus, in a scenario where the time-varying sound is considered, the phases of the interacting sound waves and thereby the interference patterns will vary.
In the range-dependent TLs for both 50 and 500 Hz in Figs. 7 and 8, the observations of interference peaks are supported by the distinct destructive interferences observed especially for downward refraction. The JEPE, BTPE, and GTPE models are all capturing these destructive interferences, however, with a phase shift between the GTPE and the remaining models at higher ground impedance. The phase shift is larger for 50 Hz, suggesting it to be caused by the significant difference in grid resolution, since a coarser grid results from the low frequency. It is further noticed how the receiver (2 m above ground) moves in and out of the shadow zone in the case of upward refraction at 500 Hz as a consequence of the earlier appearance of the shadow zone and the increasing terrain height. In flat terrain areas, an  obvious distinguishing between positions inside and outside of the shadow zone can be made. In the case of complex terrain, it is evident from Fig. 8 that special care should be taken before concluding if a receiver is in the shadow zone at certain distances from the source. Thus, one receiver closer to the source but at lower ground levels may be in a shadow zone, while a receiver further away but at higher ground levels can consequently be outside of the shadow zone. Furthermore, neglecting turbulence causes high uncertainties in the shadow zone, since the scattering of noise is not included (Wilson et al., 1996). Hence, by including the turbulence in the computations, larger sound levels are expected while the random behaviour of the velocity field in a turbulent atmosphere is expected to yield better agreements between the different PE models.
In Figs. 9 and 10, the difference in range-dependent TL, DTL ¼ TL GTPE À TL BTPE , is given for the four Gaussian hills, the three different classes of ground flow resistivity, and the different sound speed profiles. The computations are performed for frequencies of 50 and 500 Hz, respectively. In order to subtract the BTPE obtained TL from the TL computed by the GTPE, the GTPE is interpolated to the coarser grid of the BTPE. For the lowest ground flow resistivity, r ¼ 12 kPa s/m 2 , it is observed how the models agree well at the front and top of the hills while they deviate behind the hills in all refracting conditions. The steeper the hill the larger the deviation observed just behind the hill. This implies that the models are treating the potential shadow zone appearing behind the hill differently. The size of this deviation is consistent with all speed of sound profiles, suggesting that it is caused exclusively by the terrain effects in the different PE models. The deviations could further be a result of numerical noise in this region. Thus, if the low values of TL observed behind the Gaussian hills are within the level of numerical errors, the modeling results become more uncertain. As the receiver moves further away from the hill, the models converge toward a common TL. As discussed, the large apparent TL values in the shadow zone behind for instance a hill would be expected to be significantly decreased if turbulence is considered in the computations. The same deviation behind the hills is seen for the higher flow resistivities, r ¼ 250 kPa s/m 2 and r ¼ 2000 kPa s/m 2 . However, the models are not showing a similar convergence at longer distances in the case of a downward refracting atmosphere. Instead significant differences at some distances are obtained. As discussed for the undulating terrain and further shown in Fig. 11, this implies a phase shift between the reflections at the ground captured by the GTPE and by the remaining models. The difference increases for the highest ground flow resistivity of r ¼ 2000 kPa s/m 2 , and is further noticed to be the most prominent at the lowest hill height reaching TL differences up to 20 dB at 50 Hz. The differences are especially observed for downward refraction at both 50 and 500 Hz. Propagation in a no-wind atmosphere yields differences in the shadow zone behind the hills, and only smaller differences at longer distances. A similar observation is done for upward refracting propagation at 50 Hz. Yet/However for 500 Hz, the shadow zone caused by the refraction appears earlier, causing increasing deviations in the models due to numerical precision differences.
The significance of the destructive interferences seems to decrease with the steepness of the Gaussian hill. Thus, the interaction with the hill seemingly dominates the propagation through the rest of the domain. Since the GTPE shows a larger TL in the shadow zone behind the hill, the contribution from the interference patterns caused by the ground impedance and observed by the interference minima after the hill is smaller than for the BTPE and JEPE. Since the interference with the ground occurs earlier for a taller hill, the first minimum is likewise expected earlier. Thus, the minimum disappears in the more apparent shadow zone obtained with the GTPE. Eventually, for the GTPE at the steepest hill, there are no visible interferences at further distances while the BTPE and JEPE still capture interference at the bottom of the hill.
Since the BTPE and GTPE evaluate the square-root operator in the PE differently, the refraction will be treated differently in the models. Thus, it is likely that the varying treatment of the refraction causes differences in the obtained results. The significant interferences observed in this study are especially expected to be caused by the single point source representation of the wind turbine. When applied to wind turbine noise propagation, the use of distributed uncorrelated point sources along the rotor plane has previously shown good agreement with wind turbine noise measurements . The reason is that the distributed point sources average out the interference effects from the ground at the lower boundary. Thus, the distinct phase difference between the interferences observed by the BTPE and GTPE, respectively, would diffuse when introducing a multiple point source representation of the rotor.

Full spectrum TL
In addition to the range-dependent TL comparison, the numerical part of the inter-model comparison includes spectral comparisons at selected ground distances from the sources. The spectra are evaluated at the centre 1/3 octave band frequencies from 20 Hz to 2 kHz. In this analysis, two different Gaussian hills (the lowest, Hill 1, with H ¼ 20.4 m, and the steepest slope, Hill 4, with H ¼ 84.9 m) as well as the undulating terrain were selected. In this way, the effects of the phase shifted interferences observed in for instance Fig.  11 are quantified. The resulting spectra are shown as the difference between the BTPE and the GTPE at each frequency, DTLð f i Þ ¼ TL GTPE ð f i Þ À TL BTPE ð f i Þ. Furthermore, the ground impedance and atmospheric refraction are varied in a similar way as in the studies of range-dependent propagation. The resulting DTLð f i Þ for the three terrain cases are presented in Figs. 12-14. For all three terrain cases (Hill 1, Hill 4, and the undulating terrain) the largest differences throughout the spectra are seen for the highest impedance of r ¼ 2000 kPa s/m 2 at downward refraction. This agrees with the observations done for the frequencies of 50 and 500 Hz in the range-dependent study. In addition, the constant sound speed/no wind atmosphere does however introduce significant DTLð f i Þ at some receiver positions as well. The upward refraction is seen to cause large differences between the models for a wider range of frequencies. Thus, as the distance from the source increases, lower frequencies introduce a shadow zone. Comparing the spectra obtained for Hill 1 and 4, respectively, the large deviations between the PE models obtained at Hill 1 for downward refraction is observed for a wider range of frequencies as well. Hill 1 in general introduces the most significant spectral differences. It is, however, observed for both Hill 1 in Fig. 13 and Hill 4 in Fig. 14, that the position at the hill tops at 500 m ground distance yields a higher agreement between the two models at all frequencies. This is likely to be caused by the smaller observer angle at the hill tops and by the fact that any interferences caused by the ground effects are expected to occur further down the propagation path. In the undulating terrain in Fig. 12, alternations in the terrain are introduced earlier leading to earlier interferences. Thus, apparent differences between the BTPE and GTPE models are already observed at 500 m ground distance.
It is evident that the uncertainties of the PE models are highly dependent on the captured interferences. In general, the analysis of a few selected receiver points along a propagation path can cause the chosen grid points to be subject to the phase shifts observed in this study, while other neighbouring grid points might, for instance, capture an intersection of the models. Thus, the agreement of the spectra obtained by the PE models might vary significantly with the chosen receiver positions.

Effect of numerical step size
As the used step sizes differ between the GTPE and the BTPE, a convergence study of the BTPE was performed for the three atmospheric profiles described in Sec. III A and the four Gaussian hills described in Sec. III B by varying Dr ¼ ½0:1; 0:5; 1; 2; 5; 10k for 50 Hz and Dr ¼ ½2:5; 5; 10k for 500 Hz. An example is given in Fig. 15 for the lowest Gaussian hill and r ¼ 2000 kPa s/m 2 , which is later in the article shown to introduce the most significant interference patterns. The figure shows average convergence within 60.5 dB between Dr ¼ 0:1k and 2:5k for 50 Hz and between Dr ¼ 2:5k and 10k for 500 Hz. The interference patterns are also changing between the different discretizations with around 50 m delay for the 10k compared to the 0:1k at 50 Hz and around 10 m delay for 15k compared to 2k for 500 Hz. Thus, this supports the choice of a finer grid with a range step size of Dr ¼ 2:5k for computations done at 50 Hz, while Dr ¼ 10k is sufficient for 500 Hz computations.

IV. FIELD MEASUREMENTS
In this section, the validations against two different measurement campaigns, a loudspeaker and a wind turbine noise campaign, are considered. The measurement setups are briefly described in separate sections and are each followed by comparisons with the computational results. All spectral results are obtained by computations at 1/3 octave band centre frequencies.

A. Loudspeaker measurements
This section compares the BTPE and GTPE models with sound propagation measurements. The focus is initially kept on noise emitted from a loudspeaker. In this way, the source can be approximated as a point source and the source intensity is well known through microphone measurements 2 m away. The experimental method will briefly be described, and further details about the data selection and the background noise correction method can be found in Nyborg et al. (2022). A sketch of the experimental setup is shown in Fig. 16. The loudspeaker and the source microphone are mounted at the nacelle of a wind turbine with a hub height of 109 m. Seven microphones positioned in a line spanning from a distance of r ¼ 158 m to r ¼ 978 m monitored the noise as the loudspeaker was periodically turned on for 12 min and off for 12 min. The ground flow resistivity long the line of microphones was measured to be r ¼ 250 kPa s/m 2 up to a ground distance of r ¼ 270 m from the wind turbine and r ¼ 500 kPa s/m 2 for the rest of the domain. Furthermore, a meteorological mast near the wind turbine provided measurements of the wind speed, wind direction, temperature, and relative humidity. In total, 2 h of measurements are used for this comparison. The measurements of the meteorological conditions have been averaged for every 2 min, and are used to generate the speed of sound field in the computational domain. In addition, the noise measurements are evaluated in the 2 min equivalent sound, L eq2min ðf Þ, for 1/3 octave bands between 100 Hz and 2 kHz during the 2 h. This approach yields a statistical evaluation of the BTPE and GTPE models.
The sound speed profile consists of a linear temperature profile determined by the temperature gradient measured between 3 and 105 m height and a power law wind speed profile based on the measured wind speed and the shear coefficient at hub height. The wind speed is projected to the computational two-dimensional (2D) slice of the PE models spanning from the source to the receiver by the angle of the wind direction relative to the slice. This is done under the assumption that the flow is uniform since the terrain relative to the height of the loudspeaker position is approximately flat. The mean speed of sound profile, c(z), for the entire measurement period with the loudspeaker on is plotted in Fig. 17.

B. Modeling results
The results from the BTPE and GTPE models are compared with the experimentally obtained L eq2min ðf Þ in Fig. 18. In addition to the two higher fidelity models, the ISO standard model (ISO, 1997) is included. It is noted that the results from the ISO model given here, are interpolated from octave bands to 1/3 octave bands, yielding the method to not completely follow the standard use of the model. The model is therefore noted as the ISO-1/3 octave model. The measured L eq2min ðf Þ is corrected for background noise by using the newly implemented deconvolution method  and the IEC method (IEC, 2012), respectively. The signal obtained by using either method is shown in Fig.  18. Thus, for some frequencies the deconvolution method is not following the IEC procedure by yielding a background noise-corrected signal at signal-to-noise ratios lower than 3 dB.
Comparing the modeling results with the measurements, both models agree with the measured 2 min equivalent measured sound, L eq2min ðf Þ. However, the PE models are generally performing better than the more simple ISO-1/3 octave model. For all models and at all microphones, the results are starting to deviate from the measurements going beyond 800 Hz. This was discussed by Nyborg et al. (2022) to be caused by the directivity of the loudspeaker resulting in spectral dissimilarities between microphone 0 and the remaining microphones at higher frequencies. As in the work of Nyborg et al. (2022), this discrepancy is further investigated by comparing the measured/computed difference between two microphones. Consequently, the influence from microphone 0 is eliminated, and only TL effects evaluated by the models and the atmospheric absorption are compared. In Fig. 19, DL p denotes the TL difference. When evaluating the TL difference like this, the higher frequencies are better captured by the models, which supports that the directivity of the source spectrum causes the deviations seen in Fig. 18. In Fig. 19, it becomes more apparent that the ISO-1/3 octave model does not evaluate the variation with frequency to the same extent as the higher fidelity models do. Overall, the BTPE and GTPE agree reasonably with the measured TL difference. However, some significant differences are seen for the BTPE at 500 Hz in the top panel and bottom left panel and at 1 kHz in the bottom right panel of Fig. 19. As the same differences are seen in Fig. 18 for microphones 4 and 5, they are most likely caused by destructive interferences evaluated by the BTPE model.

C. Ytterberg measurements
This comparison case considers computations of propagation from a wind farm further described by Conrady (2019). Inputs in the form of altitude and sound speed profile have been calculated as described by Bolin et al. (2020). The effective speed of sound is given in the locations of the wind turbine and the microphone and has been interpolated to each grid point between the two positions. Two measurement occasions have been selected providing the 10 min equivalent measured sound, L eq10min ð f Þ, for 1/3 octave bands between 25 Hz and 2.5 kHz. The cases are referred to as cases 1 and 2. At both measurement times, the ground was covered with snow according to the weather forecast (see  for further details) and therefore the ground flow resistivity has been set to r ¼ 10 kPa s/m 2 in accordance with Embleton et al. (1983) for dry snow. Furthermore, the terrain takes the form of a hill sloping downwards in the direction from the wind turbine to the microphone with a total ground distance of r ¼ 955.5 m. Due to the low ground flow resistivity, no significant reflections at the ground are expected. Therefore, the wind turbine is modeled with a single point source at hub height, z ¼ 105 m.

D. Modeling results
The modeling results of the two Ytterberg wind turbine cases are obtained for the BTPE and GTPE model. The sound power level of the source is given by IEC measurements (IEC, 2012) at a time different from the two cases presented here. The computed and modeled L eq10min spectra are presented in Fig. 20 with the associated measured speed of sound profiles at the source and receiver positions. It should be noted, that for these computations different impedance models are used in the BTPE and GTPE, respectively. While the BTPE uses the impedance model by Delany and Bazley (1970), the GTPE continues to use the impedance model (Attenborough, 1985). The use of different impedance models is especially expected to cause differences at lower frequencies. This is seen to be the case for both cases 1 and 2 when comparing the two models, where the use of the Delany and Bazley model in the BTPE results in a better agreement with the measurements for low frequencies, especially in case 2. In both cases, the BTPE model seems to evaluate larger interferences at high frequencies than the GTPE. In particular, in case 2 for frequencies of f ! 800 Hz, these interferences are observed, while the measurements and the GTPE model provide smoother spectra. In case 1, the GTPE shows interference peaks for some frequencies as well, but still shows an overall smoother spectrum than the BTPE. Both models yield a generally better comparison to the measurements in case 2. This could, among other things, be caused by the applied sound power level being more representative in this case or a potentially higher uncertainty in the background noise in case 1.

V. CONCLUSIONS
This paper presented an inter-model comparison study of two sound propagation models: the BTPE model and the GTPE model. The study includes theoretical cases of varying complex terrains, ground impedances, and speed of sound profiles as well as validations against two measurement campaigns. In addition, the results using the JEPE model were presented for the numerical comparison cases. The aim of the paper was thus to quantify any differences between PE models in the context of acoustic propagation from wind turbines. Through the comparisons, it is evident that propagation in complex terrains with acoustically harder surfaces yields larger differences, especially for downward refraction. This was seen to be caused by the significant phase shifts of the destructive interference patterns which can be a result of the different grid spacings used in the models. The large interferences are, however, diffused when using multiple point sources distributed along in this example a wind turbine rotor. Thus, this approach reduces the difference between the PE methods and between numerical and experimental results significantly for wind turbine noise purposes. In an upward refracting atmosphere, the interferences caused by the ground surface are not as prominent. The models, therefore, agree reasonably for shorter distances and lower frequencies but start to deviate when entering the shadow zone at larger distances and higher frequencies. The shadow zone is already an area of high numerical error due to the unrealistic low sound pressure levels obtained when neglecting the turbulence and thus the scattering of sound into the shadow zone. However, when considering wind turbine noise, the energy of the sound source is expected to be significantly higher at lower frequencies, while the atmospheric absorption causes the highfrequency noise to decay faster. With this in mind, the inclusion of turbulence at the higher frequencies might be of less importance in these cases.
The models agree reasonably with the loudspeaker and wind turbine measurements for a majority of the 1/3 octave band centre frequencies. Some deviations do however occur and are possibly caused by the uncertainty of the used source strength, the uncertainty of the background noise correction, or the modeled interferences at the specific receiver position. These deviations should be considered numerical errors when applying the Parabolic Equation models to the study of wind turbine noise propagation.

ACKNOWLEDGMENTS
This work considers loudspeaker measurements funded by the Innovation Fund Denmark in the project entitled "DecoWind-Development of low-noise and cost-effective wind farm control technology" under project number 8055-00041B. Finally, the authors thank the reviewers and the associate editor of the manuscript for the constructive feedback and immense help during the editing process.