© The Author(s), Published by EDP Sciences, 2024

1 Introduction

Auralization of elevated environmental noise sources, such as aircraft, wind turbines, and drones, is a powerful tool to support the development of low-noise technologies and to conduct studies on human reactions to noise. Long-term exposure to aircraft noise can have severe consequences for health of residents near airports [1]. Therefore, the development of low-noise aircraft and flight operations is an important aspect of a sustainable aviation future.

Auralization allows simulation data to be made audible. For the application of auralizations in the development of low-noise technologies, auralizations must be highly realistic. Realistic auralizations of elevated noise sources, such as aircraft, typically consider multiple sound propagation effects, including geometrical spreading, Doppler effect, and air absorption. Auralizations can be pre-computed or generated dynamically in a real-time framework as part of an interactive virtual reality environment [2, 3]. Additionally, we argue that turbulence in the atmospheric boundary layer leads to well-audible effects, and must be considered to achieve realistic auralizations of outdoor sound propagation.

Atmospheric turbulence consists of fluctuations in wind velocity and temperature. Fluctuating wind velocity and temperature lead to fluctuations in the effective sound speed

$$c_{\mathrm{eff}}=c\left(T\right)+u,$$(1)

where c(T) is the temperature dependent speed of sound and u is the wind velocity component in the direction of sound propagation. In consequence, sound propagated through atmospheric turbulence exhibits fluctuations in the acoustical phase and amplitude.

Figure 1 presents two Airbus A320 aircraft flyovers during take-off in different meteorological conditions. The spectrogram at the top displays a flyover measured in low-turbulence conditions at 7 a.m. The spectrogram at the bottom displays a flyover that was measured during high-turbulence conditions at noon. In contrast to the first spectrogram, the second one shows noticeable vertical patterns that indicate strong amplitude fluctuations caused by atmospheric turbulence.

Figure 1 Spectrograms of measured Airbus A320 flyover during take-off at approximately 650 m above ground. Colors indicate power spectral density (PSD). Top: Low turbulence case at 7 a.m. Bottom: High turbulence case at 12 a.m.

In the past decade, modeling the effect of atmospheric turbulence in auralization has gained more attention, as it considerably impacts the audible characteristics of the acoustical phase and amplitude [4–8]. Simultaneously, the theoretical description of sound propagation in the atmospheric boundary layer has seen large improvements and was extended to describe vertical and slanted sound propagation [9].

While the model presented in [9] and validated in [10] is the best model currently available to theoretically describe turbulence induced fluctuations, it is only valid in the weak-scattering regime. The weak-scattering regime implies relatively short propagation distances and weak turbulence. Both conditions are usually violated for daytime aircraft flyover. We will show in Section 5 that the application of the theoretical model for the weak-scattering regime to aircraft flyover, often outside of the range of applicability of the model, leads to overestimated amplitude fluctuations in most meteorological daytime conditions. We will then leverage the empirical data presented in Section 3 to extend the applicability of the model to aircraft flyovers by proposing a series of modifications and establishing a semi-empirical model.

This paper is structured as follows. Section 2 briefly explains the origin of log-amplitude fluctuations and its property to saturate for sound propagation in strong-scattering regime. Existing theoretical approaches to model log-amplitude fluctuations and saturation in auralization are reviewed. Section 3 presents the collection and the processing of the empirical data. The empirical dataset includes over 5000 aircraft flyover events during take-off and landing. Positional data of the aircraft allows detailed information on the flight trajectory. In addition, meteorological data is used to characterize the state of the atmospheric boundary layer (ABL). Section 4 presents the analysis of the acoustical field measurements with regard to its dependency on the state of the ABL and signs of amplitude saturation. Section 5 presents the semi-empirical model developed in this work. The semi-empirical model considers an effective propagation length based on the boundary layer height, saturation of log-amplitude fluctuations, turbulence decay time during periods of decreasing solar radiation before sunset, and turbulence production by nocturnal low-level-jets. Section 6 discusses strengths and limitations of the proposed semi-empirical model. Section 7 demonstrates an example aircraft flyover auralization with the semi-empirical model. Section 8 concludes the article.

2 Theoretical background

2.1 Conversion of log-amplitude fluctuations to dB

In this article, theoretical formulations for variance of log-amplitude fluctuations 〈χ²〉 will be compared to the standard deviation of measured fluctuations in dB, referred to as σ_χ. To allow the comparison between theoretical log-amplitude variance and measured fluctuations, 〈χ²〉 can be rewritten as variance in dB as

$${ \sigma }_{\chi }^2=\langle {\left(20\lg \left(A/A_0\right)\right)}^2\rangle =\langle {\chi }^2\rangle \left(20\lg \right(e){)}^2,\hspace{1em}[\mathrm{dB}]$$(2)

where the angle brackets 〈〉 denote an ensemble average, A is the sound pressure amplitude with turbulence and A₀ the sound pressure amplitude in a non-turbulent atmosphere. The standard deviation of amplitude fluctuations in dB is then

$${\sigma }_{\chi }=\langle \chi \rangle \left(20\lg \left(e\right)\right)=\langle \chi \rangle 8.7\hspace{1em}\left[\mathrm{dB}\right].$$(3)

2.2 Origin and saturation of amplitude fluctuations

An introduction to the processes that lead to log-amplitude and phase fluctuations of sound propagated in a inhomogeneous medium is given in [11] and briefly summarized in the following.

Long distance sound propagation in a homogeneous medium can be described by a planar wavefront perpendicular to the sound propagation direction. Sound speed fluctuations along the wavefront diffract the sound and cause ripples in the formerly planar wavefront. For weak scattering, the wavefront curvature leads to weak focusing and defocusing of the sound, resulting into weak log-amplitude fluctuations. Log-amplitude fluctuations are most sensitive to inhomogeneities in the order of the Fresnel zone, while phase fluctuations are mainly caused by the largest scales of inhomogeneities in the atmosphere.

Strong scattering results into folding of the wave front and leads to closely spaced, correlated ray paths, which are called microrays. Interference of these microrays leads to strong focusing effects of the sound. The same effect can be observed for light causing strong focusing areas on the bottom of a swimming pool.

For very long propagation distances and strong diffraction, variance of log-amplitude fluctuations of monochromatic signals saturate at 〈χ²〉 = π²/24 corresponding to a standard deviation (STD) of sound level fluctuations of approximately 5.6 dB. Microrays become increasingly uncorrelated and cause less interference. The scintillation index (SI) describes the variability of the signal as

$$\mathrm{SI}=\frac{\langle I^2\rangle -\langle I{\rangle }^2}{\langle I{\rangle }^2},$$(4)

with I being the sound intensity.

The SI approaches 1 in full saturation. Figure 2 illustrates the behaviour of the SI with increasing range and scattering. In the partially-saturated regime, the SI is greater than 1. For more details, the reader is referred to [11] and [12].

Figure 2 Qualitative dependence of scintillation index with range. The range at which full saturation occurs is frequency dependent. Red area corresponds to behavior that is expected for aircraft flyover. Adapted from Figure 4.9 in [11].

Aircraft overflights just before landing or after takeoff can often be classified as partially saturated, as shown in Section 4.3.

2.3 Existing approaches to model log-amplitude fluctuations in auralizations

Two theoretical models have been used in the past to generate amplitude fluctuations in auralizations. A first approach was based on a Gaussian turbulence model to describe correlation functions of amplitude and phase fluctuations by [13]. The model was used in [5] to generate aircraft auralizations with turbulence induced amplitude and phase fluctuations by generating sequences of acoustic scintillations. Recently, [8] advanced the auralization approach introduced in [5] and integrated a more sophisticated model to describe amplitude and phase fluctuations published by [9]. In contrast to the preceding model, it is based on the von Kármán turbulence spectrum and considers vertical sound propagation. [8] elaborate on the differences of the two models, point out advantages and limitations, and compare results for horizontal propagation.

[8] restricted the variance of log-amplitude fluctuations to 〈χ²〉 < 0.8 for auralization purposes based on measurements presented in [14]. There, amplitude fluctuations of broadband signals emitted by a loudspeaker mounted at a wind turbine nacelle have been studied. Microphones were placed close to the ground in horizontal distances up to 1278 m. The variance of log-amplitude fluctuations 〈χ²〉 shows signs of saturation but reach values around 1 corresponding to level fluctuations of more than 8 dB.

Instead of determining a maximum level of log-amplitude fluctuations, [15] defined a saturation distance as the distance at which the energy of the coherently propagating sound wave has decayed by e⁻¹ compared to sound propagation in a homogeneous field. Wenzel proposed a model to estimate the order-of-magnitude of the saturation distance which shows rough agreements with measurements. The model was applied in [5] to approximate the saturation effect for log-amplitude fluctuations in aircraft auralization. The saturation distance depends on the correlation length L_cor, which is not straightforward to determine for slanted or vertical sound propagation.

2.4 Theoretical model for amplitude fluctuations in weak-scattering regime

