4.1 General description of the acoustic field

The three-dimensional normalized pressure fields $\unicode[STIX]{x1D6F7}$ obtained at different instants of time are displayed in figure 7(a–i). The black box represents the computational moving frame, whereas the red box indicates the physical domain of interest. The colour level ranges from $-10~\text{Pa}~\text{kg}^{-1/2}~\text{m}^{3/2}$ (blue) to $+10~\text{Pa}~\text{kg}^{-1/2}~\text{m}^{3/2}$ (red). This interval is chosen in order to visualize the various infrasonic phases, which may have very different amplitudes. Successive points on the axes $x_{1}$ , $x_{2}$ and $x_{3}$ are $100$  km distant.

Figure 7.

Three-dimensional normalized pressure field $\unicode[STIX]{x1D6F7}$ obtained for different instants of time. The black box represents the computational moving frame, whereas the red box is the physical domain of interest. The times are (a)  $t_{1}=146.79$  s, (b)  $t_{2}=216.69$  s, (c)  $t_{3}=356.49$  s, (d)  $t_{4}=566.19$  s, (e)  $t_{5}=775.89$  s, (f)  $t_{6}=1055.49$  s, (g)  $t_{7}=1265.19$  s, (h)  $t_{8}=1404.99$  s, (i)  $t_{9}=1544.79$  s.

Near the source location, a spherical wavefront is observed. However, as a result of the vertical gradients of the speed of sound and of the horizontal wind, the infrasonic wave is continuously deformed during the propagation. At the instant $t_{2}$ , partial reflections start being observed between $25$ and $55$  km altitude. They are induced by variations of the speed of sound and of the wind of length scale comparable to the acoustic wavelength. The amplitude of these reflections appears strengthened at the instant $t_{3}$ . Moreover, because of nonlinear effects, the acoustic wavefront lengthens while propagating towards the thermosphere. At the instant $t_{4}$ , the infrasonic wave is refracted back towards the Earth’s surface. In addition, part of the wavefront leaves the domain through the top boundary. The computational frame starts moving between times $t_{4}$ and $t_{5}$ . At the instant $t_{5}$ , the signal at ground level, between $x_{1}=200$  km and $x_{1}=250$  km on the $x_{1}$ axis, is a superposition of partial reflections and stratospheric returns, while a thermospheric phase is reaching the Earth’s surface at about $x_{1}=180$  km distance from the source location. Furthermore, part of the wavefront leaves the moving window through the plane $x_{2}=-L_{2}$ with a negligible level of spurious reflection. At the instant $t_{6}$ , the thermospheric phase is clearly visible on the ground at around $x_{1}=250$  km. At successive instants, the pressure field becomes more and more complex and a multitude of arrivals are recorded at the Earth’s surface. The computational box stops moving just after reaching the boundary of the physical domain, between times $t_{6}$ and $t_{7}$ . The simulation ends up once the wavefront exits the frame.

4.2 The signal in the near field

The pressure perturbation $p^{\prime }$ recorded on the $x_{1}$ axis at $x_{1}=7.92$  km and the corresponding energy spectral density ${\mathcal{E}}_{p^{\prime }}$ are plotted in figure 8(a,b). The pressure signal in the near field exhibits an N-waveform, has a maximum amplitude of $2330$  Pa and a total duration of $6$  s. The central frequency is equal to about $0.13$  Hz and corresponds to a wavelength of $2.6$  km.

Figure 8.

(a) Pressure perturbation $p^{\prime }$ recorded on the $x_{1}$ axis at $x_{1}=7.92$  km and (b) corresponding energy spectral density ${\mathcal{E}}_{p^{\prime }}$ .

4.3 Nonlinear effects in the upper atmosphere

The normalized pressure perturbations $\unicode[STIX]{x1D6F7}$ recorded along the vertical axis $x_{3}$ at three different altitudes are plotted in figure 9(a–c). At $x_{3}=50$  km, an N-shaped signal with a period of $12$  s is observed. As a result of nonlinearities, the N-wave lengthens while propagating towards the upper atmosphere, its duration being equal to about $26$ and $75$  s at $x_{3}=90$ and $130$  km altitude, respectively. A diminution of the amplitude associated with the increase of the wave period is observed as well. Moreover, at $x_{3}=130$  km, the central part of the N-wave is found to be curved. This result is in disagreement with the weakly nonlinear ray theory, according to which an N-wave should conserve its shape during the propagation. The weakly nonlinear ray theory rests upon the assumption that the amplitude of the signal $p^{\prime }$ remains at least one order of magnitude smaller than the atmospheric pressure $\bar{p}$ , that is, $p^{\prime }<0.1\,\bar{p}$ (for instance, see Whitham (Reference Whitham1974)). As illustrated in figure 10, where the maximum in time of the ratio $P_{3}(x_{3},t)=p^{\prime }(x_{1}=0,x_{2}=0,x_{3},t)/\bar{p}(x_{3})$ is plotted as a function of the altitude $x_{3}$ , the hypothesis $\max _{t}(P_{3})<0.1$ is however not verified near the source as well as in the thermosphere. As an example, $\max _{t}(P_{3})$ is equal to about $0.7$ at $x_{3}=114$  km altitude. Therefore, the nonlinear effects cannot be considered weak in the present simulation. To the best of the authors’ knowledge, a model for moderate amplitude waves in a variable atmosphere is currently not available. Qualitative insights into the understanding of the aforesaid signal distortion can nevertheless be gained through a closer analysis of the propagation of one-dimensional plane waves in a homogeneous medium at rest. First, let it be assumed that the viscous and the thermal conduction terms are negligibly small. Second, with $x_{3}$ the propagation direction, let it be supposed that the specific entropy $s$ is constant, $s=\bar{s}$ , and that the velocity $u_{3}=u_{3}^{\prime }$ is a single-valued function of the pressure $p$ . According to Whitham (Reference Whitham1974), the evolution of a right-running signal is then described by the following nonlinear partial differential equation:

(4.2) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}t}+(u_{3}+c)\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}x_{3}}=0,\end{eqnarray}$$

with $p=\bar{p}(\unicode[STIX]{x1D70C}/\bar{\unicode[STIX]{x1D70C}})^{\unicode[STIX]{x1D6FE}}$ , $c=\sqrt{\unicode[STIX]{x1D6FE}p/\unicode[STIX]{x1D70C}}$ and $u_{3}=2(c-\bar{c})/(\unicode[STIX]{x1D6FE}-1)$ . Writing the sum $(u_{3}+c)$ as a function of the pressure $p$ and expanding it for $p^{\prime }/\bar{p}\rightarrow 0$ yields

