Finite-Difference Modeling of Acoustic and Gravity Wave Propagation in Mars Atmosphere: Application to Infrasounds Emitted by Meteor Impacts

Abstract

The propagation of acoustic and gravity waves in planetary atmospheres is strongly dependent on both wind conditions and attenuation properties. This study presents a finite-difference modeling tool tailored for acoustic-gravity wave applications that takes into account the effect of background winds, attenuation phenomena (including relaxation effects specific to carbon dioxide atmospheres) and wave amplification by exponential density decrease with height. The simulation tool is implemented in 2D Cartesian coordinates and first validated by comparison with analytical solutions for benchmark problems. It is then applied to surface explosions simulating meteor impacts on Mars in various Martian atmospheric conditions inferred from global climate models. The acoustic wave travel times are validated by comparison with 2D ray tracing in a windy atmosphere. Our simulations predict that acoustic waves generated by impacts can refract back to the surface on wind ducts at high altitude. In addition, due to the strong nighttime near-surface temperature gradient on Mars, the acoustic waves are trapped in a waveguide close to the surface, which allows a night-side detection of impacts at large distances in Mars plains. Such theoretical predictions are directly applicable to future measurements by the INSIGHT NASA Discovery mission.

Appendices

Appendix 1: Analytical Solution for Acoustic Wave Propagation in a Windy Homogeneous Atmosphere with Vibrational Relaxation

In the framework of acoustic wave propagation in a windy homogeneous medium with vibrational relaxation, the state and momentum equations read

$$ \begin{aligned} (\partial_{t} + \mathbf{w}.\nabla ) p &= -\nabla .\mathbf{v}\ast \varPhi -\rho c^{2} S, \\ \rho (\partial_{t} + \mathbf{w}.\nabla )\mathbf{v} &= -\nabla p, \end{aligned} $$

(22)

where, $\forall t \in \mathbb{R}^{+}$

$$ \varPhi (t) = H(t)\ast \psi (t), $$

(23)

with $H$ the Heaviside function and $\psi$ the relaxation function of the viscoelastic medium considered.

Taking the divergence of (22)-2 yields

$$ (\partial_{t} + \mathbf{w}.\nabla )\nabla .\mathbf{v} = - \frac{1}{ \rho }\Delta p, $$

(24)

where $\Delta$ is the Laplacian operator. Thus, taking the total derivative $D_{t}$ (where $D_{t} = \partial_{t} + \mathbf{v}\nabla$) of (22)-1 yields

$$ \begin{aligned}[b] (\partial_{t} + \mathbf{w}.\nabla )^{2} p &= -\nabla .(\partial_{t} + \mathbf{w}.\nabla )\mathbf{v}\ast \varPhi - (\partial_{t} + \mathbf{w}. \nabla )S \\ &= \frac{1}{\rho }\Delta p\ast \varPhi -\rho c^{2} (\partial_{t} + \mathbf{w}.\nabla )S. \end{aligned} $$

(25)

By considering, in (25), plane wave solutions of the form $p = \tilde{p}e^{-i\omega t}$ and $\varPsi$, defined in (16), the relaxation function for a Zener mechanism (Carcione et al. 1988), we can perform the complex integration of the convolution term

$$ \begin{aligned}[b] \Delta \tilde{p} \int_{-\infty }^{\infty} e^{-i\omega t}\tilde{\varPhi }H dt &= \Delta \tilde{p} \int_{-\infty }^{\infty} e^{-i\omega t}\frac{1}{ \tau_{\epsilon }}\biggl(1 - \frac{\tau_{\epsilon }}{\tau_{\sigma }}\biggr)e^{-\frac{t}{ \tau_{\sigma }}}H dt \\ &= \Delta \tilde{p}\frac{1}{\tau_{\epsilon }}\biggl(1 - \frac{\tau_{\epsilon }}{\tau_{\sigma }}\biggr) \int_{0}^{\infty} e^{-(i\omega + \frac{1}{ \tau_{\sigma }})t}dt \\ &=\Delta \tilde{p}\frac{1}{\tau_{\epsilon }}\biggl(1 - \frac{\tau_{\epsilon }}{\tau_{\sigma }}\biggr) \frac{1 + i\omega \tau_{\epsilon }}{1 + i\omega \tau_{\sigma }} \end{aligned} $$

(26)

then, from the plane wave form in the Fourier space $(\partial_{t} = -i \omega )$, and one can recast (25) as