In the following, the theoretical model for variance of log-amplitude fluctuations in the weak-scattering regime described in [9] and [10] is introduced. Log-amplitude fluctuations due to atmospheric turbulence is defined as

$$\chi =\ln \left(A/A_0\right).$$(5)

Sound pressure p is represented in the frequency domain as its Fourier transform $\widehat{p}$ and expressed using the Rytov method as

$$\widehat{p}={\widehat{p}}_0\exp (\chi +\mathrm{i\varphi }),$$(6)

where ϕ is the phase fluctuation due to atmospheric turbulence, and ${\widehat{p}}_0$ is the transform in a non-turbulent atmosphere.

For vertical or slanted sound propagation in the direction of the coordinate l with length L = h/cos θ, the integration variable η = l/L along the path is introduced. η = 0 corresponds to the starting point of the sound propagation path, η = 1 to its end, the height is given by z = ηh. Log-amplitude variance for a source at z = 0 and a receiver at height z = h is then given by

$$\langle {\chi }^2\rangle =\frac{{\pi }^2{k}^{2}h}{2\cos \theta }{\int }_0^1\mathrm{} \mathrm{d}\eta {\int }_0^{\infty }\mathrm{} {\mathrm{\Phi }}_{\mathrm{eff}}(\eta h;\kappa )\left[1-\cos \left(\frac{\eta (1-\eta ){\kappa }^{2}L}{k}\right)\right]\kappa \mathrm{d}\kappa .$$(7)

Here, k is the acoustical wavenumber, h is the height of the receiver, κ is the turbulence wavenumber, and Φ_eff is the effective turbulence spectrum. Equation (7) assumes line-of-sight propagation.

The effective turbulence spectrum is given by:

$${\mathrm{\Phi }}_{\mathrm{eff}}(z;\kappa )=\frac{{\tilde{\mathrm{\Phi }}}_T(z;\kappa )}{T_0^2}+\frac{E(z;\kappa )}{\pi c_0^2{\kappa }^2},$$(8)

with the generalized von Kármán spectra for temperature and wind velocity fluctuations:

$${\tilde{\mathrm{\Phi }}}_T(z;\kappa )=\frac{\mathrm{\Gamma }\left(\alpha \right)}{{\pi }^{3/2}\mathrm{\Gamma }(\alpha -3/2)}\frac{{\sigma }_T^2\left(z\right)L_T^3\left(z\right)}{{\left(1+{\kappa }^2{L}_T^2\left(z\right)\right)}^{\alpha }},$$(9)

$$E(z;\kappa )=\frac{4\alpha \mathrm{\Gamma }\left(\alpha \right)}{\sqrt{\pi }\mathrm{\Gamma }(\alpha -3/2)}\frac{{\sigma }_v^2\left(z\right)L_v^5\left(z\right){\kappa }^4}{{\left(1+{\kappa }^2{L}_v^2\left(z\right)\right)}^{\alpha +1}}.$$(10)

Here, Γ(n) = (n − 1)! is the gamma function, T₀ is a reference temperature, c₀ the reference speed of sound. The ordinary von Kármán spectra is obtained by setting α = 11/6. ${\sigma }_T^2$, ${\sigma }_{\mathrm{vs}}^2$, and ${\sigma }_{\mathrm{vb}}^2$ are the variances of temperature fluctuations and of shear-produced and buoyancy-produced wind velocity fluctuations. The corresponding length scales of the fluctuations are denoted as L_T, L_vs, and L_vb.

The variances and length scales of the respective fluctuations in the atmospheric boundary layer are expressed based on Monin-Obukhov (MOST) and mixed-layer similarity theories. The effective turbulence spectrum can then be expressed in terms of height above ground z, boundary layer height z_i, and the turbulence scaling parameters, i.e. friction velocity $u_{*}$, and surface sensible heat flux Q_H. The variances are given by:

$${\sigma }_T^2\left(z\right)=\frac{4.0T_{\mathrm{*}}^2}{{\left(1-10z/L_o\right)}^{2/3}},{\sigma }_{\mathrm{vs}}^2=3.0u_{\mathrm{*}}^2,{\sigma }_{\mathrm{vb}}^2=0.35w_{\mathrm{*}}^2.$$(11)

The length scales for temperature and velocity fluctuations are given by:

$$L_T\left(z\right)=2.0z\frac{\left(1-7z/L_o\right)}{\left(1-10z/L_o\right)},L_{\mathrm{vs}}\left(z\right)=1.8z,L_{\mathrm{vb}}=0.23z_i.$$(12)

Here, $T_{*}=-QH/{\rho }_0{c}_p{u}_{*}$ is the surface layer temperature scale, $w_{\mathrm{*}}={\left(z_{i}g{Q}_H/{\rho }_0{c}_p{T}_{2\mathrm{m}}\right)}^{1/3}$ the convection velocity scale, and $L_o=-u_{\mathrm{*}}^3{T}_{\mathrm{surf}}{\rho }_0{c}_p/\left(g{\kappa }_v{Q}_H\right)$ is the Obukhov length. The parameters ρ₀, c_p, g, and κ_v = 0.4 describe the air density, specific heat, gravitational acceleration, and the von Kármán constant. T_surf is the earth surface temperature, T_2m is the mean temperature at 2 m height.

For sound propagation from an elevated source at height h_s to a receiver at height h_r, the argument ηh is replaced by (1 − η)h_s + ηh_r.

Equation (7) is used in [8] and [16] for auralization of wind farm noise.

2.5 Theoretical log-amplitude fluctuations in terms of strength and wave parameters

An alternative formulation of equation (7) is presented in [17] for horizontal sound propagation. This alternative formulation for log-amplitude fluctuations allows to determine a theoretical boundary between the weak and the strong-scattering regime. Based on the formulation for horizontal sound propagation given in [17], we derived the equivalent formulation for vertical and slanted sound propagation. The derivation is presented in the following.

The integrals over the turbulence spectra can be expressed in terms of the integral length scales ${\mathcal{L}}_T$ and ${\mathcal{L}}_v$ [17, 18]:

$$\frac{{\sigma }_T^2{\mathcal{L}}_T}{2{\pi }^2{T}_0^2}={\int }_0^{\infty }\mathrm{}{\tilde{\mathrm{\Phi }}}_T\left(\kappa \right)\kappa \mathrm{d}\kappa ,$$(13)

$$\frac{2{\sigma }_v^2{\mathcal{L}}_v}{\pi }={\int }_0^{\infty }\mathrm{}\frac{E\left(\kappa \right)}{\kappa }\mathrm{d}\kappa .$$(14)

The integral length scales for the generalized von Kármán spectra are given by

$${\mathcal{L}}_T=\frac{\sqrt{\pi }\mathrm{\Gamma }(\alpha -1)}{\mathrm{\Gamma }(\alpha -3/2)}{L}_T, {\mathcal{L}}_v=\frac{\sqrt{\pi }\mathrm{\Gamma }(\alpha -1)}{\mathrm{\Gamma }(\alpha -3/2)}{L}_v.$$(15)

Inserting equation (8) into equation (7) and multiplying out gives

$$\begin{array}{ll}\langle {\chi }^2\rangle =\frac{{\pi }^2{k}^{2}h}{2\cos \theta }& \left[{\int }_0^1\mathrm{}{\int }_0^{\infty }\mathrm{}\frac{{\tilde{\mathrm{\Phi }}}_T(\eta h,\kappa )}{T_0^2}\kappa \mathrm{d}\kappa \mathrm{d}\eta -{\int }_0^1\mathrm{}{\int }_0^{\infty }\mathrm{}\frac{{\tilde{\mathrm{\Phi }}}_T(\eta h,\kappa )}{T_0^2}\cos \left(\frac{\eta (1-\eta ){\kappa }^{2}L}{k}\right)\kappa \mathrm{d}\kappa \mathrm{d}\eta \right]+\\ & +\frac{{\pi }^2{k}^{2}h}{2\cos \theta }\left[{\int }_0^1\mathrm{}{\int }_0^{\infty }\mathrm{}\frac{E(\eta h,\kappa )}{\pi c_0^2\kappa }\mathrm{d}\kappa \mathrm{d}\eta -{\int }_0^1\mathrm{}{\int }_0^{\infty }\mathrm{}\frac{E(\eta h,\kappa )}{\pi c_0^2\kappa }\cos \left(\frac{\eta (1-\eta ){\kappa }^{2}L}{k}\right)\mathrm{d}\kappa \mathrm{d}\eta \right].\end{array}$$(16)

Splitting up the function, we can write

$$\langle {\chi }^2\rangle =\langle {\chi }_T^2\rangle +\langle {\chi }_v^2\rangle .$$(17)

For brevity and better readability, the equations are given separately for temperature and velocity fluctuations and the dependencies on z = ηh are omitted in the following. Replacing the first integral in equation (16) by equation (13), excluding ${\sigma }_T^2{\mathcal{L}}_T/2{\pi }^2{T}_0^2$ and inserting equation (15) gives:

$$\langle {\chi }_T^2\rangle =\frac{\sqrt{\pi }\mathrm{\Gamma }(\alpha -1)}{\mathrm{\Gamma }\left(\alpha -\frac{3}{2}\right)}\frac{k^{2}h}{4T_0^2\cos \theta }{\int }_0^1\mathrm{}{\sigma }_T^2{L}_T{\int }_0^{\infty }\mathrm{}\left[1-\frac{2(\alpha -1)}{{\left(1+{\kappa }^2{L}_T^2\right)}^{\alpha }}{L}_T^2\kappa \cos \left(\frac{\eta (1-\eta ){\kappa }^{2}L}{k}\right)\right]\mathrm{d}\kappa \mathrm{d}\eta .$$(18)

In the square brackets, the Gamma function cancels out because of the relationship Γ(α)/Γ(α − 1) = (α − 1). Substituting the integration variable by $\kappa =\frac{\widehat{\kappa }}{L_T}$ gives

$$\langle {\chi }_T^2\rangle =\frac{\sqrt{\pi }\mathrm{\Gamma }(\alpha -1)}{4\mathrm{\Gamma }\left(\alpha -\frac{3}{2}\right)}\frac{k^{2}h}{T_0^2\cos \theta }{\int }_0^1\mathrm{}{\sigma }_T^2{L}_T{\int }_0^{\mathrm{\infty }}\mathrm{}\left[1-\frac{2(\alpha -1)}{{\left(1+{\widehat{\kappa }}^2\right)}^{\alpha }}\widehat{\kappa }\cos \left(\eta (1-\eta ){\widehat{\kappa }}^2{D}_T\right)\right]\mathrm{d}\widehat{\kappa }\mathrm{d}\eta ,$$(19)

with the wave parameter D_T:

$$D_T\left(z\right)=\frac{L}{k{L}_T^2\left(z\right)}.$$(20)

Substituting the integration variable again by $\tilde{\kappa }={\widehat{\kappa }}^2$ gives:

$$\langle {\chi }_T^2\rangle =\frac{\sqrt{\pi }\mathrm{\Gamma }(\alpha -1)}{4\mathrm{\Gamma }\left(\alpha -\frac{3}{2}\right)}\frac{k^{2}h}{T_0^2\cos \theta }{\int }_0^1\mathrm{}{\sigma }_T^2{L}_T\left[1-H_T\left(D_T,\eta \right)\right]\mathrm{d}\eta ,$$(21)

with

$$H_T\left(D_T,\eta \right)={\int }_0^{\infty }\mathrm{}\frac{(\alpha -1)}{{\left(1+\tilde{\kappa }\right)}^{\alpha }}\cos \left(\eta (1-\eta ){\tilde{\kappa }}^2{D}_T\right)\mathrm{d}\tilde{\kappa }.$$(22)

The function H_T(D_T, η) describes the effect of diffraction on 〈χ²〉 [17]. The term

$${\mathrm{\varphi }}_T^2=\frac{\sqrt{\pi }\mathrm{\Gamma }(\alpha -1)}{2\mathrm{\Gamma }\left(\alpha -\frac{3}{2}\right)}\frac{k^{2}h}{T_0^2\cos \theta }{\sigma }_T^2{L}_T$$(23)

is called the strength parameter ${\mathrm{\varphi }}_T^2$ and corresponds to the variance of phase fluctuations in the geometrical acoustics approximation due to temperature fluctuations. It should not be confused with the turbulence spectrum of temperature fluctuations ${\tilde{\mathrm{\Phi }}}_T$ in equation (13). We decided to keep the convention of giving the turbulence spectra Φ and the strength parameter ${\mathrm{\varphi }}^{\mathbf{2}}$ the same Greek letter in order to be consistent with the literature used in this article. However, we use different fonts to make the distinction easier.

Similarly, the second part of equation (16) can be transformed into

$$\langle {\chi }_v^2\rangle =\frac{\sqrt{\pi }\mathrm{\Gamma }(\alpha -1)}{\mathrm{\Gamma }\left(\alpha -\frac{3}{2}\right)}\frac{k^{2}h}{c_0^2\cos \theta }{\int }_0^1\mathrm{}{\sigma }_v^2{L}_v\left[1-H_v(D_v,\eta )\right]\mathrm{d}\eta ,$$(24)

with

$$H_v(D_v,\eta )={\int }_0^{\mathrm{\infty }}\mathrm{}\frac{\alpha (\alpha -1)}{{\left(1+\tilde{\kappa }\right)}^{\alpha +1}}\tilde{\kappa }\cos \left(\eta (1-\eta ){\tilde{\kappa }}^2{D}_v\right)\mathrm{d}\tilde{\kappa },$$(25)

and the the wave parameter D_v:

$$D_v\left(z\right)=\frac{L}{k{L}_v^2\left(z\right)}.$$(26)

The strength parameter ${\mathrm{\varphi }}_v^2$ for wind velocity fluctuations is given by

$${\mathrm{\varphi }}_v^2=\frac{2\sqrt{\pi }\mathrm{\Gamma }(\alpha -1)}{\mathrm{\Gamma }\left(\alpha -\frac{3}{2}\right)}\frac{k^{2}h}{c_0^2\cos \theta }{\sigma }_v^2{L}_v.$$(27)

For considering log-amplitude variance caused by fluctuations in temperature, shear-produced and buoyancy produced wind velocity fluctuations, equation (16) can be expressed in terms of strength parameters ${\mathrm{\varphi }}_{T,\mathrm{vs},\mathrm{vb}}^2$, wave parameters D_T,vs,vb and the functions H_T,v:

$$\langle {\chi }^2\rangle ={\int }_0^1\mathrm{}\frac{{\mathrm{\varphi }}_T^2}{2}\left[1-H_T(D_T,\eta )\right]\mathrm{d}\eta +{\int }_0^1\mathrm{}\frac{{\mathrm{\varphi }}_{\mathrm{vs}}^2}{2}\left[1-H_v(D_{\mathrm{vs}},\eta )\right]\mathrm{d}\eta +{\int }_0^1\mathrm{}\frac{{\mathrm{\varphi }}_{\mathrm{vb}}^2}{2}\left[1-H_v(D_{\mathrm{vb}},\eta )\right]\mathrm{d}\eta .$$(28)

A similar formulation of equation (28) is given in [9] using the extinction coefficient γ. As for equation (7), sound propagation from an elevated source at height h_s to a receiver at height h_r, can be considered by replacing the argument ηh by (1 − η)h_s + ηh_r.

This representation of log-amplitude fluctuations is the basis for the new semi-empirical model presented in Section 5, as it allows to estimate the boundary between weak and strong scattering (see Sect. 2.6).

2.6 Estimation of boundary between weak and strong scattering

Equations (7), (21) and (24) are valid in the weak-scattering regime and do not consider saturation effects. In the Rytov method, 〈χ²〉 must remain less than 1. However, this value is easily exceeded using equations (7) and (28) for relatively large propagation distances in the atmospheric boundary layer.

The representation of 〈χ²〉 based on wave and strength parameters D and ϕ can be used to determine the boundary between weak and strong-scattering regime. Equation (28) is only valid in the weak-scattering regime, if the coefficient b² is not exceeded:

$${\langle {\chi }^2\rangle |}_{\begin{array}{c}\mathrm{weak}\\ \mathrm{scat}.\end{array}}<b^2.$$(29)

Different values are assumed for the coefficient b² in literature. In [17] it was set to b² = 0.5. In [10], field measurements with b² > 0.2 showed signs of strong scattering and were excluded from the model validation.

Using the inequality by equating equations (28) and (29), the maximum strength parameter ${\mathrm{\varphi }}_{\max }^2$ in the weak-scattering regime is given by [17]:

$${\mathrm{\varphi }}_{\max }^2=\frac{2b^2}{{\beta }_T\left[1-H_T\left(D_T\right)\right]+{\beta }_{\mathrm{vs}}\left[1-H_v\left(D_{\mathrm{vs}}\right)\right]+{\beta }_{\mathrm{vb}}\left[1-H_v\left(D_{\mathrm{vb}}\right)\right]}.$$(30)

The coefficients β_T, β_vs, and β_vb are determined by the ratio of the strength parameters relative to the strongest strength parameter:

$${\beta }_{T,v}=\frac{{\mathrm{\varphi }}_{T,v}^2}{\max \left({\mathrm{\varphi }}_T^2,{\mathrm{\varphi }}_v^2\right)}.$$(31)

The wave and the strength parameter can be used to represent the wave propagation regime and boundaries between weak and strong scattering in a D–ϕ diagram. Figure 3 shows a qualitative sketch of a D–ϕ diagram. For low wave parameters, log-amplitude fluctuations approach the full-saturation regime via the partially-saturated regime, which is characterized by interference of microrays and strong focusing effects.

Figure 3 Qualitative D–ϕ diagram showing unsaturated, partially saturated and fully saturated regime depending on strength parameter ϕ and wave parameter D. Adapted from Figure 4.10 in [11] and Figure 5 in [17].