(4.3) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}t}+\bar{c}\left[1+\frac{\unicode[STIX]{x1D6FD}}{\bar{\unicode[STIX]{x1D70C}}\bar{c}^{2}}\,p^{\prime }-\frac{\unicode[STIX]{x1D6FD}^{2}}{2\bar{\unicode[STIX]{x1D70C}}^{2}\bar{c}^{4}}\,{p^{\prime }}^{2}+\text{h.o.t.}\right]\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}x_{3}}=0,\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}=(\unicode[STIX]{x1D6FE}+1)/2$ is the nonlinear parameter. In the presence of shock waves, the previous hypotheses on the specific entropy $s$ and on the velocity $u_{3}$ are not verified. However, it can be shown that, at a discontinuity, the relations $s-\bar{s}=0$ and $u_{3}=u_{3}^{\prime }=2(c-\bar{c})/(\unicode[STIX]{x1D6FE}-1)$ are valid up to the third order in $p^{\prime }/\bar{p}$ when $p^{\prime }/\bar{p}$ tends to $0$ (Whitham Reference Whitham1974). Therefore, equation (4.3) can still be employed in the presence of small or moderate amplitude shocks. When only the first two terms of the Taylor series for $(u_{3}+c)$ are retained, $u_{3}+c=\bar{c}+\unicode[STIX]{x1D6FD}p^{\prime }/(\bar{\unicode[STIX]{x1D70C}}\bar{c})$ , the classical Burgers equation is recovered. Ray theory is entirely based on an extension of this latter equation to inhomogeneous media (for instance, see Sabatini et al. (Reference Sabatini, Marsden, Bailly and Gainville2016b ) and references therein). In the case of the Burgers equation, an N-wave conserves its shape as a consequence of the linearity in $p^{\prime }$ of the first-order approximation of $(u_{3}+c)$ . The central part of the waveform remains straight and only a variation of slope can be observed during the propagation. The curvature of the signals recorded in the thermosphere might be thus attributed to high-order terms in $p^{\prime }$ , which cannot be neglected beyond a certain altitude. A similar discussion of the distortion of high-amplitude N-waves can be found in Inoue & Yano (Reference Inoue and Yano1997).

Figure 9.

Normalized pressure signals $\unicode[STIX]{x1D6F7}$ recorded on the $x_{3}$ axis at (a)  $x_{3}=50$  km, (b)  $x_{3}=90$  km and (c)  $x_{3}=130$  km.

Figure 10.

Maximum of the ratio $P_{3}$ between the pressure perturbation $p^{\prime }$ on the $x_{3}$ axis and the atmospheric pressure $\bar{p}$ . The dashed line represents the conventional limit of validity of the weakly nonlinear ray theory, $\max _{t}(P_{3})=0.1$ .

Figure 11.

Spectrum ${\mathcal{E}}_{\unicode[STIX]{x1D6F7},3}$ of the signals $\unicode[STIX]{x1D6F7}$ recorded along the $x_{3}$ axis.

The energy spectral density ${\mathcal{E}}_{\unicode[STIX]{x1D6F7},3}(x_{3},f)={\mathcal{E}}_{\unicode[STIX]{x1D6F7}}(0,0,x_{3},f)$ of the signals $\unicode[STIX]{x1D6F7}$ recorded along the $x_{3}$ axis is finally plotted in figure 11 as a function of the frequency $f$ and of the altitude $x_{3}$ . The nonlinear lengthening of the N-wave leads to a shift of the spectrum towards very low frequencies. Therefore, throughout the propagation, most of the energy is carried by spectral components which are weakly affected by the viscous and the thermal conduction terms. Dissipative effects only influence the shock region, but play a minor role in the waveform.

Figure 12.

Acoustic field $\unicode[STIX]{x1D6F7}$ on the $x_{1}{-}x_{3}$ plane at (a) $t_{3}$ , (c) $t_{5}$ and (e)  $t_{7}$ . Pressure signals $p^{\prime }$ recorded on the $x_{1}$ axis at (b)  $x_{1}=100$  km, (d)  $x_{1}=230$  km and (f)  $x_{1}=360$  km. The grey lines in (a,c) represent the acoustic rays.

4.4 Description of the arrivals at ground level

Two waveguides are induced by the vertical profiles (2.1) of the speed of sound and of the wind, a stratospheric duct between the Earth’s surface and about $55$  km altitude, and a thermospheric duct between the ground and $115$  km altitude. As a consequence, different arrivals reach the Earth’s surface. The pressure perturbations $p^{\prime }$ detected on the $x_{1}$ axis at the ground stations $x_{1}=100$  km, $x_{1}=230$  km and $x_{1}=360$  km are plotted in figures 12(b), 12(d) and 12(f), respectively. In order to help the understanding of the various arrivals, the acoustic fields $\unicode[STIX]{x1D6F7}$ obtained on the $x_{1}{-}x_{3}$ plane at times $t_{3}$ , $t_{5}$ and $t_{7}$ are also displayed in figures 12(a), 12(c) and 12(e), respectively. The results obtained by ray theory (Candel Reference Candel1977) are reported as well in figure 12(a,c). At the ground station located at $x_{1}=100$  km range, refer to figure 12(b), two distinct wave packets are observed. The first arrival, labelled $I_{dtr}$ , where the subscripts $d$ and $tr$ stand respectively for diffracted and tropospheric, is not predicted by the geometrical acoustics approximation and may be attributed to a creeping wave which diffracts energy along the Earth’s surface (Pierce Reference Pierce1985). It is detected between about $t=290$  s and $t=300$  s, its maximum overpressure is equal to $30.8$  Pa and its duration is around $6$  s. No period lengthening is observed, meaning that the propagation at ground level is essentially linear. An analogous conclusion has been drawn in a recent study by de Groot-Hedlin (Reference de Groot-Hedlin2017). The second arrival, detected from $t=360$ to $460$  s, is associated with the partial reflections the wavefront undergoes while travelling within the small-scale inhomogeneous layer located between $25$ and $55$  km altitude. For this reason, it is denoted $I_{ps}$ , where the subscripts $p$ and $s$ stand for partial reflection and stratospheric. Its amplitude is about $7$  Pa. It can also be remarked that, at time $t_{1}$ , a stratospheric return, labelled $I_{s}$ , becomes visible in figure 12(a) at about $50$  km altitude. At the ground station located at $x_{1}=230$  km on the $x_{1}$ axis, three different wave packets are recorded along with the partial reflections $I_{ps}$ . The stratospheric phase $I_{s}$ is observed between $780$ and $820$  s, with an amplitude equal to about $73$  Pa. Furthermore, contrarily to the prediction of ray theory, a thermospheric wave packet $I_{t}$ is also detected between $980$ and $1100$  s, with a maximum overpressure of $3$  Pa. As discussed in § 4.7, the arrival $I_{t}$ is due to diffraction at the thermospheric caustic. At $t=t_{3}$ and $t=t_{5}$ , arch-like structures are observed between the ground and the lower thermosphere. They are associated with partial reflections occurring whenever the wavefront propagates downward or upward through the small-scale inhomogeneous layer. The resulting phase, labelled as $I_{pst}$ , is clearly visible at the ground stations located at $x_{1}=230$  km and at $x_{1}=360$  km. Finally, zooms of the stratospheric and thermospheric arrivals detected at $x_{1}=230$  km are displayed in figures 13(a) and 13(c), respectively. The corresponding energy spectral densities are reported in figures 13(b) and 13(d). Both signals exhibit U-shaped waveforms. Their durations are around $35$ and $90$  s, and their central frequencies are equal to about $0.06$ and $0.01$  Hz.

Figure 13.

Zooms of the pressure signals $p^{\prime }$ recorded on the $x_{1}$ axis at $x_{1}=230$  km. (a) Stratospheric and (c) thermospheric arrivals and (b,d) corresponding spectra.

Figure 14.

(a) Pressure $p_{1}^{\prime }$ recorded at ground level as a function of the distance $x_{1}$ and of the retarded time $t_{r}$ . Yellow dots are associated with the couples $(x_{1},t_{r})$ for which an acoustic ray arrives at the ground. In order to visualize all the infrasonic arrivals, the colour level ranges from $-10$ to $+10$  Pa. (b) Maximum in time of $p_{1}^{\prime }(x_{1},t)$ as a function of the distance $x_{1}$ . The dashed line represents the decay of a spherical wave in a free homogeneous medium. (c) Spectrum ${\mathcal{E}}_{p^{\prime },1}(x_{1},f)$ of the arrivals recorded at ground level as a function of the distance $x_{1}$ and of the frequency $f$ .