$$ (-i\omega + \mathbf{w}.\nabla )^{2} \tilde{p} = \frac{\hat{K_{c}}}{ \rho }\Delta \tilde{p} - \rho c^{2}(-i\omega + \mathbf{w}.\nabla ) \tilde{S}. $$

(27)

The tilde denotes the time Fourier transform and $\hat{K_{c}}$ is the complex bulk modulus, in the case of a single relaxation mechanism, which reads (Carcione et al. 1988)

$$ \hat{K_{c}} = \rho c^{2} \frac{1 + i\omega \tau_{\epsilon }}{1 + i \omega \tau_{\sigma }} = \rho \hat{c}^{2}, $$

(28)

where the complex velocity $\hat{c}$ reads

$$ \hat{c} = c \sqrt{\frac{1 + i\omega \tau_{\epsilon }}{1 + i\omega \tau_{\sigma }} } $$

(29)

The complex equation (27) can then be expanded as

$$ \begin{gathered} \bigl(-\omega^{2} + w_{x}^{2} \partial_{x}^{2} - 2i\omega w_{x} \partial_{x}\bigr) \tilde{p} = \frac{\hat{K_{c}}}{\rho }\bigl( \partial_{x}^{2} + \partial _{z}^{2}\bigr) \tilde{p} - \rho c^{2}(-i\omega + \mathbf{w}.\nabla ) \tilde{S}, \\ \biggl\{ \frac{\rho }{\hat{K_{c}}}\bigl(-\omega^{2} + w_{x}^{2} \partial_{x} ^{2} - 2i\omega w_{x} \partial_{x}\bigr) - \partial_{x}^{2} - \partial_{z} ^{2} \biggr\} \tilde{p} = -\frac{\rho^{2}}{\hat{K_{c}}} c^{2}(-i \omega + \mathbf{w}.\nabla )\tilde{S}, \\ \biggl\{ \frac{\rho }{\hat{K_{c}}}\bigl(\omega^{2} - w_{x}^{2} \partial_{x} ^{2} + 2i\omega w_{x} \partial_{x}\bigr) + \partial_{x}^{2} + \partial_{z} ^{2} \biggr\} \tilde{p} = \frac{\rho^{2}}{\hat{K_{c}}} c^{2}(-i \omega + \mathbf{w}.\nabla )\tilde{S}, \\ \biggl\{ \frac{1}{\hat{c}^{2}}\bigl(\omega^{2} - w_{x}^{2} \partial_{x}^{2} + 2i\omega w_{x} \partial_{x}\bigr) + \partial_{x}^{2} + \partial_{z}^{2} \biggr\} \tilde{p} = \frac{\rho c^{2}}{\hat{c}^{2}} (-i \omega + \mathbf{w}.\nabla )\tilde{S}, \\ \bigl(\partial_{x}^{2} + \partial_{z}^{2} + \hat{k}^{2} + 2i\hat{k}\hat{M} \partial_{x} - \hat{M}^{2}\partial_{x}^{2}\bigr) \tilde{p} = - \frac{ \rho c^{2}}{\hat{c}^{2}} (-i\omega + \mathbf{w}.\nabla )\tilde{S}, \end{gathered} $$

(30)

with $\hat{M}$ the Mach number, that is,

$$ \hat{M} = \frac{w_{x}}{\hat{c}} = \frac{w_{x}}{c \frac{1 + i\omega \tau_{\epsilon }}{1 + i\omega \tau_{\sigma }}}, $$

(31)

where $\hat{c}$ is the complex velocity and $k$ is the wavenumber $\hat{k} = \frac{\omega }{\hat{c}}$. The hat over $\hat{M}$ and $\hat{k}$ is introduced to differentiate from the usual Mach and wave numbers in an inviscid medium.

Finally, we will consider

$$ \tilde{S} = \frac{2i\tilde{A}\hat{K_{c}}}{\omega \rho c^{2} }e^{-i \omega t}\delta_{x} \delta_{z}, $$

(32)

where $\delta_{x}, \delta_{z}$ are the Kronecker delta symbols, and $\tilde{A}$ is the Fourier transform of the amplitude of the time-dependent source function, that is,

$$ A(t) = -2(\pi f_{0})^{2}(t-t_{0})e^{-\{\pi f_{0}(t - t_{0})\}^{2}}. $$

(33)

We use the same method as in Ostashev et al. (2005) to obtain the pressure response in this framework, that is