2.7 Further limitations of the weak-scattering model

While the limitation to the weak-scattering regime is the strongest limitation to equation (7), further limitations must be considered depending on the application.

First, the model assumes stationary ABL with well developed, continuous turbulence [19]. The ABL changes considerably before noon during the transition from nocturnal, stable boundary layers (SBL) to daytime convective boundary layers and again in the late afternoon when the ABL collapses and transitions back into a nocturnal boundary layer [20].

Second, the model is not valid for sound propagation at heights z > z_i, for stable boundary layers with Q_H < 0, and for intermittent turbulence, as the model is based on MOST and mixed-layer similarity theory. The height dependent expressions for ${\sigma }_T^2$, L_T and L_vs are based on MOST and hence only valid in the atmospheric surface layer. The surface layer is approximately the lowest 10% of the ABL [20]. Therefore, these parameters must be kept constant for z > 0.1 z_i.

Third, the model assumes isotropic turbulence. In reality, the ground blocks the vertical movement of air and leads to anisotropic behavior of the largest eddies within the ABL. An extension of the model to anisotropic turbulence is presented in [21]. The extended model was compared with measurements presented in [10] and showed clear improvements in the prediction of phase fluctuations. However, the isotropic assumption is justified for amplitude fluctuations, as these are mainly impacted by small-scale turbulence. Therefore, the isotropic model is used in this publication.

Despite these limitations, it was shown in [10, 21, 22], that measured phase and log-amplitude variances, spatial coherence functions and the probability density functions of the signal agree relatively well with the theoretical predictions, provided that the range of applicability is met, i.e. relatively short sound propagation within the ABL, convective and neutral ABL, and weak-scattering regime. In case of aircraft flyover auralization, at least one of these conditions is not fulfilled in most cases. Therefore, the purpose of this article is to present an empirical extension of the model, which can be applied to aircraft flyover auralization.

3 Field measurements

3.1 Acoustical data

3.1.1 Data collection

The acoustical data used in this publication is part of the measurement campaigns within the sonAIR projects conducted in 2014 [23] and 2022. Civil aircraft were measured during take-off and landing at Zurich Airport (ZRH). Data was collected from March to August throughout the day from 6 a.m. to 11 p.m. For the following data analysis, flyovers of aircraft of the type Airbus A319, A320, A321, A320neo, A321neo are included. The microphones were located within the approach and take-off paths at varying distances up to 20 km to the runway. The microphones were mounted on poles of 6 m or 10 m height depending on the measurement location. In total, the data considered in this publication was collected at twelve different measurement locations. The topography of the measurement locations varies from mostly flat areas to more hilly landscapes. Figure 4 shows one exemplary measurement location.

Figure 4 Mounted ultrasonic anemometer (front) and microphone (back).

The dataset encompasses 5305 flyover measurements, of which 1430 are during landing and 3875 during take-off. Figure 5 demonstrates the distances and angles of the aircraft at the moment of minimal distance between microphone and aircraft during the flyover. Blue (darker) points represent landings, red (lighter) points represent take-offs. Flyovers with distances L_min < 400 m have been discarded as the time level history of the measurements were dominated by sudden level changes due to the source directivity of the aircraft. In these cases, the turbulence-induced amplitude fluctuations could not be well separated from the effect of source directivity.

Figure 5 Distances and polar angles θ of measured aircraft flyovers for the moment of shortest distance between aircraft and microphone (θ = 0 is vertical). Blue (darker) points represent landings, Red (lighter) points represent take-offs.

Flight trajectories are available based on radar data from Zurich Airport. Since radar data is sometimes subject to uncertainties, radar data require prior processing steps in order to obtain trajectories with high accuracy suitable for de-Dopplerization. The process involves smoothing of abrupt fluctuations, corrections or rejection of incorrect trajectories, ensuring a more consistent and continuous trajectory, and adding additional information required for the data analysis, such as aircraft type. These processing steps are performed with tools integrated in the aircraft noise calculation program FLULA2, which has been developed at Empa [24].

3.1.2 Audio data post-processing

The acoustical measurements of aircraft flyovers are analyzed to yield a single number parameter characterizing the amplitude fluctuations. Figure 6 shows the processing steps of the acoustical measurements.

Figure 6 Flowchart of audio signal processing to estimate the amplitude fluctuation indicator σχ (dB) from measured aircraft flyovers

The measurements are de-Dopplerized to eliminate analysis artefacts caused by narrowband source components such as engine fan tones and wing cavity noise. The Doppler effect causes a frequency shift in these components. Frequency components entering or leaving the analyzed frequency band during the considered time window can distort the amplitude fluctuation.

After de-Dopplerization, the signals are bandpass filtered to exclude frequencies below 446 Hz and above 1120 Hz, corresponding to the lower and upper limits of the 500 Hz and 1 kHz one-third octave band respectively. The lower limit is chosen to eliminate low-frequency background noise and fluctuations in the ground effect. The higher limit is chosen to exclude birdsong. The large bandwidth reflects the broadband character of amplitude fluctuations and decreases the impact of background noise on the overall sound pressure level.

The sound pressure level time history L_p(t) of the bandpassed signal is computed as

$$L_p\left(t\right)=10 \mathrm{lo}{\mathrm{g}}_{10}\left(\frac{2}{T}{\int }_{t-T/2}^{t+T/2}\mathrm{}\frac{p^2\left(\tau \right)}{p_{\mathrm{ref}}^2}{w}_0(\tau -t)\mathrm{d}\tau \right),$$(32)

where p_ref = 2 · 10⁻⁵ Pa is the reference sound pressure and w₀ a Hann window with length T = 100 ms. To extract the level fluctuations F_p(t), the sound pressure level time history L_p(t) is detrended by substracting a weighted moving average:

$$F_p\left(t\right)=L_p\left(t\right)-\frac{2}{T^{\mathrm{\prime }}}{\int }_{t-\frac{T^{\mathrm{\prime }}}{2}}^{t+\frac{T^{\mathrm{\prime }}}{2}}\mathrm{}{L}_p\left({\tau }^{\mathrm{\prime }}\right)w_0^{\mathrm{\prime }}(t-\tau \mathrm{\prime })\mathrm{d}\tau \mathrm{\prime }.$$(33)