In order to characterize the different arrivals, the pressure perturbation $p_{1}^{\prime }(x_{1},t)=p^{\prime }(x_{1},0,0,t)$ recorded on the $x_{1}$ axis is plotted in figure 14(a) as a function of the distance $x_{1}$ and of the retarded time $t_{r}=t-x_{1}/\bar{c}(0)$ . This diagram allows comparison of the arrival times of the various phases with the one of a wave propagating at ground level with a speed equal to $\bar{c}(0)$ . Yellow dots are associated with the couples $(x_{1},t_{r})$ for which an acoustic ray is detected. The maximum in time of the function $p_{1}^{\prime }$ is illustrated in figure 14(b). The direct phase $I_{dtr}$ is recorded without delay, meaning that its propagation speed is close to $\bar{c}(0)$ . Its maximum overpressure decreases with range, more rapidly than the amplitude of a spherical wave in a homogeneous medium. As discussed in § 4.5, this decay is associated with the negative gradient of the speed of sound in the troposphere, which deviates the acoustic energy towards the upper atmosphere. The direct phase $I_{dtr}$ represents the dominant contribution up to $x_{1}=130$  km. The partial reflections $I_{ps}$ are detected all along the $x_{1}$ axis, with a delay which diminishes with the distance $x_{1}$ . The amplitude of $I_{ps}$ is of the order of $1$  Pa between $x_{1}=0$  km and $x_{1}=70$  km, reaches a local minimum at $x_{1}=70$  km and then increases again up to values of about $15$  Pa at $x_{1}=140$  km. In addition, a phase jump of $\unicode[STIX]{x03C0}$ is observed at $x_{1}=70$  km. The stratospheric arrival $I_{s}$ is recorded for $x_{1}\gtrsim 100$  km, whereas, according to ray theory, it can be detected only for $x_{1}\gtrsim 212$  km. The stratospheric phase represents the most important contribution to the pressure signal for $x_{1}\gtrsim 120$  km. Its amplitude remains higher than $10$  Pa up to at least $x_{1}=450$  km, with two local maxima of $75$ and $65$  Pa at $x_{1}=230$  km and $x_{1}=275$  km, respectively. The thermospheric phase $I_{t}$ is detected at ground level for $x_{1}\gtrsim 180$  km, while the geometrical-acoustics approximation predicts no arrival for $x_{1}\lesssim 280$  km. The amplitude of the signal $I_{t}$ is much lower than the maximum overpressure of the stratospheric wave packet, mainly as a result of the nonlinear lengthening in the upper atmosphere.

The spectrum ${\mathcal{E}}_{p^{\prime },1}(x_{1},f)={\mathcal{E}}_{p^{\prime }}(x_{1},0,0,f)$ of the various arrivals recorded at ground level is finally displayed in figure 14(c) as a function of the distance along the $x_{1}$  axis and of the frequency $f$ . Up to $x_{1}=100$  km, the main contribution is due to the direct phase $I_{dtr}$ , which has an N-waveform, leading to a broadband spectrum ${\mathcal{E}}_{p^{\prime },1}$ . On the contrary, the energy spectral density ${\mathcal{E}}_{p^{\prime },1}$ of the arrivals beyond $x_{1}=100$  km is essentially contained in the frequency range $[0,0.3]\,\text{Hz}$ .

4.5 Analysis of the direct arrival $I_{dtr}$

As previously mentioned, the direct arrival $I_{dtr}$ is not predicted by the geometrical acoustics. Because of the negative vertical gradient of the speed of sound, rays are deviated towards the upper atmosphere and do not reach the ground up to $x_{1}=212$  km distance from the source location. From a physical point of view, the effect of an inhomogeneous atmosphere can be better understood with the help of an analogy proposed by Pierce (Reference Pierce1985) and Berry & Daigle (Reference Berry and Daigle1988) among others. The sound field above a flat surface in a quiescent and upward refracting atmosphere, i.e. with a negative vertical gradient of the speed of sound, can be considered equivalent to the acoustic field above a curved convex surface in a homogeneous medium, as illustrated in figure 15. In the latter case, the rays emanating from the source are straight lines and, for geometrical reasons, cannot penetrate into the region below the tangent to the curved ground, which is called shadow zone. Hence, in the shadow side of the sound field, no ray can be detected and, according to the geometrical theory, the amplitude of the pressure perturbation $p^{\prime }$ should be nil. The direct arrival $I_{dtr}$ is thus induced by diffraction and is interpreted as a creeping wave propagating along the $x_{1}$  axis.

Figure 15.

Equivalence between (a) the sound field above a flat ground in an upward refracting atmosphere and (b) the acoustic field above a convex surface in a homogeneous medium.

In order to gain more insight into how acoustic energy penetrates in the shadow zone and, in particular, on the Earth’s surface, an analysis of the phase $I_{dtr}$ is carried out in this paragraph. As formerly stated, propagation at ground level is essentially linear, excepting in the source region. Furthermore, viscous and thermal conduction effects play a minor role in the troposphere. Besides, to simplify the forthcoming developments, the mean atmosphere is modelled as an effective quiescent medium with the speed of sound and the density given by $\bar{c}_{eff}=\bar{c}+\bar{u}_{1}$ and $\bar{\unicode[STIX]{x1D70C}}_{eff}=\bar{\unicode[STIX]{x1D70C}}\bar{c}^{2}/\bar{c}_{eff}^{2}$ , respectively. As demonstrated by Godin (Reference Godin2002), such an approximation is valid for low angles of propagation with respect to the horizontal axis, and is consequently justified for a creeping wave. Under the above hypotheses, the acoustic field is axisymmetric and the normalized pressure perturbation $\unicode[STIX]{x1D6F7}=p^{\prime }/\sqrt{\bar{\unicode[STIX]{x1D70C}}_{eff}}$ can be described by the wave equation (Bergmann Reference Bergmann1946; Pierce Reference Pierce1990)

(4.4) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D6F7}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D71A}^{2}}+\frac{1}{\unicode[STIX]{x1D71A}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D71A}}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D6F7}}{\unicode[STIX]{x2202}x_{3}^{2}}+\left[\frac{1}{2\bar{\unicode[STIX]{x1D70C}}_{eff}}\frac{\text{d}^{2}\bar{\unicode[STIX]{x1D70C}}_{eff}}{\text{d}x_{3}^{2}}-\frac{3}{4\bar{\unicode[STIX]{x1D70C}}_{eff}^{2}}\left(\frac{\text{d}\bar{\unicode[STIX]{x1D70C}}_{eff}}{\text{d}x_{3}}\right)^{2}\right]\unicode[STIX]{x1D6F7}-\frac{1}{\bar{c}_{eff}^{2}(x_{3})}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D6F7}}{\unicode[STIX]{x2202}t^{2}}=\frac{G}{\sqrt{\bar{\unicode[STIX]{x1D70C}}_{eff}}},\end{eqnarray}$$

where $\unicode[STIX]{x1D71A}=\sqrt{x_{1}^{2}+x_{2}^{2}}$ and $G$ is a function of the source terms $\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D70C}}$ and $\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D70C}e_{t}}$ in (2.3). Since the characteristic acoustic wavelength $\unicode[STIX]{x1D706}\simeq 2.6$  km is smaller than the typical length scale for variations of the atmospheric density ( $\simeq 8.4$  km), the vertical derivatives of the function $\bar{\unicode[STIX]{x1D70C}}_{eff}$ are henceforth neglected. Taking the Fourier–Bessel transform