$$ \tilde{p} = \frac{iA}{2(1-\tilde{M})^{3/2}} \biggl( H_{0}^{1}(\xi ) - \frac{i \tilde{M}\cos\beta }{\sqrt{1 - \tilde{M}^{2}\sin {\beta }^{2}}}H^{1} _{1}(\xi ) \biggr) e^{\frac{-i\tilde{k}\tilde{M}R\cos {\beta }}{1- \tilde{M}^{2}}}, $$

(34)

where $(H^{1}_{i})_{i=1,2}$ are the Hankel functions of first kind, and $\xi$ reads

$$ \xi = \frac{\tilde{k}R\sqrt{1 - \tilde{M}^{2}\sin {\beta }^{2}}}{1- \tilde{M}^{2}}. $$

(35)

Using Hankel function asymptotic (Abramowitz and Stegun 1964) by considering the approximation $|\tilde{k}R| >> 1$, yields

$$ \tilde{p} = \frac{A(\sqrt{1 - \tilde{M}^{2}\sin {\beta }^{2}}- \tilde{M}\cos {\beta })}{\sqrt{2\pi \tilde{k}R}(1-\tilde{M}^{2})(1- \tilde{M}^{2}\sin {\beta }^{2})^{3/4}}e^{\frac{i}{1-\tilde{M}^{2}}(\sqrt{1 - \tilde{M}^{2}\sin {\beta }^{2}}-\tilde{M}\cos {\beta })\tilde{k} R + \frac{i\pi }{4}}. $$

(36)

Appendix 2: Numerical Implementation of ADE-PML for a Medium Subject to Gravity

This boundary condition corresponds to mesh stretching in the boundary layer which requires to recast the derivatives in the new stretched set of complex-valued coordinates (Komatitsch and Martin 2007). We do this using ADE-PML memory variables. In this new reference frame (19) read

$$ \begin{aligned} \partial_{t} p & = -\mathbf{w}.\tilde{\nabla } p - \rho c^{2} \tilde{\nabla } . \mathbf{v} -\tilde{\nabla }p_{0} \mathbf{v} + e_{1}, \\ \partial_{t} \rho_{p} &= -\mathbf{w}.\tilde{\nabla } \rho_{p} - \tilde{\nabla }.(\rho \mathbf{v}), \\ \rho \partial_{t}\mathbf{v} &= - \rho \bigl\{ (\mathbf{v}.\tilde{ \nabla }) \mathbf{w} + (\mathbf{w}.\tilde{\nabla })\mathbf{v} \bigr\} + \tilde{ \nabla } . \boldsymbol{\Sigma } + \tilde{\mathbf{g}}\rho_{p}, \\ \partial_{t} e_{1} &= -\frac{1}{\tau_{\sigma }}\biggl[\biggl(1 - \frac{\tau_{ \sigma }}{\tau_{\epsilon }}\biggr) + e_{1}\biggr], \end{aligned} $$

(37)

with

$$ (\boldsymbol{\Sigma })_{ij} = -p\delta_{ij} + \mu \biggl( \tilde{\partial }_{j}( \mathbf{v}+\mathbf{w})_{i} + \tilde{ \partial }_{i}(\mathbf{v}+ \mathbf{w})_{j} - \frac{2}{3}\delta_{ij}\tilde{\nabla }.\mathbf{v}\biggr) + \eta_{V}\delta_{ij}\tilde{\nabla }.\mathbf{v}. $$

Note that we wrote the pressure equation (37)-1 in terms of the ambient pressure $p_{0}$ rather than in terms of its value coming from hydrostatic equilibrium (i.e. in the interior domain $\nabla p_{0} = \rho \mathbf{g}$). Indeed, as we stretch coordinates the hydrostatic equilibrium equation is also modified as

$$ \tilde{\nabla } p_{0} = \rho \mathbf{g} . $$

(38)

In this new reference frame derivatives will read, for any unknown $\phi$

$$ \tilde{\partial }_{x} \phi = \frac{1}{\kappa_{x}}\partial_{x} \phi + Q _{x}^{\phi } ; \qquad \tilde{\partial }_{z} \phi = \frac{1}{\kappa_{z}}\partial_{z}\phi + Q_{z}^{\phi } , $$

(39)

where the memory variable is denoted $Q_{x,z}^{\phi }$. Memory variables obey a first-order in time differential equation for $x$-derivatives and $z$-derivatives. We will only give details for the $x-$derivative since it is similar for the $z$-derivatives. We consider the four time sub-steps $(t^{n,i})_{i=1,4}$ as $t^{n,i} = \Delta t (n+\nu_{i})$ for the four RK4 stages where $\nu_{i} = 0, \nu_{2} = 0.5, \nu_{3} = 0.5$ and $\nu_{4} = 1$. Then the memory variable $Q_{x}^{\phi }$ obey the following differential equation