The moving average is computed by convolution of L_p(t) with a Hann window $w_0^{\text{'}}$ of length T′ = 15 s. Detrending the data removes the effect of the source directivity of the aircraft and the slowly changing, non-fluctuating part of the sound pressure due to changing propagation distance.

Finally, from F_p the standard deviation of log-amplitude fluctuations σ_χ (dB) is determined for the time frame from 10 s prior to 10 s after the moment of shortest distance between microphone and aircraft. The time frame length of 20 s allows to capture amplitude fluctuations with a modulation frequency of 0.25 Hz and higher while considering only the period of shortest distance to the microphone.

As the amplitude fluctuations develop during sound propagation through the turbulent atmosphere, the sound pressure signal is not backpropagated to the source by compensating for the propagation effects air absorption and geometrical spreading. The purpose of the de-Dopplerization is the removal of artefacts caused by tones. The effect of the de-Dopplerization on σ_χ is negligible.

Additionally, the modulation spectra of amplitude fluctuations are analyzed. For this, the time history of amplitude modulations is bandpass filtered before computing the standard deviation σ_χ. This represents an additional step between block 4 and block 5 in Figure 6. Thus, the standard deviation is obtained for the modulation bands from 0.25 Hz to 8 Hz in octave bands.

The analysis procedure has been verified by generating and analyzing amplitude modulated white noise signals. The presented processing of amplitude modulated signals can successfully extract the standard deviation and the modulation spectrum of the signal.

3.2 Meteorological data

Table 1 lists the meteorological parameters, that are considered in the analysis of amplitude modulations.

Meteorological information for the respective measurement location are available from the hourly reanalysis dataset ERA5 provided by the European Centre for Medium-Range Weather Forecasts [25]. The ERA5 dataset has a spatial resolution of approximately 31 km × 31 km. Additionally, information is available from the numerical weather prediction model COSMO-2 for 2014 and COSMO-1 for 2022, provided by the Swiss Federal Office of Meteorology and Climatology, MeteoSwiss [26]. COSMO-2 and COSMO-1 have a spatial resolution of 2 km × 2 km and 1 km × 1 km, respectively.

To characterize turbulence in the atmospheric boundary layer, the parameters temperature at 2 m height T_2m, surface temperature T_surf, boundary layer height z_i, surface sensible heat flux Q_H, and friction velocity $u_{*}$ are included in the analysis. For Q_H the ERA5 parameter “mean surface sensible heat flux” is chosen. Mean values for a given hour represent the mean of the preceding hour. Thus, the value is shifted by 30 min. The data is then linearly interpolated to the moment of flyover.

The ERA5 model can not resolve local meteorological differences caused by topography because of its coarse spatial resolution. Further, the temporal resolution of one hour is low considering the strong variability of wind related parameters and the short duration of an aircraft flyover. The boundary layer height z_i is expected to show only low spatial variability and to change slowly and continuously. Therefore, linear interpolation of the hourly data is sufficient. However, it is known that very low boundary layer heights during stable nocturnal boundary layers are not very accurate [27].

To investigate the validity of the temporal interpolation of $u_{*}$ and Q_H, in-situ turbulence measurements have been conducted with an ultrasonic anemometer at selected measurement locations, as shown in Figure 4. At these locations, surface sensible heat flux Q_H and friction velocity $u_{*}$ have been determined using the eddy-covariance method [28]. Figures 7a and 7b show comparisons of measured and estimated values of Q_H and $u_{*}$. The measured values are averages over a measuring period of 30 min. Measurements of Q_H show a good agreement with the shifted and interpolated data from ERA5. Measured friction velocity shows an approximate agreement with $2u_{*,\mathrm{ERA}5}$. The origin of the factor 2 is not entirely clear. The low spatial resolution might indeed not well represent the local turbulence. A different possible explanation is that the height of the ultrasonic anemometer was not sufficient to be always located within the surface layer of the ABL. Further, the expected error for the measured turbulent fluxes of momentum is estimated to be 10–50% for an average period of 1 hour, while the expected error for the turbulent flux of sensible heat is only 5% ([28], Chapter 7).

Figure 7 Comparison of measured data to interpolated hourly ERA5 dataset. (a) Surface sensible heat flux QH. (b) Friction velocity.

4 Analysis of empirical data

4.1 Development of amplitude fluctuations in the course of the day

Figure 8 gives an overview over the changing magnitude of amplitude fluctuations in the course of the day. Values for σ_χ range between 1 dB and 4 dB with a mean of 2.2 dB. These values are clearly below the theoretical limit of 〈χ²〉 = π²/24 corresponding to σ_χ ≈ 5.6 dB. The low values are a possible consequence of the broadband character of the signal. It is expected that broadband signals lead to lower log-amplitude fluctuations than narrowband signals, because the interference of the microrays is affected by the signal bandwidth [11]. Measurements of pulse propagation on the ocean surface showed that the scintillation index SI decreases with increasing signal bandwidth [29]. As a consequence, broadband signals approach saturation much slower than narrow-band signals [30].

Figure 8 Standard deviation σχ (dB) as a variable of time of the day. Colors indicate (a): the propagation path within the ABL LABL. Dark dots represent low LABL, light dots indicate high LABL. (b): u100m, wind speeds at 100 m height. Dark dots suggest presence of a low level jet.

The colorbar in Figure 8a represents the propagation path within the ABL L_ABL. Dark points represent short L_ABL and are associated with shallow, stable boundary layers. Light (yellow) dots represent long L_ABL which are associated with high and strongly convective, well-mixed boundary layers. The values of σ_χ clearly exhibit a diurnal cycle related to the atmospheric boundary layer. The colorbar in Figure 8b represents wind at 100 m height in the ERA5 model.

Figure 9 shows the standard deviation σ_χ as a variable of the parameters given in Table 1. The color indicates the time of the day in hours.

Figure 9 Standard deviation σχ (dB) depending on parameters given in Table 1. Colors indicate time of the day in hours.

During morning hours, σ_χ has on average lower values than during the rest of the day. In Figure 9 morning hours are represented by dark blue dots. Low values of $u_{*}$, Q_H, and z_i indicate shallow, weakly stable boundary layers in which turbulence is suppressed. During the late morning hours, the sun starts heating up the ground, resulting into growing heat flux Q_H from the warm ground towards the cooler air. Until noon, the stable boundary layer has developed into a turbulent convective boundary layer with growing atmospheric boundary layer heights. As the turbulent kinetic energy (TKE) in the atmosphere grows, σ_χ increases.

During afternoon hours, the boundary layer has reached its maximum height while σ_χ remains high. The boundary layer height z_i can exceed the altitude of the aircraft, resulting into high L_ABL (long propagation distances within the ABL) and hence high σ_χ. These cases are associated with the lightest dots in Figure 8a.

During early evening hours around 6 p.m., σ_χ and z_i still exhibit high values, although the intensity of the sun has decreased by now. Orange dots in Figure 9 approximately represent this period of the day. In Figures 9b and 9d can be noticed that σ_χ and z_i are high, while Q_H tends towards zero. The reason for the time delay is the time it takes the largest turbulent eddies to decay into smaller ones until they finally dissipate. The turbulence decay time $t^{\mathrm{*}}=\frac{z_i}{w_{\mathrm{*}}}$ is largest shortly before sunset [31]. In these periods of the day, Q_H is not a reliable predictor for σ_χ.

During sunset, a clear transition from light to dark dots in Figure 8a represents the boundary layer collapse, as no more large turbulent eddies are created by upward convection of warm air. The inverted heat flux from the warm air towards the cool ground suppresses turbulence. A shallow nocturnal stable boundary layer establishes. Usually, turbulent kinetic energy in the stable boundary layer is mainly caused by surface roughness leading to shear produced wind velocity fluctuations.

However, atmospheric turbulence and hence σ_χ can remain high after sunset. During the collapse of the convective boundary layer, a residual layer forms above the stable boundary layer up to the height of the convective boundary layer of the day. The residual can still exhibit strong turbulence, often due to nocturnal low-level jets (LLJ) [27, 32]. LLJ can appear at the beginning of the night on top of the stable boundary layer at heights between several tens and several hundreds of meters and are decoupled from wind at the ground. In fact, Figure 8b shows high values of σ_χ after 8 p.m., which are associated with high wind speeds at 100 m height in the ERA5 dataset. Further, the interaction of gravity waves between stably stratified layers and local topography can lead to production of turbulent kinetic energy [33].

4.2 Influence of geometry and meteorology

Friction velocity, shown in Figure 9a, shows a weak correlation with σ_χ. However, σ_χ quickly levels off for $u_{*}\gtrsim 0.05\mathrm{ }\mathrm{m}/\mathrm{s}$. Similarly, in Figure 9b, σ_χ levels off for $Q_H\gtrsim 100\mathrm{W}/{\mathrm{m}}^2$ and, in Figure 9c, for z_i ≳ 1 km. Figure 9d does not reveal a clear dependency between fluctuations and the mean propagation distance L_avg. Temperature difference between 5 cm and 2 m above ground (not shown) did not show to be a reliable predictor for σ_χ.

However, the sound propagation path within the ABL (L_ABL) was found to have a linear relationship with σ_χ as demonstrated in Figure 9e. This relationship combines information encapsulated in z_i and L_avg revealing a clear dependency on the sound propagation distance that is actually exposed to turbulent kinetic energy. Therefore, sound propagation above the boundary layer height has a significantly reduced impact on the development of amplitude fluctuations. Figure 9f presents the dependence on the ratio z_i/h_s of boundary layer height to sound source height. Levels quickly increase for z_i/h_s → 1 and stay constant if z_i exceeds the altitude of the aircraft.

4.3 D–ϕ diagram

Figure 10 presents a D–ϕ diagram of the flyover events for f = 1 kHz. For each data point, representing a single flyover, the wave parameter and the strength parameter correspond to the sum of their contributions from temperature and both wind velocity fluctuations: D = D_T + D_vs + D_vb and $\mathrm{\varphi }={\mathrm{\varphi }}_T+{\mathrm{\varphi }}_{\mathrm{vs}}+{\mathrm{\varphi }}_{\mathrm{vb}}$. The height dependence of ϕ and D are considered by integration of equations (20), (26), (23), and (27) along the propagation path:

$$D_{T,v}\left(h\right)=\frac{h}{k\cos \theta }{\int }_0^1\mathrm{}\frac{1}{L_{T,v}^2\left(\eta h\right)} \mathrm{d}\eta ,$$(34)

$${\mathrm{\varphi }}_T^2\left(h\right)=\frac{\sqrt{\pi }\mathrm{\Gamma }(\alpha -1)}{2\mathrm{\Gamma }(\alpha -3/2)}{k}^2\frac{h}{\cos \theta }{\int }_0^1\mathrm{}\frac{{\sigma }_T^2\left(\eta h\right)L_T\left(\eta h\right)}{T_0^2}\mathrm{d}\eta ,$$(35)

$${\mathrm{\varphi }}_v^2\left(h\right)=\frac{2\sqrt{\pi }\mathrm{\Gamma }(\alpha -1)}{\mathrm{\Gamma }(\alpha -3/2)}{k}^2\frac{h}{\cos \theta }{\int }_0^1\mathrm{}\frac{{\sigma }_v^2{L}_v\left(\eta h\right)}{c_0^2}\mathrm{d}\eta .$$(36)

Figure 10 D–ϕ diagram of aircraft flyover events for f = 1 kHz. Colored circles correspond to D and ϕ for each measured aircraft flyover. The colors indicate the measured standard deviation of amplitude fluctuations, σχ (dB). The black “–” signs represent ϕmax for each flyover, indicating the theoretical boundary between weak and strong scattering calculated by equation (30) with b2 = 0.5.

The values ϕ_max of each flyover event, calculated by equation (30), are indicated by the black “–” signs. ϕ_max estimates the theoretical boundary between the weak and strong-scattering regime for each measured flyover event. For a first estimate we chose b² = 0.5, corresponding to the value in [17]. In 3734 flyover measurements, corresponding to 70.4% of all measurements, ϕ exceeds ϕ_max for f = 1 kHz. Comparison to Figure 3 shows, that the measured amplitude fluctuations can be attributed to the partially-saturated regime.

4.4 Modulation spectrum

Figure 11 presents the normalized modulation spectrum of amplitude fluctuations. The largest modulations are observed in the octave bands 0.5 Hz and 1 Hz. Only events with σ_χ > 2.5 dB are considered, for which the results have the best signal-to-noise ratio. The modulation spectrum of each flight event is normalized by its overall standard deviation (STD) σ_χ. The vertical errorbars represent the STD of the normalized values. Fluctuations in the 4 Hz and 8 Hz octave band are very likely to be dominated by fluctuations in background noise.

Figure 11 Normalized modulation spectrum of amplitude fluctuations in octave bands. The modulation spectrum for each flight event is normalized by its overall standard deviation σχ (dB). Errorbars represent the STD of the normalized values.

4.5 Frequency dependency in empirical data

In order to test the frequency dependency of log-amplitude fluctuations, the data collected in 2022 is additionally analyzed for two distinct frequency ranges. The data is analyzed according to Section 3.1.2, this time applying smaller bandwidths of the bandpass filters. The lower frequency band encompasses the 500 Hz and the 630 Hz one-third octave band. The higher frequency band encompasses the 800 Hz and the 1000 Hz one-third octave bands.

Figure 12 compares the standard deviation of amplitude fluctuations of the two frequency bands. The results show no systematic difference between the two frequency bands. Theoretically, the relation between 〈χ²〉 and the acoustic wavenumber k for the Kolmogorov spectrum is 〈χ〉 ∼ k^7/12 ([19], Eq. 7.115). The Kolmogorov spectrum can be used instead of the von Kármán spectrum for log-amplitude variance, as the Kolmogorov spectrum is a good approximation of the von Kármán spectrum for the relevant turbulence wavenumbers. A possible reason is, that both frequency bands are located in the partially-saturated regime.

Figure 12 Comparison of standard deviation of log-amplitude fluctuations for two frequency bands. σχ,lowfreq (dB) is bandpass filtered for 446–707 Hz, corresponding to the one-third octave frequency bands with center frequencies 500 Hz and 630 Hz. σχ,highfreq (dB) is bandpass filtered for 708–1122 Hz, corresponding to the one-third octave frequency bands with center frequencies 800 Hz and 1000 Hz.

5 Semi-empirical model for amplitude fluctuations

5.1 Model description

The theoretical model described in Section 2.4 is taken as a basis and modified leveraging the empirical data presented in Sections 3 and 4. The variance of log-amplitude fluctuations 〈χ²〉 is described in the semi-empirical model by:

$$\widehat{\langle {\chi }^2\rangle }=\langle {\chi }_T^2\rangle \frac{1}{1+\frac{{\mathrm{\varphi }}_T^2}{{\mathrm{\varphi }}_{\max }^2}}+\langle {\chi }_{\mathrm{vs}}^2\rangle \frac{1}{1+\frac{{\mathrm{\varphi }}_{\mathrm{vs}}^2}{{\mathrm{\varphi }}_{\max }^2}}+\langle {\chi }_{\mathrm{vb}}^2\rangle \frac{1}{1+\frac{{\mathrm{\varphi }}_{\mathrm{vb}}^2}{{\mathrm{\varphi }}_{\max }^2}}.$$(37)

Here, $\langle {\chi }_{T,\mathrm{vs},\mathrm{vb}}^2\rangle$ are the log-amplitude fluctuations due to fluctuations in temperature, shear-produced and buoyancy-produced wind velocity given by equations (7) and (28) if only the respective source of turbulence is considered. ${\mathrm{\varphi }}_{T,\mathrm{vs},\mathrm{vb}}^2$ are the strength parameters for fluctuations in temperature, shear-produced and buoyancy-produced wind velocity given by equations (23) and (27). ${\mathrm{\varphi }}_{\max }^2$ is the maximum strength parameter in the weak-scattering regime, given by equation (30).

The semi-empirical model is based on the following five modifications of equation (28):

5.1.1 Consideration of turbulence decay time

First, at conditions of decreasing intensity of sun, where Q_H decreases much faster than both z_i and the TKE in the ABL, a temporal offset for the value of Q_H is used, based on the turbulence decay time $t^{*}=z_i/w_{*}$. For a gradual decrease of surface sensible heat flux Q_H typical for afternoon hours, the kinetic turbulent energy in the ABL was simulated in [31]. The TKE remains constant for a time period $t^{\text{'}}\le t^{*}$. For $t^{*}< t^{\text{'}}\lesssim 3t^{*}$ the TKE decreases slightly, for $t^{\text{'}}>4t^{*}$ the TKE shows a decay proportional to t′⁻². Hence, for Q_H the following condition is introduced:

$$Q_H\left(t\right)=\{\begin{array}{ll}{Q}_H\left(t\right)& \mathrm{for} \frac{\mathrm{d}{Q}_H}{\mathrm{d}t}\ge 0\\ Q_H(t-3t^{\mathrm{*}})& \mathrm{otherwise}.\end{array}$$(38)

5.1.2 Neglect buoyancy-produced velocity fluctuations in SBL

Second, if Q_H ≤ 0, the term for buoyancy-produced velocity fluctuations in the effective turbulence spectrum, equation (37), is set to 0:

$$\begin{array}{ccc}{\sigma }_{\mathrm{vb}}^2=0& \mathrm{for}& Q_H<0.\end{array}$$(39)

For Q_H ≈ 0, the state of the ABL is near neutral and TKE is mainly generated by wind shear. Consequently, the nominator in equation (10) tends towards 0 for buoyancy-produced wind velocity fluctuations. For Q_H < 0, the state of the ABL is stable. Turbulence is suppressed by the negative surface sensible heat flux. The model presented in Section 2.4 is not valid for SBL. Therefore, SBL are approximated by neglecting the term for buoyancy-produced wind velocity fluctuations in equation (37).

5.1.3 Low level jet approximation

Third, generation of TKE at the top of the boundary layer by nocturnal low-level-jets is approximated by replacing the friction velocity $u_{*}$ with

$${\tilde{u}}_{\mathrm{*}}=(0.1u_{100\mathrm{m}}+u_{\mathrm{*}})/2,$$(40)

if the following conditions are fulfilled: 0.1u_100m > 3 $u_{*}$ and Q_H < 10 W/m². The model described in Section 2.4 does not consider TKE generation by nocturnal low-level-jets, as it is based on MOST and mixed-layer similarity theory. The chosen approximation is not an exact representation of the turbulence mechanism, but it is able to approximate the results of the field measurements.

5.1.4 Effective propagation distance

Fourth, the propagation distance L is replaced by the effective propagation distance L_eff based on L_ABL and the ratio of z_i to the height h_s of the aircraft:

$$L_{\mathrm{eff}}=\{\begin{array}{ll}{L}_{\mathrm{ABL}}{\left(\frac{z_i}{h_s}\right)}^{-q\left(1-\frac{z_i}{h_s}\right)}& \mathrm{for} h_s>z_i\\ L& \mathrm{otherwise},\end{array}$$(41)

with the empirical coefficient q = 0.6. The dependency of σ_χ on the distance propagated within the ABL instead of the actual propagation distance was discussed in Section 4. The top of stable boundary layers is often lower than the altitude h_s of an aircraft. The effective propagation distance is approximated as a combination of L_ABL and the ratio z_i/h_s. For z_i ≈ h_s, the correction is minimal. For z_i/h_s → 0, L_eff is much greater than L_ABL.

The empirical coefficient q determines how much of the sound propagation path located above the ABL is considered as propagation distance contributing to the generation of log-amplitude fluctuations.

This correction targets two sources of TKE outside of the boundary layer. If the sound source is located above the boundary layer, sound waves propagate through the residual layer located on top of the SBL, which is not entirely turbulence free [27] and therefore contributes to the generation of amplitude fluctuations. Further, shear-produced turbulence is not only created close to the ground. An additional source of TKE during the day is the entrainment zone, which is the transition between the convective boundary layer and the air above [20]. Air from the residual layer or the free troposphere is entrained in the boundary layer, creating an entrainment flux of temperature and momentum.

In the semi-empirical model, the strength and the wave parameters are evaluated for the effective propagation distance L_eff = h_eff/cos θ by replacing h by h_eff in equations (34)–(36).

5.1.5 Considering partial saturation

The representation of 〈χ²〉 in terms of strength and wave parameters allows to describe the boundary between weak and strong-scattering regime using the parameter ϕ_max (see Eq. (30) in Sect. 2.6). For ϕ_T,v ≪ ϕ_max, the propagated sound is only weakly scattered by atmospheric turbulence and log-amplitude fluctuations are well described by equations (7) and (28), as shown in [10] by validation measurements. For ϕ_T,v ≥ ϕ_max, the model overpredicts log-amplitude fluctuations. Therefore, we use ${\mathrm{\varphi }}_{\max }$ to scale $\langle {\chi }_{T,v}^2\rangle$ by:

$${\widehat{\langle {\chi }_{T,v}^2\rangle }}_{=}\langle {\chi }_{T,v}^2\rangle \frac{1}{1+\frac{{\mathrm{\varphi }}_{T,v}}{{\mathrm{\varphi }}_{\max }}}.$$(42)

Equation (42) and the coefficient b in equation (30) are empirically determined, with

$$b^2=1/6.$$(43)

In [17] the coefficient was set to b² = 1/2, which is higher than the empirically determined value. However, a direct comparison might be misleading, as the empirical data is broadband while [17] present theory for mono-frequency signals.

5.2 Model result

Figure 13 compares the diurnal cycles of the measured and modeled amplitude fluctuations.

Figure 13 Left: Measured σχ,meas (dB). Middle: Semi-empirical model for σχ,model (dB). Right: Difference of measured and modeled amplitude fluctuations σχ,meas − σχ,,model. The color scale indicates the effective propagation distance Leff.

The broadband character of the signal is approximated by taking the average of linearly spaced frequencies in the frequency range 446–1120 Hz. As mentioned in Section 4, it is expected that broadband signals lead to lower log-amplitude fluctuations than narrowband signals, because the interference of the microrays is affected by the signal bandwidth. These effects are not considered in this article and should be subject of further model extensions.

5.3 Determination of empirical coefficients

The coefficients b and q were determined based on the empirical dataset. Smaller values for b result into overall lower model results, as the boundary between weak and strong scattering is decreased. For b → 1, model results increase.

The coefficient q = 0.6 in the exponent of equation (41) determines the slope of the diurnal cycle. For q → 0, σ_χ becomes very low for small boundary layer heights. Especially cases with shallow stable boundary layers then result into very low amplitude fluctuations σ_χ. For q → 1, cases with relatively shallow boundary layers but high traveling altitude of the aircraft are overestimated.

Figure A1 in the Appendix showcases results of a parameter study for the model parameters q and b.

6 Model evaluation

6.1 Morning hours

The semi-empirical model shows the largest uncertainties for stable boundary layers in the morning hours. Predictions during morning and night are lower than the fluctuations measured during the field measurements. During the field measurements, σ_χ rarely fell below the limit of 1 dB. Turbulence-free atmospheric conditions occur only in rare cases. Usually, intermittent turbulence in the residual layer above the stably stratified nocturnal boundary layer is still able to cause fluctuations in the acoustical amplitude. However, the absence of very low σ_χ can be partly attributed to unavoidable background noise in the measurements, i.e. sound level fluctuations due to sound sources other than the aircraft. Therefore, overfitting of the model to the data should be avoided, especially in weak turbulence conditions. Additionally, these situations are often associated with Q_H < 0, which is outside of the range of applicability for MOST and mixed layer similarity theory. Similarly, the validations of the theoretical model in [10] showed clear deviations for measurements during early morning hours, while the model agreed well with measurements during daytime. With these considerations in mind, the values of the coefficient b in equation (30) and of q in the exponent of equation (41) have been determined based on the measurement data.

6.2 Comparison to weak-scattering model

Figure 14 compares the the measurements with the results of the theoretical model given by equation (7).

Figure 14 Comparison of measured and modeled σχ (dB). Left: Measured values. Right: Theoretical weak-scattering model according to equation (7). In most aircraft flyover cases, the limits of applicability of the weak-scattering model are violated.

Using the theoretical model to estimate the amplitude fluctuations for the aircraft flyover conditions leads to a mean absolute error (MAE) of 2.9 dB in comparison to the measurements. In most cases, the meteorological conditions are beyond the limits of applicability, i.e. sound propagation above the ABL, strong scattering, or stable boundary layers. The semi-empirical model has a MAE of 0.4 dB. The semi-empirical model, incorporating the modifications presented in equations (38)–(43), significantly improves the prediction of amplitude fluctuations.

6.3 Expected uncertainty of the model

The accuracy of the semi-empirical model is very good considering the inherent uncertainties in the parameters describing atmospheric turbulence. The standard deviation of $u_{*}$ and Q_H is often of the same order of magnitude as the value itself (Tab. 1 in [10]). To estimate the expected differences between two flight events in almost identical conditions, we determined the difference in σ_χ between any two consecutive flights at the same measurement site with less than 10 min time difference. For all measurement sites with sufficient data, the differences are normally distributed around 0 dB with a standard deviation of 0.35–0.47 dB (not shown). This standard deviation is of the order of the MAE of the semi-empirical model. Further adaptations of the model to decrease the MAE are likely to lead to overfitting of the model to the data.

6.4 Spatial resolution of meteorological data

According to the analysis, the ERA5 data was found to be a better input for the model compared to COSMO-1 and COSMO-2, which have considerably higher spatial resolutions. Using local meteorological input data did slightly increase the mean absolute error of the semi-empirical model. This can be explained by local meteorological effects caused by topography. Local meteorological data is not the most suitable characterization of sound propagation with relatively long distances through the atmospheric boundary layer in non-ideally flat environments. Instead, the ERA5 model’s low spatial resolution demonstrated a greater ability to predict amplitude fluctuations induced by turbulence. This enhances the model’s usability since all input data can be obtained from an open-access online databases, eliminating the need for in-situ measurements and enabling the modeling of past events.

6.5 Limitations

The semi-empirical model is based on data collected in meteorological conditions without rainy events and without extreme weather conditions like temperatures below 0 °C and wind speeds higher than 7 m/s. Further, the dataset contains only flyover events with propagation paths between L_avg = 570 m and L_avg = 3250 m.

The semi-empirical model does not consider the gradual approach to the full-saturation regime, as depicted in Figure 2. For much larger distances, e.g. time frames before and after the considered flyover of 20 s or much higher altitudes of the aircraft, the intensity of the fluctuation is expected to decrease. Therefore, the semi-empirical model might not give plausible results for approaching and departing aircraft if relatively long time frames are considered. The semi-empirical model targets the flyover.

The model is based on the analysis of moving sound sources. This limits the generalization of the modulation spectrum to similar applications. Stationary sources are not characterized by the modulation spectrum shown in Figure 11. The generation of sound by the interaction of the moving sound source with the atmospheric turbulence was neglected in the analysis. Further studies are necessary to allow a separation between fluctuations generated at the source and generated by propagation through turbulence. In the present analysis, only narrow-body airliners were considered.

7 Application to auralization

To demonstrate the use of the proposed semi-empirical model for amplitude fluctuations described in Section 5, two flyovers of an aircraft were synthesized under different atmospheric conditions. The synthesis methods are described in [7, 34]. The emission and flight path data of a departing Airbus A320-214 are taken from a previous study [7] (denoted as validation case A). Here, buzzsaw and fan tonal noise emission are attenuated by 10 dB and 5 dB, respectively, to better match the measurements and provide a more realistic impression of the source sound characteristics. The receiver is located close to the noise certification point for takeoffs and at h_r = 10 m above grassy ground. During direct flyover, the aircraft is at an altitude of h_s = 700 m above ground. The flight trajectory of the synthesized flights approximately matched those of the real-world comparison cases. Two effects of atmospheric turbulence on sound propagation were simulated. Turbulence-induced coherence loss in ground effect was modeled as described in [6] using a time-variant decorrelation filter. Amplitude fluctuations were modeled by a time-variant high shelving filter with variable cutoff frequency and gain [34]. The cutoff frequency of the shelving filter and its gain function were steered by the semi-empirical model.

Because both, the distance L and the source height h_s, vary during the flyover, also the cutoff frequency changes over time. During the flyover event, the cutoff frequency is computed as the lowest frequency for which the inequality

$$\max ({\mathrm{\varphi }}_T^2,{\mathrm{\varphi }}_{\mathrm{vs}}^2,{\mathrm{\varphi }}_{\mathrm{vb}}^2)\ge {\mathrm{\varphi }}_{\max }^2$$(44)

is true.

Two different meteorological conditions (Condition I and Condition II) were chosen from the dataset reflecting realistic situations. The situations characterize a moderately weak-turbulence condition during morning and a strong-turbulence condition at noon. Spectrograms of the measurements are presented in Figure 1 in Section 1. Parameters are given in Table 2. The corresponding mean amplitude fluctuations for the time frames from 10 s prior to 10 s after the moment of shortest distance between microphone and aircraft, determined with the semi-empirical model, are σ_χ = 1.1 dB in Condition I and σ_χ = 2.7 dB in Condition II. Because of the distance dependence of the amplitude fluctuation, σ_χ is higher at the beginning and end of the synthesized flyover and lowest in the moment of direct flyover.

Figure 15 shows the resulting spectrograms of the syntheses for these two conditions. It can be observed that the strength of the vertical pattern caused by the fluctuations are very similar to those in the measurements shown in Figure 1.

Figure 15 Synthesized aircraft flyover events for two sets of meteorological parameters shown in Table 1 corresponding to the meteorological conditions of the examples in Figure 1.

Audio files are provided on Zenodo [35]. The provided audio files are listed in Table 3. The baseline synthesis is the auralization without application of the amplitude modulation filter.

The difference of the fluctuation intensity between the low-turbulence and high-turbulence case is clearly perceivable in the measurements as well as in the syntheses. In the meteorological condition II, the fluctuations are very prominent. Comparing the high-turbulence synthesis (1 ) to the high-turbulence measurement (2 ) shows a very good agreement of the fluctuation’s sound characteristics. For the low-turbulence synthesis (3 ) and low-turbulence measurements (4 ), the turbulent fluctuations are considerably softer than in the high-turbulence condition. The characteristics of the low-turbulence fluctuations are well matched in the low-turbulence synthesis. Comparing the low-turbulence condition synthesis to the turbulence-free baseline synthesis (5 ) reveals that the low fluctuation is still perceivable and does increase the realism of the synthesis compared to the measurement. A listening experiment is planned for the future to confirm the chosen modifications.

8 Conclusion

A semi-empirical model was presented to predict amplitude fluctuations in aircraft sound propagated through atmospheric turbulence. The model is based on the theoretical model by Ostashev and Wilson [9] and on empirical data. The empirical dataset comprises a large number of aircraft flyovers in a wide range of meteorological conditions, information on the respective flight trajectories and meteorological data. The analysis of the empirical data revealed a strong dependence of σ_χ, the standard deviation of amplitude fluctuations, on the length of the propagation path within the atmospheric boundary layer, denoted as L_ABL.

The semi-empirical model incorporates five main modifications compared to the theoretical model. First, considering turbulence decay times allows predictions of amplitude fluctuations in the afternoon, when the decrease in solar intensity and hence surface sensible heat flux otherwise leads to a mismatch between turbulence scaling parameters and the turbulent kinetic energy in the atmosphere. Second, buoyancy-produced velocity fluctuations are neglected for stable boundary layers with Q_H < 0 W/m². Third, turbulence generation by nocturnal low-level jets is approximated based on wind speed at 100 m height. Fourth, the propagation distance L is replaced by an effective propagation distance based on the boundary layer height in relation to the source height. Fifth, frequency- and distance dependence partial saturation of log-amplitude fluctuations is considered based on the theoretical boundary between weak and strong scattering. The boundary is determined using a representation of log-amplitude variance based on strength and wave parameters.

Applying the modifications to the theoretical model decreased the mean absolute error from 2.9 dB to 0.4 dB for the considered dataset. In contrast to the purely theoretical weak-scattering model, the semi-empirical model does account for a saturation effect of amplitude fluctuations. In the weak-scattering regime, i.e. weak turbulence or relatively short propagation distances, the semi-empirical model coincides with the theoretical model by [9].

The analysis showed that local meteorological conditions can be inappropriate to characterize sound propagation through atmospheric turbulence from an elevated sound source at relatively large distances. The model’s error slightly increased when using the numerical weather model COSMO-1 with a 1 km × 1 km grid resolution compared to using the coarser ERA5 model, which has a 31 km × 31 km resolution. The use of the ERA5 dataset, which provides open access to data from several decades in the past, enhances the model’s usability.

The application of the semi-empirical model for aircraft auralization was demonstrated for flyovers of a departing airliner at about 700 m above ground. Two flyovers were synthesized under different atmospheric conditions. Turbulent amplitude fluctuations were modeled by time-variant filters allowing to model the distance dependence of the turbulence effect. The meteorological conditions considered in the syntheses are based on real-world flyover events. The syntheses demonstrated a high level of agreement for both the low-turbulence and high-turbulence condition. The audio files are provided on an open-access online repository [35].

Listening experiments are planned for the near future to provide a rigorous assessment of the semi-empirical model and its application to aircraft auralization.

The semi-empirical model will facilitate the modeling of sound propagation through atmospheric turbulence for elevated sound sources like aircraft, wind turbines and drones. Listening experiments using auralizations of outdoor sound sources will benefit from improved realism. In addition, the semi-empirical model can be used to estimate the measurement uncertainty of measured aircraft noise, especially for maximum noise levels.

Acknowledgments

The authors appreciate that SWISS International Airlines supplied their flight data recorder data for this study. We thank Zurich airport for providing radar data, Lothar Bertsch from the German Aerospace Center DLR for providing the aircraft noise emission data, and the colleagues at Empa who conducted the field measurements. We thank the anonymous reviewers for their very helpful comments and feedback. In particular, we would like to thank the reviewer who provided us with very useful feedback and references on saturation of log-amplitude fluctuations.

Funding

This research was supported by the international research project Localization and Identification Of moving Noise sources (LION) in which the Swiss contribution was funded by the Swiss National Science Foundation (SNSF) by Grant No. 185530.

Conflicts of interest

The authors declare that there is no conflict of interest.

Data availability statement

The sound files associated with this article are available on Zenodo, under the reference [35].

Appendix

Figure A1 shows the effect of changing the empirical coefficients b in equations (43) and (30), and q in the exponent of equation (41).

Figure A1 Parameter study for empirical coefficients b and q of the semi-empirical model.

References

All Tables

All Figures

	Figure 1 Spectrograms of measured Airbus A320 flyover during take-off at approximately 650 m above ground. Colors indicate power spectral density (PSD). Top: Low turbulence case at 7 a.m. Bottom: High turbulence case at 12 a.m.
In the text

	Figure 2 Qualitative dependence of scintillation index with range. The range at which full saturation occurs is frequency dependent. Red area corresponds to behavior that is expected for aircraft flyover. Adapted from Figure 4.9 in [11].
In the text

	Figure 3 Qualitative D–ϕ diagram showing unsaturated, partially saturated and fully saturated regime depending on strength parameter ϕ and wave parameter D. Adapted from Figure 4.10 in [11] and Figure 5 in [17].
In the text

	Figure 4 Mounted ultrasonic anemometer (front) and microphone (back).
In the text

	Figure 5 Distances and polar angles θ of measured aircraft flyovers for the moment of shortest distance between aircraft and microphone (θ = 0 is vertical). Blue (darker) points represent landings, Red (lighter) points represent take-offs.
In the text

	Figure 6 Flowchart of audio signal processing to estimate the amplitude fluctuation indicator σχ (dB) from measured aircraft flyovers
In the text

	Figure 7 Comparison of measured data to interpolated hourly ERA5 dataset. (a) Surface sensible heat flux QH. (b) Friction velocity.
In the text

	Figure 8 Standard deviation σχ (dB) as a variable of time of the day. Colors indicate (a): the propagation path within the ABL LABL. Dark dots represent low LABL, light dots indicate high LABL. (b): u100m, wind speeds at 100 m height. Dark dots suggest presence of a low level jet.
In the text

	Figure 9 Standard deviation σχ (dB) depending on parameters given in Table 1. Colors indicate time of the day in hours.
In the text

	Figure 10 D–ϕ diagram of aircraft flyover events for f = 1 kHz. Colored circles correspond to D and ϕ for each measured aircraft flyover. The colors indicate the measured standard deviation of amplitude fluctuations, σχ (dB). The black “–” signs represent ϕmax for each flyover, indicating the theoretical boundary between weak and strong scattering calculated by equation (30) with b2 = 0.5.
In the text

	Figure 11 Normalized modulation spectrum of amplitude fluctuations in octave bands. The modulation spectrum for each flight event is normalized by its overall standard deviation σχ (dB). Errorbars represent the STD of the normalized values.
In the text

Figure 12 Comparison of standard deviation of log-amplitude fluctuations for two frequency bands. σχ,lowfreq (dB) is bandpass filtered for 446–707 Hz, corresponding to the one-third octave frequency bands with center frequencies 500 Hz and 630 Hz. σχ,highfreq (dB) is bandpass filtered for 708–1122 Hz, corresponding to the one-third octave frequency bands with center frequencies 800 Hz and 1000 Hz.

In the text

	Figure 13 Left: Measured σχ,meas (dB). Middle: Semi-empirical model for σχ,model (dB). Right: Difference of measured and modeled amplitude fluctuations σχ,meas − σχ,,model. The color scale indicates the effective propagation distance Leff.
In the text

	Figure 14 Comparison of measured and modeled σχ (dB). Left: Measured values. Right: Theoretical weak-scattering model according to equation (7). In most aircraft flyover cases, the limits of applicability of the weak-scattering model are violated.
In the text

	Figure 15 Synthesized aircraft flyover events for two sets of meteorological parameters shown in Table 1 corresponding to the meteorological conditions of the examples in Figure 1.
In the text

	Figure A1 Parameter study for empirical coefficients b and q of the semi-empirical model.
In the text