(4.5) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D712}}(k_{\unicode[STIX]{x1D71A}},x_{3},\unicode[STIX]{x1D714})=\int _{-\infty }^{+\infty }\int _{0}^{+\infty }\unicode[STIX]{x1D712}(\unicode[STIX]{x1D71A},x_{3},t){\mathcal{J}}_{0}(k_{\unicode[STIX]{x1D71A}}\unicode[STIX]{x1D71A})\text{e}^{+\text{i}\unicode[STIX]{x1D714}t}\unicode[STIX]{x1D71A}\,\text{d}\unicode[STIX]{x1D71A}\,\text{d}t,\end{eqnarray}$$

with $\unicode[STIX]{x1D712}=\unicode[STIX]{x1D6F7}$ or $\unicode[STIX]{x1D712}=G$ and where ${\mathcal{J}}_{0}$ is the zeroth-order Bessel function of the first kind, the depth-separated Helmholtz equation is obtained:

(4.6) $$\begin{eqnarray}\frac{\text{d}^{2}\hat{\unicode[STIX]{x1D6F7}}}{\text{d}x_{3}^{2}}+\left(\frac{\unicode[STIX]{x1D714}^{2}}{\bar{c}_{eff}}-k_{\unicode[STIX]{x1D71A}}^{2}\right)\hat{\unicode[STIX]{x1D6F7}}=\frac{{\hat{G}}}{\sqrt{\bar{\unicode[STIX]{x1D70C}}_{eff}}}.\end{eqnarray}$$

The normalized pressure perturbation $\unicode[STIX]{x1D6F7}$ is then retrieved by the inverse transform

(4.7) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}(\unicode[STIX]{x1D71A},x_{3},t)=\frac{1}{4\unicode[STIX]{x03C0}}\iint _{-\infty }^{+\infty }\unicode[STIX]{x1D6F7}(k_{\unicode[STIX]{x1D71A}},x_{3},\unicode[STIX]{x1D714}){\mathcal{H}}_{0}^{(1)}(k_{\unicode[STIX]{x1D71A}}\unicode[STIX]{x1D71A})\text{e}^{-\text{i}\unicode[STIX]{x1D714}t}k_{\unicode[STIX]{x1D71A}}\,\text{d}k_{\unicode[STIX]{x1D71A}}\,\text{d}\unicode[STIX]{x1D714},\end{eqnarray}$$

where ${\mathcal{H}}_{0}^{(1)}$ is the zeroth-order Hankel function of the first kind. By means of the residue theorem, it is possible to demonstrate that integral (4.7) can be approximated by a linear combination of the eigenfunctions $\hat{\unicode[STIX]{x1D6F7}}$ of the differential equation (Jensen et al. Reference Jensen, Kuperman, Porter and Schmidt2011)

(4.8a ) $$\begin{eqnarray}\frac{\text{d}^{2}\hat{\unicode[STIX]{x1D6F7}}}{\text{d}x_{3}^{2}}+\left(\frac{\unicode[STIX]{x1D714}^{2}}{\bar{c}_{eff}}-k_{\unicode[STIX]{x1D71A}}^{2}\right)\hat{\unicode[STIX]{x1D6F7}}=0,\end{eqnarray}$$

completed with the rigid-wall constraint

(4.8b ) $$\begin{eqnarray}\frac{\text{d}\hat{\unicode[STIX]{x1D6F7}}}{\text{d}x_{3}}=0,\quad x_{3}=0~\text{km},\end{eqnarray}$$

and with the radiation condition

(4.8c ) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D6F7}}\sim \text{e}^{\text{i}k_{3}^{\infty }x_{3}},\quad k_{3}^{\infty }=\lim _{x_{3}\rightarrow +\infty }\sqrt{\frac{\unicode[STIX]{x1D714}^{2}}{\bar{c}_{eff}(x_{3})}-k_{\unicode[STIX]{x1D71A}}^{2}},\quad x_{3}\rightarrow +\infty .\end{eqnarray}$$

The real part of the wavenumber $k_{3}^{\infty }$ is required to be positive or nil in order for the function $\hat{\unicode[STIX]{x1D6F7}}$ to represent an outgoing wave as $x_{3}\rightarrow +\infty$ . Since the diffracted arrival $I_{dtr}$ is mainly influenced by the negative gradient of the temperature in the troposphere, the function $\bar{c}_{eff}$ is here defined as