$$ Q_{x}^{\phi }\bigl(t^{n,i}\bigr) = b_{x,i} Q_{z}^{\phi }\bigl(t^{n,1}\bigr) + a_{x,i} \partial_{x}\phi \bigl(t^{n,i}\bigr) , $$

(40)

where

$$ a_{x,i} = -d_{x}\frac{\Delta t\nu_{i}}{\kappa_{x}\kappa_{x}(1+\Delta t \theta \nu_{i}(\frac{d_{x}}{\kappa_{x}}+b_{x}))} ; \qquad b_{x} = \frac{ \kappa_{x}-\theta \Delta t\nu_{i}(d_{x} + \kappa_{x}\alpha_{x})}{ \kappa_{x}+\theta \Delta t\nu_{i}(d_{x} + \kappa_{x}\alpha_{x})} , $$

(41)

with

$$ \begin{gathered} \tau_{x} = \frac{1}{\frac{d_{x}}{\kappa_{x}}+\alpha_{x}}; \qquad \kappa _{x} = 1 + ( \kappa_{max}-1)^{m}; \\d_{x} = d_{0} \biggl(\frac{x}{L}\biggr)^{N}; \qquad \alpha_{x} = \alpha_{max}\biggl[1-\biggl(\frac{x}{L}\biggr)^{p} \biggr] , \end{gathered} $$

(42)

where the $\theta$ parameter indicates if the computation is performed explicitly ($\theta = 0$), semi-implicitly ($\theta = 0.5$) or implicitly ($\theta = 1$). Here this parameter is taken as $\theta = 0.5$. Since we are mainly interested in damping acoustic waves (most of the waves reaching the boundary will be acoustic waves) we take the same parameters as in Martin et al. (2010), meaning that

$$ \begin{gathered} N = 2; \qquad m = 2; \qquad p = 1; \qquad \alpha_{max} = \pi f_{0}; \\ \kappa_{max} = 7;\qquad d_{0} = - \frac{(N+1)c_{max,x}log(R_{c})}{2L} , \end{gathered} $$

(43)

where $f_{0}$ is the dominant frequency of the source, $c_{max,x}$ the maximum effective velocity along $x$ (sum of background flow velocity and adiabatic sound speed) and $R_{c}$ the theoretical reflection coefficient taken as $R_{c} = 0.0001$. The same formulation holds for derivatives along $z$.

Appendix 3: “Image” Boundary Conditions

This boundary condition consists in duplicating the physical domain and sources at ghost points located below the bottom boundary in order to obtain a wave interference ensuring zero vertical velocity at the boundary (Morse and Ingard 1968, Sect. 7.4). One enforces that the pressure is equal on each side of the boundary

$$ p(-\mathbf{x}, t) = p(\mathbf{x}, t), \quad \forall \mathbf{x} \in \varOmega \backslash \varGamma $$

(44)

and the sign of gravity and wind is changed in the ghost domain as

$$ \begin{aligned} g_{y}(-\mathbf{x}) &= - g_{y}(\mathbf{x}), \\ w_{x}(-\mathbf{x}) &= -w_{x}(-\mathbf{x}), \quad \forall \mathbf{x} \in \varOmega \backslash \varGamma . \end{aligned} $$

(45)

while all other unknowns are duplicated from the real domain to the ghost one

$$ \begin{aligned} \mu (-\mathbf{x}) &= \mu (\mathbf{x}), \eta_{V}(-\mathbf{x}) = \eta _{V}(\mathbf{x}), \\ \tau_{\epsilon }(-\mathbf{x}) &= \tau_{\epsilon }(\mathbf{x}), \tau_{\sigma }(-\mathbf{x}) = \tau_{\sigma }(\mathbf{x}), \quad \forall \mathbf{x} \in \varOmega \backslash \varGamma \end{aligned} $$

(46)

As a consequence, the pressure being equal on both sides of the boundary, the pressure gradient in the momentum equation is zero at the boundary. Also, setting opposite gravity directions above and below the boundary leads to zero gravity at the surface, hence the gravity term in the momentum equation vanishes. Therefore, at the boundary, the momentum equation is verified since the zero velocity condition is imposed. Finally, we impose the boundary to be motionless and we choose

$$ \begin{aligned} v_{x}(x, -z) &= -v_{x}(x, z), \\ v_{z}(x, -z) &= -v_{z}(x, z), \quad \forall z \in \varOmega. \end{aligned} $$

(47)