(4.8d ) $$\begin{eqnarray}\bar{c}_{eff}(x_{3})=\left\{\begin{array}{@{}ll@{}}\bar{c}(x_{3})+\bar{u}_{1}(x_{3}),\quad & x_{3}\in [0,x_{3}^{\infty }]\\ \bar{c}(x_{3}^{\infty })+\bar{u}_{1}(x_{3}^{\infty }),\quad & x_{3}>x_{3}^{\infty },\end{array}\right.\end{eqnarray}$$

Figure 16.

(a) Effective speed of sound profile $\bar{c}_{eff}$ . (b) Wavenumbers $k_{\unicode[STIX]{x1D71A}}$ for $f=0.2$  Hz; the red square corresponds to the least attenuated mode.

where $x_{3}^{\infty }=15.75$  km is the altitude of the first minimum of the sum $\bar{c}(x_{3})+\bar{u}_{1}(x_{3})$ . The resulting profile $\bar{c}_{eff}(x_{3})$ is displayed in figure 16(a). In line with this definition, the above radiation condition is replaced by the following boundary constraint (Jensen et al. Reference Jensen, Kuperman, Porter and Schmidt2011):

(4.8e ) $$\begin{eqnarray}\frac{\text{d}\hat{\unicode[STIX]{x1D6F7}}}{\text{d}x_{3}}-\text{i}k_{3}^{\infty }\hat{\unicode[STIX]{x1D6F7}}=0,\quad k_{3}^{\infty }=\sqrt{\frac{\unicode[STIX]{x1D714}^{2}}{\bar{c}_{eff}(x_{3}^{\infty })}-k_{\unicode[STIX]{x1D71A}}^{2}},\quad x_{3}=x_{3}^{\infty }.\end{eqnarray}$$

Problem (4.8) represents a differential eigensystem and can be solved through a pseudospectral collocation method (Bertin et al. Reference Bertin, Millet and Bouche2014; Sabatini & Bailly Reference Sabatini and Bailly2015). As an illustration, the eigenvalues $k_{\unicode[STIX]{x1D71A}}$ obtained for $f=0.2$  Hz are plotted in figure 16(b). As suggested by Pierce (Reference Pierce1985), along the $x_{1}$ axis, integral (4.7) is dominated by wavenumbers $k_{\unicode[STIX]{x1D71A}}$ close to $k_{0}=\unicode[STIX]{x1D714}/\bar{c}_{eff}(0)$ . Moreover, since the imaginary parts of the eigenvalues $k_{\unicode[STIX]{x1D71A}}$ are strictly positive, the corresponding contributions, which are proportional to ${\mathcal{H}}_{0}^{(1)}(k_{\unicode[STIX]{x1D71A}}r)$ , decay with distance. Consequently, at large distances, only the least attenuated mode is likely to be observed. The phase speed $v_{\unicode[STIX]{x1D719}}=\unicode[STIX]{x1D714}/\text{Re}(k_{\unicode[STIX]{x1D71A}})$ and the imaginary part $\unicode[STIX]{x1D6FC}=\min _{k_{\unicode[STIX]{x1D71A}}}(\text{Im}(k_{\unicode[STIX]{x1D71A}}))$ of the least attenuated mode are depicted as functions of the frequency $f$ in figure 17(a). Globally, they both depend weakly on the pulsation $\unicode[STIX]{x1D714}$ . For $f>0.05$  Hz, the velocity $v_{\unicode[STIX]{x1D719}}$ is approximately equal to $\bar{c}_{eff}(0)=340~\text{m}~\text{s}^{-1}$ , and the attenuation factor $\unicode[STIX]{x1D6FC}$ increases only slightly as $f$ augments and has an average value of about $\unicode[STIX]{x1D6FC}_{mean}=2.85\times 10^{-2}~\text{km}^{-1}$ in the range $[0.05,0.2]$  Hz. As the wind velocity is nil at ground level, the least attenuated mode propagates with a celerity equal to the speed of sound on the Earth’s surface, $\bar{c}(0)$ , which confirms the result of the previous paragraph. Moreover, by virtue of the above considerations, the absolute value of the pressure perturbation $p^{\prime }$ along the $x_{1}$ axis and in the far range is found to behave as

(4.9) $$\begin{eqnarray}p^{\prime }(x_{1},0,0,t)\propto \frac{1}{\sqrt{x_{1}}}\text{e}^{-\unicode[STIX]{x1D6FC}_{mean}x_{1}},\end{eqnarray}$$

since ${\mathcal{H}}_{0}^{(1)}(k_{\unicode[STIX]{x1D71A}}x_{1})\sim 1/\sqrt{x_{1}}$ for $k_{\unicode[STIX]{x1D71A}}x_{1}\gg 1$ . Hence, conforming to earlier remarks, the amplitude of the phase $I_{dtr}$ weakens more rapidly than that of a spherical wave in a homogeneous medium, for which $p^{\prime }(x_{1},0,0,t)\propto 1/x_{1}$ . To conclude, the maximum in time of the function $p_{1}^{\prime }$ is again illustrated in figure 17(b) and compared to the result given by expression (4.9). A very good agreement is observed between numerical and theoretical predictions for $x_{1}>20$  km.

Figure 17.

(a) Phase speed (black line) and imaginary part (blue line) of the least attenuated mode. (b) Maximum in time of the function $p_{1}^{\prime }(x_{1},t)$ as a function of the distance $x_{1}$ , compared to the result given by formula (4.9) and to the behaviour of a spherical wave in a homogeneous medium.

4.6 Analysis of the partial reflections

The sensitivity of ground recordings to the speed of sound and wind fluctuations of length scale of the same order as the acoustic wavelength has received particular attention in recent years (Chunchuzov et al. Reference Chunchuzov, Kulichkov, Popov, Waxler and Assink2011, Reference Chunchuzov, Kulichkov, Popov and Hedlin2014, Reference Chunchuzov, Kulichkov, Popov, Perepelkin, Vasilev, Glushkov and Firstov2015; Bertin et al. Reference Bertin, Millet and Bouche2014). In particular, Chunchuzov et al. (Reference Chunchuzov, Kulichkov, Popov and Hedlin2014, Reference Chunchuzov, Kulichkov, Popov, Perepelkin, Vasilev, Glushkov and Firstov2015) demonstrated the feasibility of retrieving the vertical structure of the atmosphere from the reflected signals, whose energy spectral densities are closely related to the characteristic wavenumbers of the inhomogeneities. In order to gain more insight into the relationship between the wavenumber spectrum of the atmospheric fluctuations and the Fourier transform of the partial reflections, a more detailed analysis of the phase $I_{ps}$ is now carried out on a simplified model.

4.6.1 Reflection of a spherical wave by a one-dimensional speed-of-sound inhomogeneity

The configuration schematically illustrated in figure 18 is considered. A spherical wavefront is emitted in an infinite domain $Ox_{1}x_{2}x_{3}$ by a point source placed at the origin $O$ and propagates through a one-dimensional speed-of-sound inhomogeneity confined between the altitudes $x_{3wb}$ and $x_{3wt}$ ; the medium is at rest and the speed of sound away from the inhomogeneity is equal to $\bar{c}_{\infty }$ . The interaction between the wavefront and the speed-of-sound fluctuation generates a transmitted wave and a reflected wave. In order to simplify the analytical developments, nonlinear, viscous and thermal conduction effects are henceforth neglected. Moreover, since the large-scale variations of the atmospheric medium induced by the gravitational acceleration do not contribute to the reflected wave, the gradient $\text{d}\bar{p}/\text{d}x_{3}$ is assumed nil and the ambient pressure $\bar{p}$ is set equal to $\bar{p}_{0}$ . Under these hypotheses, the pressure perturbation $p^{\prime }$ is described by the wave equation (Bergmann Reference Bergmann1946; Pierce Reference Pierce1990)

(4.10) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}p^{\prime }}{\unicode[STIX]{x2202}x_{1}^{2}}+\frac{\unicode[STIX]{x2202}^{2}p^{\prime }}{\unicode[STIX]{x2202}x_{2}^{2}}+\bar{\unicode[STIX]{x1D70C}}(x_{3})\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{3}}\left(\frac{1}{\unicode[STIX]{x1D70C}(x_{3})}\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}x_{3}}\right)-\frac{1}{\bar{c}^{2}(x_{3})}\frac{\unicode[STIX]{x2202}^{2}p^{\prime }}{\unicode[STIX]{x2202}t^{2}}=-F(t)\unicode[STIX]{x1D6FF}(\boldsymbol{x}),\end{eqnarray}$$

Figure 18.

Propagation of a spherical wave through a small-scale speed-of-sound inhomogeneity.

where $F(t)$ indicates the temporal envelope of the point source and $\bar{\unicode[STIX]{x1D70C}}(x_{3})=\unicode[STIX]{x1D6FE}\bar{p}_{0}/\bar{c}^{2}(x_{3})$ . As shown in appendix B, for $x_{3}\ll x_{3w}$ , the incident sound field $p_{i}^{\prime }(\boldsymbol{x},t)$ is represented by

(4.11) $$\begin{eqnarray}p_{i}^{\prime }(\boldsymbol{x},t)=\frac{F(t-\Vert \boldsymbol{x}\Vert /\bar{c}_{\infty })}{4\unicode[STIX]{x03C0}\Vert \boldsymbol{x}\Vert }.\end{eqnarray}$$

Accordingly, the temporal shape of the source is advected in space at the speed $\bar{c}_{\infty }$ and the maximum overpressure decays as $1/\Vert \boldsymbol{x}\Vert$ . Moreover, the temporal Fourier transform $\hat{p}_{r}^{\prime }(\boldsymbol{x},\unicode[STIX]{x1D714})$ of the reflected wave $p_{r}^{\prime }(\boldsymbol{x},t)$ ,

(4.12) $$\begin{eqnarray}\hat{p}_{r}^{\prime }(\boldsymbol{x},\unicode[STIX]{x1D714})=\int _{-\infty }^{+\infty }p_{r}^{\prime }(\boldsymbol{x},t)\text{e}^{\text{i}\unicode[STIX]{x1D714}t}\,\text{d}t,\end{eqnarray}$$

is given by the integral

(4.13) $$\begin{eqnarray}\hat{p}_{r}^{\prime }(\boldsymbol{x},\unicode[STIX]{x1D714})=\frac{\hat{F}(\unicode[STIX]{x1D714})}{4\unicode[STIX]{x03C0}}\iint _{-\infty }^{+\infty }{\displaystyle \frac{\text{i}{\mathcal{R}}_{c}\text{e}^{+\text{i}[k_{1}x_{1}+k_{2}x_{2}-k_{3}(k_{1},k_{2},\unicode[STIX]{x1D714})x_{3}]}}{2\unicode[STIX]{x03C0}k_{3}(k_{1},k_{2},\unicode[STIX]{x1D714})}}\,\text{d}k_{1}\,\text{d}k_{2},\end{eqnarray}$$

where $k_{1}$ , $k_{2}$ and $k_{3}=\sqrt{(\unicode[STIX]{x1D714}/\bar{c}_{\infty })^{2}-(k_{1}^{2}+k_{2}^{2})}$ represent the wavenumber components, $\hat{F}$ is the Fourier transform of $F(t)$ and ${\mathcal{R}}_{c}$ is the reflection coefficient of the inhomogeneous layer. The acoustic field $\hat{p}_{r}^{\prime }$ can be interpreted as a sum of monochromatic plane waves whose amplitude is proportional to both the incident spectrum $\hat{F}$ and the reflection coefficient ${\mathcal{R}}_{c}$ . The value of ${\mathcal{R}}_{c}$ depends on the speed of sound $\bar{c}(x_{3})$ , on the incidence angle $\unicode[STIX]{x1D711}=\text{atan}(k_{3}/\sqrt{k_{1}^{2}+k_{2}^{2}})$ with respect to the horizontal plane and on the pulsation $\unicode[STIX]{x1D714}$ . According to Lekner (Reference Lekner1987), the reflection coefficient ${\mathcal{R}}_{c}$ is defined as the limit

(4.14) $$\begin{eqnarray}{\mathcal{R}}_{c}\equiv \lim _{x_{3}\rightarrow -\infty }r_{c}(x_{3})\text{e}^{2\text{i}\tilde{k}_{3}(x_{3})x_{3}},\end{eqnarray}$$

where the wavenumber $\tilde{k}_{3}$ is given by

(4.15) $$\begin{eqnarray}\tilde{k}_{3}^{2}(x_{3})\equiv \frac{\unicode[STIX]{x1D714}^{2}}{\bar{c}^{2}(x_{3})}-k_{1}^{2}-k_{2}^{2}=\frac{\unicode[STIX]{x1D714}^{2}}{\bar{c}^{2}(x_{3})}-\frac{\unicode[STIX]{x1D714}^{2}}{\bar{c}_{\infty }^{2}}\cos ^{2}(\unicode[STIX]{x1D711}),\end{eqnarray}$$

and the function $r_{c}(x_{3})$ satisfies the Riccati problem

(4.16) $$\begin{eqnarray}\frac{\text{d}r_{c}}{\text{d}x_{3}}+2\text{i}\tilde{k}_{3}r_{c}-\frac{\bar{\unicode[STIX]{x1D70C}}}{2\tilde{k_{3}}}\frac{\text{d}(\tilde{k}_{3}/\bar{\unicode[STIX]{x1D70C}})}{\text{d}x_{3}}(1-r_{c}^{2})=0,\end{eqnarray}$$

with the condition

(4.17) $$\begin{eqnarray}\lim _{x_{3}\rightarrow +\infty }r_{c}=0.\end{eqnarray}$$

Since the inhomogeneous layer is assumed confined in the interval $[x_{3wb},x_{3wt}]$ , $r_{c}(x_{3}^{\infty })=0$ for $x_{3}^{\infty }>x_{3wt}$ and the above equation can be integrated from $x_{3}=x_{3}^{\infty }$ to $x_{3}=x_{3}^{-\infty }<x_{3wb}$ through any stable time-integration algorithm.

4.6.2 Reflection coefficient of a Morlet wavelet

The reflection coefficient is now analysed for a speed of sound $\bar{c}(x_{3})$ defined as the sum of a constant $\bar{c}_{\infty }$ and an elementary Morlet wavelet of amplitude $\unicode[STIX]{x1D716}_{\bar{c}}\bar{c}_{\infty }$ , of wavelength $\unicode[STIX]{x1D706}_{w}$ and located at altitude $x_{3w}$ :

(4.18) $$\begin{eqnarray}\bar{c}(x_{3})=\bar{c}_{\infty }[1+\unicode[STIX]{x1D716}_{\bar{c}}\text{e}^{-z^{2}/2}\cos (2\unicode[STIX]{x03C0}z)],\quad \text{with}~z=\frac{x_{3}-x_{3w}}{\unicode[STIX]{x1D706}_{w}}.\end{eqnarray}$$

Hereafter, the parameter $\unicode[STIX]{x1D716}_{\bar{c}}$ will be considered small. The modulus $|{\mathcal{R}}_{c}|$ obtained with $\unicode[STIX]{x1D716}_{\bar{c}}\bar{c}_{\infty }=5~\text{m}~\text{s}^{-1}$ is plotted in figure 19 as a function of the normalized pulsation $\unicode[STIX]{x1D714}\unicode[STIX]{x1D706}_{w}/\bar{c}_{\infty }$ and of the sinus of the incidence angle $\unicode[STIX]{x1D711}$ . The function $|{\mathcal{R}}_{c}|$ is different from zero only in specific frequency bands and its maxima approximately satisfy the well-known Bragg law (Blanc-Benon, Dallois & Juvé Reference Blanc-Benon, Dallois and Juvé2001)

(4.19) $$\begin{eqnarray}2\unicode[STIX]{x1D706}_{w}\sin (\unicode[STIX]{x1D711})=l\unicode[STIX]{x1D706},\end{eqnarray}$$

where $l$ is an integer number and $\unicode[STIX]{x1D706}=2\unicode[STIX]{x03C0}\bar{c}_{\infty }/\unicode[STIX]{x1D714}$ is the incident wavelength. The modulus $|{\mathcal{R}}_{c}|$ increases as $\unicode[STIX]{x1D711}$ tends to $0$ , and is close to $1$ for incidences $\unicode[STIX]{x1D711}$ lower than about $10^{\circ }$ . In this case, the associated plane wave is totally reflected. Moreover, in each frequency band, there exists an angle $\unicode[STIX]{x1D711}_{t}$ , called intromission angle, for which ${\mathcal{R}}_{c}$ is nil and the signal is totally transmitted. As an example, $\unicode[STIX]{x1D711}_{t}$ is equal to $45^{\circ }$ for $l=1$ and to about $57.3^{\circ }$ for $l=2$ .

Figure 19.

Reflection coefficient ${\mathcal{R}}_{c}$ for the speed of sound profile (4.18) as a function of the normalized pulsation $\unicode[STIX]{x1D714}\unicode[STIX]{x1D706}_{w}\bar{c}_{\infty }$ and of the sinus of the incidence angle $\sin (\unicode[STIX]{x1D711})$ . The dashed lines represent the well-known Bragg law $2\unicode[STIX]{x1D706}_{w}\sin (\unicode[STIX]{x1D711})=l\unicode[STIX]{x1D706}$ , for $l=1$ , $l=2$ and $l=3$ .

According to Lekner (Reference Lekner1987), under the hypothesis of weak reflection, a first approximation of the function ${\mathcal{R}}_{c}$ for a generic inhomogeneity $\bar{c}^{\prime }=\bar{c}-\bar{c}_{\infty }$ is given by the integral

(4.20) $$\begin{eqnarray}{\mathcal{R}}_{c}\simeq \int _{-\infty }^{+\infty }\left[\frac{\text{i}\unicode[STIX]{x0394}\tilde{k}_{3}^{2}}{2k_{3}}+\frac{1}{2}\frac{1}{\bar{\unicode[STIX]{x1D70C}}}\frac{\text{d}\bar{\unicode[STIX]{x1D70C}}}{\text{d}x_{3}}\right]\text{e}^{2\text{i}k_{3}x_{3}}\,\text{d}x_{3},\quad \text{with}~\unicode[STIX]{x0394}\tilde{k}_{3}^{2}(x_{3})=\frac{\unicode[STIX]{x1D714}^{2}}{\bar{c}^{2}(x_{3})}-\frac{\unicode[STIX]{x1D714}^{2}}{\bar{c}_{\infty }^{2}}.\end{eqnarray}$$

For weak inhomogeneities ( $\bar{c}^{\prime }/\bar{c}_{\infty }\ll 1$ ), the density gradient can be rewritten as $(\text{d}\bar{\unicode[STIX]{x1D70C}}/\text{d}x_{3})/\bar{\unicode[STIX]{x1D70C}}=-2(\text{d}\bar{c}^{\prime }/\text{d}x_{3})/\bar{c}_{\infty }$ . As a result, integral (4.20) can be further simplified as

(4.21) $$\begin{eqnarray}{\mathcal{R}}_{c}\simeq \text{i}\frac{1}{k_{3}}\frac{\unicode[STIX]{x1D714}^{2}}{\bar{c}_{\infty }^{2}}(2\sin ^{2}(\unicode[STIX]{x1D711})-1)\int _{-\infty }^{+\infty }\frac{\bar{c}^{\prime }(x_{3})}{\bar{c}_{\infty }}\text{e}^{2\text{i}k_{3}x_{3}}\,\text{d}x_{3}.\end{eqnarray}$$

This formula allows relating the reflection coefficient ${\mathcal{R}}_{c}$ to the spectrum of the speed of sound fluctuations; ${\mathcal{R}}_{c}$ is indeed proportional to the Fourier transform of the inhomogeneity $\bar{c}^{\prime }/\bar{c}_{\infty }$ . Furthermore, because of the factor $(2\sin ^{2}(\unicode[STIX]{x1D711})-1)$ , the function ${\mathcal{R}}_{c}$ vanishes for $\unicode[STIX]{x1D711}=\unicode[STIX]{x1D711}_{t}=45^{\circ }$ . Hence, the existence of an intromission angle $\unicode[STIX]{x1D711}_{t}$ is not peculiar to the profile (4.18) and is strictly linked to the density variations $(\text{d}\bar{\unicode[STIX]{x1D70C}}/\text{d}x_{3})/\bar{\unicode[STIX]{x1D70C}}$ , which, under the assumption $\text{d}\bar{p}/\text{d}x_{3}=0$ , are only induced by the gradient of the speed of sound $\text{d}\bar{c}^{\prime }/\text{d}x_{3}$ . For $\unicode[STIX]{x1D711}=\unicode[STIX]{x1D711}_{t}$ , the effect of the fluctuations of $\bar{\unicode[STIX]{x1D70C}}$ is exactly compensated by the effect of the speed-of-sound inhomogeneity.

To conclude, for the specific case of the profile (4.18), it can be shown that the reflection coefficient ${\mathcal{R}}_{c}$ can be written as ${\mathcal{R}}_{c}\equiv {\mathcal{R}}\text{e}^{2\text{i}k_{3}\unicode[STIX]{x1D706}_{w}}$ , where the auxiliary function ${\mathcal{R}}$ is approximately equal to

(4.22) $$\begin{eqnarray}{\mathcal{R}}\simeq \text{i}\frac{\sqrt{2\unicode[STIX]{x03C0}}\,\unicode[STIX]{x1D716}_{\bar{c}}}{k_{3}\unicode[STIX]{x1D706}_{w}}\frac{\unicode[STIX]{x1D714}^{2}\unicode[STIX]{x1D706}_{w}^{2}}{\bar{c}_{\infty }^{2}}(2\sin ^{2}(\unicode[STIX]{x1D711})-1)\text{e}^{-2k_{3}^{2}\unicode[STIX]{x1D706}_{w}^{2}-2\unicode[STIX]{x03C0}^{2}}\cosh (4\unicode[STIX]{x03C0}k_{3}\unicode[STIX]{x1D706}_{w}).\end{eqnarray}$$

This function takes pure imaginary values and exhibits a jump of $\unicode[STIX]{x03C0}$ at the intromission angle $\unicode[STIX]{x1D711}=\unicode[STIX]{x1D711}_{t}$ . Moreover, the modulus $|{\mathcal{R}}_{c}|=|{\mathcal{R}}|$ increases linearly with the amplitude of the speed-of-sound inhomogeneity.

4.6.3 Amplitude of the wave reflected by a Morlet wavelet along the $x_{1}$ axis

By using the auxiliary variable ${\mathcal{R}}$ , the reflected sound field (4.13) can be recast as follows:

(4.23) $$\begin{eqnarray}\hat{p}_{r}^{\prime }(\boldsymbol{x},\unicode[STIX]{x1D714})=\frac{\hat{F}(\unicode[STIX]{x1D714})}{4\unicode[STIX]{x03C0}}\iint _{-\infty }^{+\infty }{\displaystyle \frac{\text{i}{\mathcal{R}}\text{e}^{+\text{i}[k_{1}x_{1}+k_{2}x_{2}-k_{3}(k_{1},k_{2},\unicode[STIX]{x1D714})(x_{3}-2x_{3w})]}}{2\unicode[STIX]{x03C0}k_{3}(k_{1},k_{2},\unicode[STIX]{x1D714})}}\,\text{d}k_{1}\,\text{d}k_{2}.\end{eqnarray}$$

This integral can be drastically simplified through the method of the stationary phase (Frisk Reference Frisk1994). More specifically, with reference to figure 20, it can be shown that, at a given recording station $S_{R}$ , the function $\hat{p}_{r}^{\prime }(\boldsymbol{x},\unicode[STIX]{x1D714})$ can be approximated as

(4.24) $$\begin{eqnarray}\hat{p}_{r}^{\prime }(\boldsymbol{x},\unicode[STIX]{x1D714})\simeq \hat{F}(\unicode[STIX]{x1D714}){\mathcal{R}}(\unicode[STIX]{x1D711}_{st},\unicode[STIX]{x1D714})\frac{\text{e}^{-\text{i}\unicode[STIX]{x1D714}r_{\star }/\bar{c}_{\infty }}}{4\unicode[STIX]{x03C0}r_{\star }},\end{eqnarray}$$

where $\unicode[STIX]{x1D711}_{st}=\text{atan}(x_{3w}/(\sqrt{x_{1}^{2}+x_{2}^{2}}/2))$ and $r_{\star }=\Vert \overrightarrow{O^{\prime }S_{R}}\Vert$ . The variable $\unicode[STIX]{x1D711}_{st}$ represents the angle formed with the horizontal plane by the line connecting the source $O$ with the intersection point between the incident wavefront and the symmetry axis $x_{3}=x_{3w}$ of the speed-of-sound profile; the term $r_{\star }$ is the distance between the point $O^{\prime }$ , located at the altitude $x_{3}=2x_{3w}$ , and the recording station $S_{R}$ . From a physical point of view, equation (4.24) states that the reflected field $\hat{p}_{r}^{\prime }$ consists of a spherical wave emanating from the image source $O^{\prime }$ , whose amplitude is proportional to the incident spectrum $\hat{F}(\unicode[STIX]{x1D714})$ and to the reflection coefficient ${\mathcal{R}}$ for a plane wave having an incidence angle equal to the specular angle $\unicode[STIX]{x1D711}=\unicode[STIX]{x1D711}_{st}$ . Expression (4.24) represents the geometrical approximation of the reflected field, according to which the reflected energy is mainly concentrated in the direction of the specular angle. The modulus $|\hat{p}_{r}^{\prime }|$ is finally given by the relation

(4.25) $$\begin{eqnarray}|\hat{p}_{r}^{\prime }(\boldsymbol{x},\unicode[STIX]{x1D714})|\simeq \left|\frac{\hat{F}(\unicode[STIX]{x1D714}){\mathcal{R}}(\unicode[STIX]{x1D711}_{st},\unicode[STIX]{x1D714})}{4\unicode[STIX]{x03C0}r_{\star }}\right|.\end{eqnarray}$$

Since the reflection coefficient $|{\mathcal{R}}|$ is different from zero only in particular frequency bands, the above formula states that the inhomogeneous layer acts as a selective filter which reflects only specific wavelengths of the incident spectrum.

Figure 20.

Propagation of a spherical wave through a small-scale inhomogeneity. At the receiver $S_{R}$ , the main contribution to the reflected wave comes from the specular angle $\unicode[STIX]{x1D711}_{st}$ , formed with the $x_{1}$ axis by the ray connecting the point source $O$ and the intersection between the incident spherical wave and the symmetry axis of the inhomogeneity.

4.6.5 Nonlinear effects on the partial reflections

As a result of the nonlinear steepening and lengthening, the spectrum of the incident wave $p_{i}^{\prime }$ evolves while the wavefront travels from the source location to the small-scale inhomogeneity. Oppositely, the reflected wave $p_{r}^{\prime }$ has generally a small amplitude, so that its propagation towards the Earth’s surface can be considered linear. Accordingly, supposing that the Fourier transform of the incident signal does not change significantly as this latter crosses the small-scale fluctuation, formula (4.25) remains valid as long as the function $\hat{F}(\unicode[STIX]{x1D714})$ is interpreted as the incident spectrum at the inhomogeneity location.

Figure 21.

(a) Wavefronts obtained by ray theory at two different instants of time. (b) Angle $\unicode[STIX]{x1D711}_{st}$ between the geometrical wavefront and the reference axis $x_{3}=40$  km as a function of the distance $x_{1}$ of the receiver. (c) Spectrum of the arrival $I_{ps}$ as a function of the distance $x_{1}$ and of the frequency $f$ .

4.7 Diffraction at the thermospheric caustic

As previously mentioned, diffraction at caustics plays a significant role in atmospheric infrasound propagation. In order to demonstrate the importance of this phenomenon, a more detailed study of the thermospheric caustic is finally carried out.

A caustic is defined within the framework of geometrical acoustics and corresponds to the locus of points where the ray-tube area vanishes, or, equivalently, to the locus of points where two infinitely adjacent rays intersect. The projection on the $x_{1}{-}x_{3}$ plane of the thermospheric caustic associated with the profiles of the speed of sound and the wind employed in the present simulation is plotted in grey in figure 22. The normalized pressure field $\unicode[STIX]{x1D6F7}$ obtained at $t_{5}$ is displayed as well, along with the wavefront computed by ray theory, which is represented by the green and black curves. The caustic divides the acoustic field into two regions, an ensonified zone between the grey lines, also called a lit region, and a shadow zone. In the context of the geometrical-acoustics approximation, the amplitude of the pressure field should be different from zero only in the lit region. However, as shown in figure 22, part of the acoustic energy penetrates in the shadow zone by diffraction. For this reason, the thermospheric signal $I_{t}$ is first recorded at about $x_{1}=180$  km on the Earth’s surface (red point in figure 22), whereas, according to ray theory, no thermospheric arrival should be detected between the source location and $x_{1}=280$  km (black point in figure 22) corresponding to the intersection between the caustic and the axis $x_{1}$ .

Figure 22.

Normalized pressure field $\unicode[STIX]{x1D6F7}$ obtained on the $x_{1}{-}x_{3}$ plane at $t=775.89$  s. The grey line represents the thermospheric caustic, whereas the green and black lines constitute the wavefront obtained by ray theory.

The energy spectral density ${\mathcal{E}}_{\unicode[STIX]{x1D6F7},\unicode[STIX]{x1D702}}$ of the signals recorded along the axis $\unicode[STIX]{x1D702}$ , perpendicular to the caustic at the point $P_{c}$ (cf. figure 22), is illustrated in figure 23(a) as a function of the distance $\unicode[STIX]{x1D702}$ and of the frequency $f$ . The spectrum of the diffracted field is mainly concentrated around $f=0.01$  Hz. Its behaviour along the axis $\unicode[STIX]{x1D702}$ for $f=0.01$  Hz is plotted in figure 23(b). In the shadow zone, for $\unicode[STIX]{x1D702}>0$ (at the right of the red dashed line), the energy spectral density ${\mathcal{E}}_{\unicode[STIX]{x1D6F7},\unicode[STIX]{x1D702}}$ exhibits an exponential-like decay. This result can be qualitatively explained by the linear theory of diffraction, according to which the amplitude of a spectral component of the pressure field, having frequency equal to $f$ , decreases as $\text{Ai}(\tilde{\unicode[STIX]{x1D702}})\sim \tilde{\unicode[STIX]{x1D702}}^{-1/4}\text{e}^{-2/3\tilde{\unicode[STIX]{x1D702}}^{3/2}}$ , where $\text{Ai}$ is the Airy function of the first kind and $\tilde{\unicode[STIX]{x1D702}}$ is a non-dimensional variable proportional to $\unicode[STIX]{x1D702}f^{2/3}$ .

The dependence of $\tilde{\unicode[STIX]{x1D702}}$ on the frequency $f$ also allows one to argue that the distance $x_{1}$ on the ground at which the thermospheric arrival $I_{t}$ is first detected augments with the amplitude of the source. Indeed, as the source energy increases, the lengthening of the N-wave is enhanced by the nonlinear effects in the upper atmosphere. As a result, the central frequency of the signal reaching the caustic decreases, the decay of the Airy function slows down and the penetration depth in the shadow zone grows.

To conclude, the wavefront obtained by ray theory is also represented in figure 22 by the green and black lines. Within the framework of geometrical acoustics, it can be shown that the waveform on the green side of the wavefront is given by the Hilbert transform of the waveform on the black side, at least in the linear regime (Pierce Reference Pierce1985). This result allows one to qualitatively understand the origin of the U-shaped signature of the thermospheric arrival (see figure 13 c), which corresponds to the Hilbert transform of the N-shaped waveform generated by the impulsive source.

Figure 23.

(a) Energy spectral density ${\mathcal{E}}_{\unicode[STIX]{x1D6F7},\unicode[STIX]{x1D702}}(\unicode[STIX]{x1D702},f)$ of the signals recorded along the axis $\unicode[STIX]{x1D702}$ , normal to the caustic at the point $P_{c}=(192.2~\text{km},0~\text{km},83.25~\text{km})$ . (b) Energy spectral density ${\mathcal{E}}_{\unicode[STIX]{x1D6F7},\unicode[STIX]{x1D702}}(\unicode[STIX]{x1D702},f)$ for $f=0.01$  Hz.