# Full waveform inversion of infrasound: sensitivity kernels from the adjoint method along with nonlinear inversion tests

Gérier Solène, Martin Roland, Garcia Raphaël F.

► To cite this version:

Gérier Solène, Martin Roland, Garcia Raphaël F.. Full waveform inversion of infrasound: sensitivity kernels from the adjoint method along with nonlinear inversion tests. Geophysical Journal International, 2026, 244 (3), 10.1093/gji/ggaf472. hal-05510792

HAL Id: hal-05510792

https://hal.science/hal-05510792v1

Submitted on 14 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Full waveform inversion of infrasound: sensitivity kernels from the adjoint method along with nonlinear inversion tests

S. Gerier<sup>1</sup>, R. Martin<sup>2</sup> and R. F. Garcia<sup>1</sup>

## SUMMARY

Atmospheric models are based on various types of geophysical data, including lidar and radar. Infrasounds, acoustic waves that can propagate over large distances, have not yet been used in atmospheric models, although they provide valuable information. Besides their sensitivity to atmospheric phenomena such as gravity waves, infrasound also presents the advantage of being omnipresent. Previous studies explored the use of infrasound packet arrival properties for model estimation. However, properties such as arrival times present less information than full waveforms. We aim here to investigate, for the first time, the sensitivity of a full infrasound waveform to model parameters and to use these sensitivities in an inverse problem to recover atmospheric structure. For this purpose, infrasound propagation is modelled by Euler equations (i.e. Navier-Stokes equations in the absence of attenuation effects), and discretization is carried out here using the finite-differences method. Waveform sensitivity to atmospheric parameters is computed through the adjoint method via a novel and optimized double checkpointing-based procedure and validated by comparison with a small perturbation method. As an illustration, these sensitivity kernels are computed for the idealized case of an explosion in Finland, recorded by a station of the Comprehensive Nuclear-Test-Ban Treaty (CTBT). These first results demonstrate the high sensitivity of infrasound waveforms to the atmospheric perturbations generated by gravity waves. Moreover, the sensitivity kernels of infrasound waveforms allow us to recover the variations of model parameters by solving an inverse problem. To demonstrate this capability, full waveform nonlinear inversions are performed using the Limited Broyden-Fletcher-Goldfarb-Shanno method: wind and sound speed profiles are inverted for a test case with idealized conditions and a synthetic data set. These estimates of infrasound sensitivity kernels are closing a knowledge gap that allows the use of infrasound full waveforms to constrain atmospheric models.

**Key words:** Inverse theory; Numerical modelling; Infrasound; Waveform inversion.

## 1 INTRODUCTION

Atmospheric models are of particular interest in many fields. An example is the monitoring of explosions using infrasound by the International Monitoring System (IMS). Atmospheric models that accurately describe the temperature, wind and acoustic wave speed of the atmosphere are required to accurately predict the propagation of acoustic waves and the source localization, since the propagation is directly impacted by atmospheric conditions (Drob *et al.* 2009). Furthermore, numerical weather predictions for forecasting tropospheric conditions are improved if the state of the upper atmosphere is better described due to coupled troposphere-stratosphere effects (Charlton *et al.* 2004; Gerber *et al.* 2009). This additional argument justifies research into improving high-altitude models. We are thus interested in the problem of estimating the atmospheric properties specified by a vertical profile on a regional scale to fully explain a set of infrasound observations.

The problem of atmospheric model estimation has been tackled for many data sets: rocket sonde, radar, Lidar, etc. (Blom & Marcillo 2017; Blanc *et al.* 2018; Eyre *et al.* 2022). The last decade was marked by significant improvements in the data assimilation capabilities for the lower, middle and upper atmospheres (Drob 2019). However, the background state described by the atmospheric model is subject to biases at altitudes where the observations become sparser (Simmons *et al.* 2020; Vera Rodriguez *et al.* 2020; Laloyaux *et al.* 2022).

<footnote n="1"><sup>1</sup>Institut Supérieur de l'Aéronautique et de l'Espace (ISAE-SUPAERO), Université de Toulouse, 10 Ave M. Pelegrin, 31400 Toulouse, France. E-mail: solene.gerier@isae-supaero.fr</footnote>

<footnote n="2"><sup>2</sup>Laboratoire de Géosciences Environnement Toulouse GET, UMR CNRS 5563, Observatoire Midi-Pyrénées, Université de Toulouse, 14 avenue Edouard Belin, 31400 Toulouse, France</footnote>
Infrasound waves can provide complementary information to characterize the middle atmosphere through remote sensing. Infrasound are acoustic waves below 20 Hz and are produced by natural sources (microbaroms, quakes, volcanoes, meteorites) as well as by anthropogenic sources (mine explosions, nuclear explosions) (Campus & Christie 2009). These waves can propagate over long distances due to their low attenuation and the presence of natural atmospheric waveguides (Drob et al. 2003), which are formed by vertical variations of temperature and winds in the stratosphere and the thermosphere.

The sensitivity of the infrasound to the atmospheric structure has been demonstrated. Donn & Rind (1971) relate the temporal variability of infrasound to temperature and wind seasonal and diurnal variations. Other studies show the seasonal variation of infrasound observations (Garcés et al. 2004; Le Pichon et al. 2006; Fricke et al. 2014).

Similarly to seismic waves, infrasound, observed at some stations, has integrated along its path the properties of the medium of propagation, such as temperature and wind of the atmosphere. Moreover, since infrasound signals propagating in the stratospheric waveguide are refracted from altitude near the stratopause to the ground (Waxler et al. 2011; Marty 2019), observations of infrasound are sensitive to the atmospheric structure near the turning point, where other direct measurements are difficult and sparse (Blom & Marcillo 2017; Assink et al. 2019). These different elements show that infrasound observations are a valuable source of information to complete the characterization of the atmospheric structure.

The use of infrasound measurements to reconstruct the atmospheric model is an inverse problem. A detailed review of infrasound-based inverse problems for probing atmospheric structure has been performed by Assink et al. (2019).

Some works are distinguished by the parametrization of the atmospheric model. Seismic studies generally decompose their domain into unit cells, such as a grid. Parameters are defined in each cell, such as in Köhn et al. (2012) and Monteiller et al. (2015). Le Pichon et al. (2005) used scaled Gaussian correction factors and Drob et al. (2009) introduced empirical orthogonal functions to describe the variations of the model along the altitude. The two above-mentioned approaches enable us to model the atmospheric structure with a few parameters and reduce the space of the solution for the inverse problem. Furthermore, constraints on the parameters can be added, as in Vera Rodriguez et al. (2020), which bounded the values of the physical properties with respect to the variations between members of an ensemble atmospheric model.

The infrasound used to solve the inverse problem in recent studies can come from different sources. Le Pichon et al. (2005) and Assink et al. (2014) studied infrasound from volcanoes. Infrasound from explosions, such as the explosions generated by the ammunition destruction at the Hukkakero site in Finland, is often used because these signals can be well characterized (Amezcua et al. 2020; Vera Rodriguez et al. 2020; Vorobeva et al. 2023). Microbaroms are also considered as sources in the inversion studies (Lalande et al. 2012; Vorobeva et al. 2024). In particular, in Arrowsmith et al. (2013), Blom & Marcillo (2017) and Vorobeva et al. (2024), no assumptions are made regarding the source location.

Furthermore, identifying the model properties has been solved in different ways. Drob et al. (2009) and Lalande et al. (2012) used least-squares formalism to evaluate the distance between the output of the forward modelling and the data, and minimize this distance by an iterative procedure. Vera Rodriguez et al. (2020) aims to minimize a cost function, defined as a least-square or any other ad-hoc metric, and uses a heuristic solver to search for the solution in their constrained solution space. Grid search methods in combination with Bayesian statistics have been used to assess the posterior probability of stratospheric and thermospheric conditions (Assink et al. 2013, 2014). Exploration of solution space using Monte-Carlo methods was done by Froment (2023) and Froment et al. (2024) for the Martian Planetary Boundary Layer (PBL) and allows the authors to estimate the atmospheric state. However, the Martian PBL is 10 km, and applying the technique to Earth's atmosphere would require much more computational time and could become prohibitive in terms of computational resources in some cases. Stochastic-based machine learning methods have also been used to estimate the stratospheric polar vortex strength (Vorobeva et al. 2024). Amezcua et al. (2020, 2024) used a data assimilation framework and applied an ensemble Kalman filter for the analysis part.

In addition, these works also distinguished themselves by the data and the associated forward modelling. Amezcua et al. (2020) used a geometric formulation of the wind to estimate crosswind from the back azimuth. Drob et al. (2009) and Lalande et al. (2012) did atmospheric tomography using traveltimes and exploiting the acoustic geometry approximation. In particular, Lalande et al. (2012) used a linearized version of the acoustic geometry approximation. A full waveform modelling of the pressure has been implemented by Chunchuzov et al. (2015) to retrieve the fine-scale structure of wind and identify gravity waves. The waveform is estimated using the nonlinear geometrical acoustics approximation (Chunchuzov et al. 2013).

To the best of our knowledge, the linearized Euler equations, which enable the full waveform of infrasound to be modelled, are not used in infrasound-based inverse problems. Nevertheless, the linearized Euler equations have the advantage of providing full waveforms and not just arrival times. If we refer to seismology, which is a similar domain, seismic waves are affected, as well as acoustic waves, by scattering effects, and therefore waveforms integrate more information than traveltimes (Fichtner 2011). Hence, full waveform inversion has thus been developed to image the Earth's structure in seismology (Liu & Tromp 2008; Virieux & Operto 2009; Peter et al. 2011; Schafer et al. 2014; Monteiller et al. 2015; Komatitsch et al. 2016).

To consider effects such as infrasound scattering in the inverse problem, we aim at applying the full waveform inversion concept designed in seismology. This concept aims at finding the model that best explains the waveform under the constraint of state equations. It is based on the adjoint-state method for probing the sensitivity of waveforms to the model (Tarantola 1986; Tromp et al. 2005; Liu & Tromp 2008), and was developed in the textbook of Tarantola (2005) for geophysical applications.
Sensitivity kernels have also been developed in ocean acoustics with methods and assumptions similar to those used in solid Earth seismology (Skarsoulis & Cornuelle 2004; Hermand 2006; Aulanier et al. 2013; Zhao & Wang 2015; Cornuelle & Skarsoulis 2022). In particular, the effect of ocean currents is neglected, and the ocean properties are assumed not to vary during acoustic wave propagation.

To the best of our knowledge, few studies in atmospheric acoustics have used the adjoint method. More broadly in other areas, we can cite the work of Lemke & Sesterhenn (2016) and Stein et al. (2019) in the sound reinforcement domain, or Spieser & Bailly (2018, 2020) and Spieser et al. (2020), who worked on the establishment of an adjoint of the linearized Euler equations. They are then able to use the reciprocity principle to compute the mean flow refraction effect on sound emitted by sources within a jet. However, the authors did not address the notion of sensitivity kernels with respect to model parameters. Adapting the adjoint method derived from seismology or ocean acoustics to infrasound requires a specific treatment due to the presence of wind in the linearized Euler equations, which introduces asymmetry in the equations.

The optimization framework can be solved using gradient methods. Here, many options are available. Different strategies are presented in Nocedal & Wright (1999). Here, we use a line search framework. In seismology, full waveform inversion relies on the Steepest Descent, (non-)Conjugate Gradient, Newton, Gauss–Newton, or Quasi-Newton methods to iteratively reduce the cost function defined as the misfit between the true measured data and the data computed for the current model (Tromp 2019). In our study, the misfit improvement direction is acquired by the Limited Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) method, one of the most popular descent directions based on Quasi-Newton methods (Nocedal & Wright 1999). As proof of concept, the algorithm is tested under ideal conditions, where receivers are located on the ground and at different altitudes, the infrasound waveform data set is idealized, and the source time function is assumed to be perfectly known.

In this paper, we aim at developing a method to obtain the sensitivity of acoustic waveforms to an atmospheric model and to use this information to solve the inverse problem. We parametrize the atmospheric structure with a 1-D model specifying the temperature, wind, pressure and density conditions. A forward modelling operator is solved by a non-centred finite-difference approximation of the linearized Euler equations. Pressure variations resulting from infrasound waves constitute our observation data set. The inverse problem is solved using the adjoint-state method to compute the sensitivity of the data to the model, and optimization is performed by the descent direction of the L-BFGS formula, defined in Nocedal & Wright (1999), corresponding to a Quasi-Newton method adapted for large problems.

The main contribution of this paper is the development of sensitivity kernels applied to infrasound waveforms. This includes establishing adjoint equations for the inverse problem under consideration and adapting the non-centred finite-difference schemes for wind advection terms to achieve a more stable and accurate resolution of the adjoint equations. The process for obtaining the infrasound wavefield throughout the domain backward-in-time (i.e. from the final time to the initial time) has been revised: we propose a novel and optimized double checkpointing-based method to calculate the infrasound wavefield backward-in-time efficiently with limited memory.

This article is divided into four parts: First, we describe the state problem, which models wave propagation in the atmosphere for a given current model. Secondly, we focus on adjoint theory for the computation of sensitivity kernels. For illustrative purposes, sensitivity kernels are computed for simplified and realistic cases in the context of explosions in the atmosphere. Thirdly, we explain how to solve the inverse problem and present inversion results for two idealized test cases. Finally, we discuss the results obtained for simplified and realistic cases, and conclude on the potential of the different methods used to solve the inverse problem.

## 2 STATE PROBLEM

In this section, we detail the state problem, also called the forward problem, which describes wave propagation for a given atmospheric model. The state problem produces an estimation of the observations used for the adjoint problem and, more generally, for the inverse problem. The assumptions, equations and their numerical resolution are highlighted here, whereas definitions of physical quantities are in Table A1 in Appendix A. The definition of the state problem (and its numerical version) is important because it is necessary to write the inverse problem.

### 2.1 Assumptions

Here, we list all the assumptions used to model the propagation of an infrasonic wave.

1) The fluid is compressible and Newtonian.

2) The fluid is a continuum field.

3) The fluid is an ideal gas. Consequently, the following relation links the wave speed $c$, the pressure $p$ and the density $\rho$: $c(x, t) = \sqrt{\gamma R T(x, t)} = \sqrt{\gamma \frac{p(x, t)}{\rho(x, t)}}$ for any point $x$ and time $t$.

4) No viscosity, no thermal diffusion, nor any other attenuation mechanism is considered.

5) A physical quantity $q$ of the fluid (such as density, pressure and velocity) can be approximated by the sum of a background stationary state $q_0$ and a perturbation $q'$:

$$q(x, t) = q_0(x) + q'(x, t)$$

$q'$ is considered small compared to the background state $q_0$: $q'(t) \ll q_0, \forall t > 0$.
6) The background state is vertically stratified. Thus, each quantity $q_0(\mathbf{x}) = q_0(x_1, x_2) = q_0(0, x_2)$, $\forall x_2$. $\mathbf{x} = (x_1, x_2)^T$ is a point of $\Omega$ where $x_1$ is the horizontal range and $x_2$ is the altitude.

7) The wind speed is written as $\mathbf{v}_0$. Thus, it corresponds to the background contribution of the particle velocity $\mathbf{v}$ of the fluid. The horizontal wind is $\mathbf{v}_0^1$. The vertical component of the wind speed $\mathbf{v}_0^2$ is supposed to be zero. The wind speed is $\mathbf{v}_0(\mathbf{x}) = \mathbf{v}_0(x_1, x_2) = (\mathbf{v}_0^1(0, x_2), 0)^T$.

8) An explosion occurs in the fluid at a given time. This phenomenon is modelled by the addition of material and energy to the domain, without the presence of a force that modifies the particle velocity. We suppose no external contribution in terms of the momentum.

9) There is no disturbance at time $t = 0$. The perturbations do not reach the atmospheric boundary before $T_{\max}$. At the ground boundary $\Gamma_1$, the velocity perturbation and wind speed are set to zero, which defines the Dirichlet conditions for modelling infrasound propagation.

The previous assumptions allow to use a simplified form of the Navier–Stokes equations to describe the infrasound propagation, called Euler equations (Landau & Lifshitz 1987; Batchelor 2000; Rieutord 2016). Indeed, in contrast to the Navier–Stokes equations, the Euler equations do not take into account any attenuation effects (neither viscosity nor diffusivity nor thermal dissipation terms).

In contrast to Tarantola (1984) and Gauthier (1986), who solved the inverse problem by considering acoustic equations, we consider the linearized Euler equations, which include wind effects. Wind introduces asymmetry (i.e. the properties of the atmospheric fluid vary with direction and sense). This asymmetry has a non-negligible effect on the adjoint method, and we wish to verify that this method still works in the presence of wind.

## 2.2 Linearized Euler equations

Assuming small disturbances, we can linearize the Euler equations, which is a standard approach in acoustics and where perturbations are distinguished from background quantities.

Using the expression of density $\rho$, particle velocity $\mathbf{v}$ and pressure $p$ given by assumption (H5) in the Euler equations (cf. eq. presented in Appendix B1), we have:

$$ 
\left\{ 
\begin{aligned}
\partial_t \rho' + \rho_0 \nabla \cdot \mathbf{v}' + \mathbf{v}_0 \cdot \nabla \rho' + \mathbf{v}' \cdot \nabla \rho_0 - f_\rho &= 0 \\
\rho_0 \partial_t \mathbf{v}' + (\rho_0 \mathbf{v}_0 \cdot \nabla) \mathbf{v}' + \nabla p' + (\rho_0 \mathbf{v}' \cdot \nabla) \mathbf{v}_0 - \rho' \mathbf{g} &= 0 \\
\partial_t p' + \gamma \rho_0 \nabla \cdot \mathbf{v}' + \mathbf{v}_0 \cdot \nabla p' + \mathbf{v}' \cdot \nabla p_0 - f_p &= 0
\end{aligned}
\right. 
\qquad (1) 
$$

with $f_p = c_0^2 f_\rho$, and the wave speed $c_0 = \frac{\gamma p_0}{\rho_0}$.

Using assumptions (H9), the eq. (1) are respectively associated with the following initial conditions and boundary conditions:

$$ 
\left\{ 
\begin{aligned}
\rho'(\mathbf{x}, 0) &= 0 \quad \forall \mathbf{x} \in \Omega \\
\mathbf{v}'(\mathbf{x}, 0) &= 0 \quad \forall \mathbf{x} \in \Omega \\
p'(\mathbf{x}, 0) &= 0 \quad \forall \mathbf{x} \in \Omega
\end{aligned}
\right. 
\qquad (2) 
$$

$$ 
\left\{ 
\begin{aligned}
\rho'(\mathbf{x}, t) &= 0 \quad \forall \mathbf{x} \in \Gamma \setminus \Gamma_1 \quad \forall t \in [0, T_{\max}] \\
\mathbf{v}'(\mathbf{x}, t) &= 0 \quad \forall \mathbf{x} \in \Gamma \quad \forall t \in [0, T_{\max}] \\
p'(\mathbf{x}, t) &= 0 \quad \forall \mathbf{x} \in \Gamma \setminus \Gamma_1 \quad \forall t \in [0, T_{\max}]
\end{aligned}
\right. 
\qquad (3) 
$$

In the following, we will refer to the system of eq. (1) by $\mathcal{L}(\mathbf{u}; \mathbf{m}) = 0$, where $\mathbf{m} = (\rho_0, \mathbf{v}_0^1, \rho_0)^T$ is called ‘model’ and describes the background state, and $\mathbf{u} = (\rho', \mathbf{v}', p')^T$ is called the ‘state variable’ and describes the wave perturbations. $\mathcal{L}(\mathbf{u}; \mathbf{m}) = 0$ are called the state equations.

Appendix B2 explains the discretization and resolution of the system of eq. (1).

## 2.3 Inverse problem definition

To solve the inverse problem, we aim at finding the model that best fits the observations, under the assumption that the model, describing the background over the space, and the observations satisfy the system of eq. (1).

Let $\chi$, such that:

$$ 
\mathbf{u} = \begin{pmatrix} \rho' \\ \mathbf{v}' \\ p' \end{pmatrix} \mapsto \chi(\mathbf{u}) = \frac{\|\mathcal{A}(p' - p^{\text{obs}})\|_{D}^{2}}{2\|p^{\text{obs}}\|_{D}^{2}}. 
\qquad (4) 
$$

with $p^{\text{obs}}$ the pressure observed at receivers. $\| \cdot \|_D : f \to \sqrt{\int_{\Omega} \int_{0}^{T_{\max}} f^2 dt d\Omega}$ is the norm over the space of $D$ (i.e. space of the data), and $A : t \to A(t) \in \mathbb{R}$ is the square root of the apodization function, defined from the Hann window as follows:

$$ 
A^2(t) = \begin{cases} 
0.5 \left(1 - \cos\left(\frac{\pi t}{t_0}\right)\right) & \text{if } t < t_0, \\
0.5 \left(1 - \cos\left(\frac{\pi (T_{\max} - t)}{t_0}\right)\right) & \text{if } t > T_{\max} - t_0, \\
1 & \text{otherwise.}
\end{cases} 
\qquad (5) 
$$
$\chi$ measures the normalized discrepancy between the observed and simulated waves (obtained by solving eq. 1) recorded on a given set of receivers. We call $\chi$ the misfit. The discrepancy between observations and synthetic waveforms is modulated by the function $A: t \mapsto \mathbb{R}$. The signal is apodized at the start and end of the temporal domain to avoid some numerical instabilities and spurious high-frequency signals (as in Xie et al. 2014).

Furthermore, it is assumed here that $\mathbf{u}$ is acquired by solving $\mathcal{L}(\mathbf{u}; \mathbf{m}) = 0$.

The inverse problem is then written mathematically as the following minimization problem:

$$ \min_{\mathbf{m}} \chi(\mathbf{u}) \quad \text{s.t.} \quad \mathcal{L}(\mathbf{u}; \mathbf{m}) = 0. \qquad (6) $$

Note that although the equations have been linearized with respect to the disturbance variables $\mathbf{u}$, the propagation equations remain nonlinear with respect to the model $\mathbf{m}$. In other words, the linearized Euler-based wave equation operator relating the model to the data is nonlinear. As a result, the objective function itself is nonlinear with respect to model $\mathbf{m}$. Consequently, the optimization problem (6) is nonlinear.

# 3 SENSITIVITY KERNELS COMPUTATIONS BY ADJOINT METHOD

## 3.1 Adjoint state theory and sensitivity kernels

The inverse problem, expressed by eq. (6), can be formulated as a least-squares problem (Tarantola 1987) and can be solved locally and iteratively by the gradient method. The derivative of the misfit $\chi$ with respect to $\mathbf{m}$, also called Fréchet derivatives, is the sensitivity kernel. Fréchet derivatives were calculated in Lalande et al. (2012) to identify the best atmospheric model that explains pressure records. The misfit of their inverse problem is similar to our study, and forward modelling differs because it is based on ray tracing.

As with conventional derivatives, the sensitivity kernels express the sensitivity to change of the function $\chi$ with respect to its input $\mathbf{m}$. Here, the derivative of $\chi$ cannot be directly computed, in particular because of the constraint and the implicit function relating $\mathbf{m}$ to $\mathbf{u}$. However, the adjoint-state method allows for an efficient computation of the sensitivity kernel.

The adjoint-state method adds the constraint $\mathcal{L}(\mathbf{u}; \mathbf{m}) = 0$ in the misfit through the augmented Lagrangian method. The problem is then expressed as

$$ \min_{\mathbf{m}, \mathbf{u}, \lambda} \tilde{\chi}(\mathbf{m}, \mathbf{u}, \lambda) \quad \text{with} \quad \tilde{\chi}(\mathbf{m}, \mathbf{u}, \lambda) = \chi(\mathbf{u}) + \langle \lambda, \mathcal{L}(\mathbf{u}; \mathbf{m}) \rangle_D \qquad (7) $$

where $\lambda$ is called the adjoint variable, and $\lambda = (\lambda_1, \lambda_2, \lambda_3)^T$, with $\lambda_1, \lambda_3 \in \mathbb{R}$, $\lambda_2 \in \mathbb{R}^2$. The scalar product $\langle f, g \rangle_D$ is defined as the integral over space and time of the product of the functions $f$ and $g$.

In contrast to the $\chi$ function called ‘misfit’, we call the $\tilde{\chi}$ function a ‘cost function’, because this is the one we are effectively minimizing.

The optimality conditions associated with the problem (7) impose that each component of the gradient of the cost function $\tilde{\chi}$, evaluated at the optimal point, is zero. Under these conditions, the optimal values of $\lambda$ and $\mathbf{u}$ can be identified.

Imposing that $\partial_{\lambda} \tilde{\chi}$ the derivative of $\tilde{\chi}$ with respect to $\lambda$ is zero implies that $\mathcal{L}(\mathbf{u}; \mathbf{m}) = 0$. Hence, $\mathbf{u}$ can be found by solving the state problem for the current model $\mathbf{m}$. We note this solution $\mathbf{u}_m$.

$\partial_{\mathbf{u}} \tilde{\chi}$ the derivative of $\tilde{\chi}$ with respect to each component of $\mathbf{u}$ (that is, with respect to $\rho'$, $\mathbf{v}'$ and $p'$) is defined as follows:

$$ \langle \partial_{\mathbf{u}} \tilde{\chi}(\mathbf{m}, \mathbf{u}, \lambda), \delta \mathbf{u} \rangle_D = \langle \partial_{\mathbf{u}} [\chi(\mathbf{u}) + \langle \lambda, \mathcal{L}(\mathbf{u}; \mathbf{m}) \rangle_D], \delta \mathbf{u} \rangle_D, $$

where $\delta \mathbf{u} = (\delta \rho', \delta \mathbf{v}', \delta p')^T$.

Using the definition of the scalar product, the optimality condition can be written as follows:

$$ \langle \partial_{\mathbf{u}} \tilde{\chi}(\mathbf{m}, \mathbf{u}, \lambda), \delta \mathbf{u} \rangle_D = \langle \delta \mathbf{u}, \mathcal{L}^{\dagger}(\lambda; \mathbf{m}) \rangle_D = 0, \qquad (8) $$

with $\mathcal{L}^{\dagger}$ defining a set of equations whose $\lambda(\mathbf{x}, T_{\max} - t) = (\rho^{\dagger} \mathbf{v}^{\dagger} p^{\dagger})^T(\mathbf{x}, t)$ is the unknown. These equations are called ‘adjoint equations’ and are developed as follows $\forall \mathbf{x} \in \Omega, \forall t \in [0, T_{\max}]$

$$ \left\{ \begin{array}{l} \partial_t \rho^{\dagger} - \mathbf{v}_0 \partial_{x_1} \rho^{\dagger} - \mathbf{v}^{\dagger} \cdot \mathbf{g} = 0 \\ \partial_t \rho^{\dagger} \mathbf{v}^{\dagger} - \rho_0 \nabla \rho^{\dagger} - \rho_0 \mathbf{v}_0^{-1} \partial_{x_1} \mathbf{v}^{\dagger} + \rho_0 (\nabla \mathbf{v}_0^{-1}) \mathbf{v}^{\dagger -1} - \nabla(\gamma p_0 p^{\dagger}) + p^{\dagger} \nabla p_0 = 0, \\ \partial_t p^{\dagger} - \nabla \cdot \mathbf{v}^{\dagger} - \mathbf{v}_0^{-1} \partial_{x_1} p^{\dagger} = f_{\rho}^{\dagger} \end{array} \right. \qquad (9) $$

where $f_{\rho}^{\dagger}(\mathbf{x}, t) = \mathcal{A}^2(T_{\max} - t) \frac{\langle p'(\mathbf{x}, T_{\max} - t) - p^{\text{obs}}(\mathbf{x}, T_{\max} - t) \rangle_{\mathcal{M}}}{\|p^{\text{obs}}(\mathbf{x}, T_{\max} - t)\|^2} \Big|_{\mathbf{x}=\mathbf{x}_r}$ is the adjoint source term. Moreover, the presence in the adjoint source term function of the Hann apodization function, defined by eq. (5), allows to apodize the signal at the start and end of the time interval to avoid numerical instabilities (similarly as in Xie et al. 2014).

The equations of the system (9) are associated with the initial and boundary conditions expressed, respectively, as follows:

$$ \left\{ \begin{array}{l} \rho^{\dagger}(\mathbf{x}, 0) = 0 \quad \forall \mathbf{x} \in \Omega, \\ \mathbf{v}^{\dagger}(\mathbf{x}, 0) = 0 \quad \forall \mathbf{x} \in \Omega, \\ p^{\dagger}(\mathbf{x}, 0) = 0 \quad \forall \mathbf{x} \in \Omega. \end{array} \right. \qquad (10) $$
$$ \left\{ \begin{array}{ll} \mathbf{v}^{\dagger 2}(\mathbf{x}, t) = 0 & \forall \mathbf{x} \in \Gamma_1, \forall t \in [0, T_{\max}], \\ \rho^{\dagger}(\mathbf{x}, t) = p^{\dagger}(\mathbf{x}, t) = \mathbf{v}^{\dagger 1}(\mathbf{x}, t) = \mathbf{v}^{\dagger 2}(\mathbf{x}, t) = 0 & \forall \mathbf{x} \in \Gamma \setminus \Gamma_1, \forall t \in [0, T_{\max}]. \end{array} \right. \qquad (11) $$

with $\Gamma_1$ the boundary corresponding to the ground interface.

Detailed calculations of the sensitivity kernels are provided in Appendix C1.

The solution of the system of eq. (9) associated with the conditions given by eqs (10) and (11) is denoted $\lambda_m$.

The operator $\mathcal{L}$, characterizing this inverse problem, is not self-adjoint (i.e. the system of adjoint eq. 9 is different from the system of state eq. 1). Different reasons explain this: considering wind, gravity and a stratified model makes the problem directional and asymmetric. Moreover, considering a non-compressible fluid, which implies a non-constant density in time, introduces also asymmetry in the problem, and explains why the forward and adjoint are not self-adjoint. Besides, in absence of gravity and motion of the background (no wind), and assuming constant densities, we can retrieve the purely acoustic wave equation from the linearized Euler equations. Yet, the wave equation is self-adjoint as it has been shown in seismology (Fichtner et al. 2006). Therefore, having a non-self-adjoint problem leads to greater complexity for solving the adjoint wavefield.

From a specific model $\mathbf{m}$ and knowing $\lambda$ and $\mathbf{u}$, the derivative of $\tilde{\chi}$ with respect to each component of $\mathbf{m}$ can be evaluated. It can be written as

$$ \langle \partial_m \tilde{\chi}(\mathbf{m}, \mathbf{u}, \lambda), \delta \mathbf{m} \rangle_M = \left\langle \partial_m (\langle \mathcal{L}(\mathbf{u}; \mathbf{m}), \lambda \rangle_D), \delta \mathbf{m} \right\rangle_M. \qquad (12) $$

Note that the scalar product $\langle \cdot, \cdot \rangle_M$ in eq. (12) is defined for the space in which the model $\mathbf{m}$ is defined. Hence, it corresponds to the integral over the second component $x_2$.

We introduce the notation $\mathcal{K}$, such that $\partial_m(\langle \mathcal{L}(\mathbf{u}; \mathbf{m}), \lambda \rangle_D) = \mathcal{K}(\mathbf{m}, \mathbf{u}, \lambda)$. $\mathcal{K}$ is called the sensitivity kernel and is defined as follows:

$$ \mathcal{K} = \int_{x_1} \int_{0} T_{\max} \begin{pmatrix} K_{\rho_0} \\ K_{v_0^1} \\ K_{p_0} \end{pmatrix} dt dx_1. \qquad (13) $$

with

$$ K_{\rho_0} = \rho^{\dagger} \nabla \cdot \mathbf{v}' - \partial_{x_2}(\rho^{\dagger} \mathbf{v}'^2) - \mathbf{v}^{\dagger} \cdot \partial_t \mathbf{v}' + \mathbf{v}^{\dagger} \cdot (\mathbf{v}_0^1 \partial_{x_1} \mathbf{v}') + \mathbf{v}^{\dagger} \mathbf{v}'^2 \partial_{x_2} \mathbf{v}_0^1 + p^{\dagger} \frac{\gamma \rho_0}{\rho_0^2} f_{\rho}, \qquad (14) $$

$$ K_{v_0^1} = \rho^{\dagger} \partial_{x_1} \rho' + \mathbf{v}^{\dagger} \cdot (\rho_0 \partial_{x_1} \mathbf{v}') - \partial_{x_2} (\rho_0 \mathbf{v}'^2 \mathbf{v}'^1) + p^{\dagger} \partial_{x_1} p', \qquad (15) $$

$$ K_{\rho_0} = p^{\dagger} \gamma \nabla \cdot \mathbf{v}' - \partial_{x_2} (p^{\dagger} \mathbf{v}'^2) - p^{\dagger} \frac{\gamma}{\rho_0} f_{\rho}. \qquad (16) $$

To simplify the notation of the eqs (14)–(16), we do not specify the evaluation in time and space of the different variables. $\rho'$, $\mathbf{v}'$, $p'$, $f_{\rho}$ are evaluated at $T_{\max} - t$ and $\rho^{\dagger}$, $\mathbf{v}^{\dagger}$, $p^{\dagger}$ at $t$. Consequently, the sensitivity kernel corresponds to a correlation between the adjoint variable $\lambda$ and the state variable $\mathbf{u}$.

Appendix C2.1 details the computation steps to obtain the expressions of $K_{\rho_0}$, $K_{v_0^1}$ and $K_{\rho_0}$.

The sensitivity kernel is valuable information to solve the problem (7), but also the problem (6), since $\partial_m \tilde{\chi}(\mathbf{m}, \mathbf{u}_m, \lambda_m)$ is equal to $\partial_m \chi(\mathbf{u}_m)$.

## 3.2 Model parametrization

To simplify sensitivity kernel computations, we chose to express the set of eq. (1) as a function of density, wind speed and pressure, and therefore our model is also described by these parameters. However, it is one parametrization among several others, and other options are available. Some parametrizations are more favourable than others because they allow to obtain an overall good solution more quickly than with other parametrizations, since the wave is more sensitive to some parameters than others. This was studied in the seismic case in Köhn et al. (2012). The choice of the parametrization is thus critical.

Here, we aim to express the sensitivity kernels for the parametrization $(\ln(\rho_0), \mathbf{v}_0^1, c_0)$. Indeed, the logarithm of density is preferable to density since, given that density changes exponentially with altitude, the variability of density is much greater than that of wind and wave speed. The logarithm resizes the values of the density to achieve a better balance between the model variables.

Exploiting the relation between $\rho_0$, $p_0$ and $c_0$ (i.e. $c_0 = \sqrt{\gamma \frac{\rho_0}{\rho_0^2}}$ provided in assumption H3), we obtain the following relations:

$$ \hat{K}_{\ln(\rho_0)} = \rho_0 K_{\rho_0} + p_0 K_{p_0}, \quad \hat{K}_{v_0^1} = K_{v_0^1}, \quad \hat{K}_{c_0} = 2 \frac{2c_0 \rho_0}{\gamma} K_{p_0}, \qquad (17) $$

where the letter $K$ stands for the sensitivity kernels in the first parametrization $(\rho_0, \mathbf{v}_0^1, p_0)$ and $\hat{K}$ stands for the sensitivity kernels in the parametrization $(\ln(\rho_0), \mathbf{v}_0^1, c_0)$.
## 3.3 Numerical implementation

### 3.3.1 Numerical grid

The adjoint state and the sensitivity kernels can be computed under their discrete form. Concerning the adjoint state $\lambda(x, t) = (\rho^\dagger \mathbf{v}^\dagger p^\dagger)^T(x, T_{\max} - t)$, it is discretized in the same manner as $\mathbf{u}(x, t) = (\rho' \mathbf{v}' p')^T(x, t)$ (see Appendix B2). Thus, $\forall(i, j, n) \in [0, L_1] \times [0, L_2] \times [0, N]$ (where the notation $[\cdot, \cdot]$ denotes integer-only intervals):

(i) $\rho^\dagger$ and $p^\dagger$ on point $(i, j - \frac{1}{2})$ at time $(N - n)\Delta t$;

(ii) the first component of $\mathbf{v}^\dagger$ on point $(i - \frac{1}{2}, j - \frac{1}{2})$ at time $(N - n + 1/2)\Delta t$;

(iii) the second component of $\mathbf{v}^\dagger$ on point $(i, j)$ at time $(N - n + 1/2)\Delta t$.

Concerning the sensitivity kernels, $K_{\rho_0}$ and $K_{\rho_0}$ are defined at $(i, j - \frac{1}{2}, n + 1/2)$, and $K_{v_0^1}$ is defined at $(i - 1/2, j - \frac{1}{2}, n + 1/2)$, $\forall(i, j, n) \in [0, L_1] \times [0, L_2] \times [0, N]$. In the following, for the sake of clarity, we use the notation $K_q(i, j, n) = K_q(i\Delta x_1, j\Delta x_2, n\Delta t)$ with $q \in \{\rho_0, v_0^1, p_0\}$.

### 3.3.2 Numerical scheme for solving the adjoint equations

The numerical schemes to solve the adjoint equations are the same schemes used for the state equations (see Appendix B2). However, special attention must be paid to the advection terms. In adjoint equations, these terms are of the form $-v_0^1 \partial_{x_1}$. For these terms, the approximation of the derivative with respect to $x_1$ is performed according to eq. (B5) if $v_0 > 0$ and according to eq. (B6) if $v_0 \le 0$.

### 3.3.3 Boundary conditions for solving the adjoint equations

The ground boundary condition is treated differently compared to the state problem. Information on the adjoint state at the ground boundary is limited: only $v^2$ is constrained at the interface with the ground. To compute the vertical derivatives evaluated at the first points of the grid just above the altitude $z = 0$ km (i.e. at $j = 1$ or $j = 1/2$), the fourth-order scheme requires the definition of ghost points. However, we do not have conditions to define a value for these points. Instead, we use a second-order scheme to compute the derivative with points belonging only to the domain. Thus, the derivatives of a quantity $q$ with respect to $x_2$ defined at the point $(i, 1/2)$ or at the point $(i, 1)$ are:

$$ \partial_{x_2} q(i, 1/2) = \frac{q_{i,1} - q_{i,0}}{\Delta x_2}, \qquad (18) $$

$$ \partial_{x_2} q(i, 1) = \frac{-3q_{i,1} + 4q_{i,2} - q_{i,3}}{2\Delta x_2}. \qquad (19) $$

For atmospheric boundary conditions, we use Auxiliary Differential Equation Perfectly Matched Layers (ADE-PML) based outer boundary conditions (Martin et al. 2010) to absorb efficiently the outgoing adjoint wavefield and mimic an open half-space atmosphere.

### 3.3.4 Numerical scheme for determining the sensitivity kernel

Let the model $\mathbf{m}$, the state variable $\mathbf{u}_\mathbf{m}$ such that $L(\mathbf{u}_\mathbf{m}; \mathbf{m}) = 0$ and the state variable $\lambda_\mathbf{m}$ satisfying eq. (9).

Gradients of the cost function are computed for $(\mathbf{m}, \mathbf{u}_\mathbf{m}, \lambda_\mathbf{m})$, according to eqs (13) to (16), as follows:

$$ \partial_{\rho_0 j} \tilde{\chi} = \Delta x_1 \Delta x_2 \Delta t \sum_i \sum_n K_{\rho_0}(i, j, n), \qquad (20) $$

$$ \partial_{v_0^j} \tilde{\chi} = \Delta x_1 \Delta x_2 \Delta t \sum_i \sum_n K_{v_0^1}(i, j, n), \qquad (21) $$

$$ \partial_{\rho_0 j} \tilde{\chi} = \Delta x_1 \Delta x_2 \Delta t \sum_i \sum_n K_{\rho_0}(i, j, n). \qquad (22) $$

In particular, we write $\partial_{q_j} \tilde{\chi} = \sum_i K_q^{2D}(i, j)$ for $q \in \{\rho_0, v_0^1, p_0\}$, with $K_q^{2D}$ defined such as:

$$ K_q^{2D}(i, j) = \Delta x_1 \Delta x_2 \Delta t \sum_n K_q(i, j, n). \qquad (23) $$

$K_q^{2D}$ corresponds to the sensitivity kernels of the infrasound waveforms to quantity $q$ for a 2-D model. This term will be specifically calculated in the Sections 3.4.2 and 3.5. The sensitivity kernels of a 2-D model are easier to interpret than those of a 1-D model because we can directly interpret it in terms of ray path of infrasound and its waveforms, whereas this information is integrated along the x-axis for the 1-D version.

Finally, the discretized model $\mathbf{m}$ is such that:

$$ \mathbf{m} = (\rho_{01}, \dots, \rho_{0L_2}, \mathbf{v}_0^1, \dots, \mathbf{v}_{0L_2}^1, p_{01}, \dots, p_{0L_2})^T \in \mathbb{R}^{3L_2}, \qquad (24) $$
and the derivative of $\tilde{\chi}$ with respect to $m$ is:

$$ \partial_m \tilde{\chi}(m, u_m, \lambda_m) = \left( \frac{\partial \tilde{\chi}}{\partial \rho_{01}}, \dots, \frac{\partial \tilde{\chi}}{\partial \rho_{0L_2}}, \frac{\partial \tilde{\chi}}{\partial v_{01}^1}, \dots, \frac{\partial \tilde{\chi}}{\partial v_{0L_2}^1}, \frac{\partial \tilde{\chi}}{\partial p_{01}}, \dots, \frac{\partial \tilde{\chi}}{\partial p_{0L_2}} \right)^T (m, u_m, \lambda_m) \in \mathbb{R}^{3L_2}. \quad (25) $$

### 3.3.5 Boundary conditions for computing the sensitivity kernel

The information in the atmospheric buffer areas of the PML is not relevant, due to the numerical and non-physical nature of the PML process in these absorbing outer layers. Hence, we define the sensitivity kernels within these zones by extension through invariant continuity. Besides, no specific treatment is required for the ground boundary.

### 3.3.6 Double checkpointing

Since the sensitivity kernels are obtained by a correlation calculation between the infrasound wavefield and the adjoint wavefield, it is necessary to know the solution of the forward problem backward-in-time, that is, from $T_{max}$ to 0 s.

The backward-in-time solution is not easy to obtain when memory resources are limited. It can even be excessively expensive in terms of computing resources and computing time, especially when the direct problem considers a large $\Omega$ propagation domain. However, it is important to speed up this calculation, as it is necessary for calculating sensitivity kernels, which are evaluated many times when solving the inverse problem.

Some techniques can be used to limit the number of calculations. Liu & Tromp (2006) and Tromp et al. (2008) use the backpropagation technique. However, computing the solution through the backpropagation process is not numerically stable in the presence of attenuation: due mainly to intrinsic non-causal and non-reversible principles, backpropagation will then amplify the wave amplitudes and the discretization errors in this case. Since PML are boundary layers that attenuate the wave, this technique is not suitable here. To overcome this numerical drawback, Komatitsch et al. (2016) developed the parsimonious storage technique, whereby they use part of the memory to limit the number of computations. However, this method still requires a lot of memory space, especially as infrasound simulations require finer grids than seismic wave simulations.

We thus propose here to use the principle of the parsimonious storage technique, while still limiting memory space and avoiding drastically increasing the number of calculations. For this, we define, implement and use a more innovative and optimized technique, called the Double checkpointing, further developed in Gerier (2024).

In this approach, we first solve the forward problem by making backups (see step 1 in Fig. 1) in the same way as parsimonious storage. Each backup is a checkpoint and can then be used as a starting point for finding the solution at a later time by solving the forward problem.

Next, we can consider the subproblem of getting the solution in the reverse time direction between the instants of two backups. To solve it, we can compute the forward problem between these two instants and make backups, known as ‘local’ checkpoints (see step 2 in Fig. 1). These so-called local checkpoints serve as a starting point for solving the forward problem up to the instant of interest (see step 3 in Fig. 1). Note that once local checkpoints are used they can be removed from storage step by step.

There are two levels of checkpoints in this method: global and local. We thus call this method double checkpointing. Compared to the single checkpointing method, this double checkpointing method has the advantage of being much faster while using a number of backups almost similar to that of the single checkpointing strategy.

## 3.4 Sensitivity kernel validation with small perturbation method

One way to validate the computation of the sensitivity kernels is through comparison with the sensitivity kernels acquired by the small perturbation method (presented in Nocedal & Wright 1999 as finite-difference derivative approximations). As a matter of fact, the sensitivity kernels of infrasound waveforms to the parameters correspond to the gradient of the misfit $\chi$ and, thus, the value of the sensitivity kernel associated with the parameter $m_i$ is the derivative of $\chi$ with respect to $m_i$, $\forall i \in [1, M]$.

The first derivative of $\chi$ with respect to $m_i$ is given by the small perturbation method, corresponding to a finite difference derivative calculated at the first order as follows:

$$ \frac{\partial \chi}{\partial m_i}(m) = \frac{\chi(m_1, \dots, m_{i-1}, m_i + \delta m_i, m_{i+1}, \dots, m_M) - \chi(m)}{\delta m_i}. \quad (26) $$

We compute the sensitivity kernel using the adjoint state method and the small perturbation method for the same configuration.

We propose to compare the sensitivity kernels of the two methods in 2-D to better understand and interpret the sensitivity of infrasound waveforms. The model $m$ defined in eq. (24) is replaced by

$$ m^{2D} = (\rho_{011}, \dots, \rho_{0ij}, \dots, \rho_{0L_1L_2}, \mathbf{v}_0^{11}, \dots, \mathbf{v}_0^{1ij}, \dots, \mathbf{v}_0^{1L_1L_2}, p_{011}, \dots, p_{0ij}, \dots, p_{0L_1L_2})^T. \text{ The number of model parameters is then } 3 \times L_1 \times L_2. \text{ The small perturbation method gives an approximation of } \frac{\partial \tilde{\chi}}{\partial m^{2D}}. \text{ The first } L_1L_2 \text{ components correspond to } K_{\rho_0}^{2D}. \text{ The next } L_1L_2 \text{ components correspond to } K_{v_0^1}^{2D} \text{ and the last components to } K_{\rho_0}^{2D}. $$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Illustration of the double checkpointing technique, with preliminary steps (1 and 2) and acquisition of the solution at time $t = n\Delta t$.</figcaption></figure>

### 3.4.1 Configuration

The configuration is defined by a 2-D model of size 20 km × 20 km. The model for which we obtain the observations is named ‘target model’, and is homogeneous, with a horizontal wind of 10 m.s⁻¹. The prior model is equal to the target model, except that (1) the wave speed is 10 per cent greater than that of the target model and (2) the horizontal wind is set to zero. The source is defined in time by the Gaussian first derivative wavelet with a dominant frequency of 0.1 Hz and it occurs at the point (4 km, 5 km). A receiver is located at (13 km, 5 km). The ground is a reflective surface, whereas the other boundaries are perfectly matched layers that are designed to absorb any outgoing wave signal.

### 3.4.2 Results

The sensitivity of the infrasound waveforms relative to the atmospheric model is calculated in the configuration described above, using the adjoint method and the small perturbation method.

Figs 2(b), (e) and (h) show the sensitivity kernels of infrasound waveforms on a 2-D grid respect to pressure, density and horizontal wind, respectively. Here, they are represented by the red and blue colours, which indicate where the waveform is sensitive to the model, and if, respectively, the parameters of the model should be increased or decreased to improve the misfit.

Two areas are particularly illuminated and correspond to the two different paths of infrasound (direct path and path through a reflection on the ground). The colour in these areas means whether we underestimate pressure or wind, or we overestimate density.

Moreover, in this simple example, we can see the shape of a banana doughnut, found in seismology (Tromp *et al.* 2014) for similar cases.

Concentric ellipses around the source and sensors can be seen in Figs 2(b), (e) and (h), with a lower intensity. These ellipses suggest that other ray trajectories can explain the observed wave trains. Indeed, a ray emitted at the source can be scattered at a point that belongs to the ellipse, and the ray can then reach the receiver at the time of the observation. These ellipses also appeared in the sensitivity kernels computed in seismology by Tromp *et al.* (2014).

Furthermore, on the one hand, as in Tromp *et al.* (2014) we can note that the scalar quantities have a sensitivity kernel that is symmetrical with respect to the perpendicular bisector of the source–receiver segment. On the other hand, vector quantities such as wind imply a non-symmetric sensitivity kernel (Fig. 2h) since these quantities are directional.

Although the sensitivity kernels are consistent with those of the seismology, we validate them with the small perturbation method whose Figs 2(a), (d) and (g) illustrate the result. Panels (a) and (b) are similar, as are panels (d) and (e) or panels (g) and (h), illustrating that both methods give the same shape of sensitivity kernel.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Comparison of the sensitivity kernels of infrasound waveforms acquired by different methods. Panels (a), (d) and (g) show the sensitivity kernels with respect to pressure, density and wind, respectively, and given by the small perturbation method. Panels (b), (e) and (h) show the sensitivity kernels with respect to pressure, density and wind, respectively, and given by the adjoint state method. Panels (c), (f) and (i) are the absolute differences of the sensitivity kernels with respect to pressure, density and wind, respectively. Colour scale is saturated at 10 per cent of the maximum sensitivity (not including the area around the source and receiver). Panel (j) shows the pressure signals of data and synthetics. Panel (k) shows the difference of time-series presented in panel (j).</figcaption></figure>

Figs 2(c), (f) and (i) complete the analysis showing the absolute difference of the sensitivity kernels acquired by the two methods.

The finite difference method and the adjoint state method give similar results, with a relative error generally lower than 10 per cent. However, the error is greater than 10 per cent for two cases. First, large errors are located around the source and receiver. We attribute this difference to the mathematical and physical singularity at the source and the receiver (the latter being also the point source of the adjoint wavefield). Indeed, the discretization of the direct problem and the small perturbation method to calculate the derivative of $\tilde{\chi}$ are approximation methods that do not allow the physics of the wave to be modelled adequately at these singularities. Besides, the error is greater than 10 per cent when the sensitivity kernel changes sign. We can explain these errors, since the small perturbation method is a first-order approximation in space.

## 3.5 Sensitivity kernel application in the case of a real explosion: the Hukkakero site case

Sensitivity kernels were computed in the previous section on small test cases. We extend the analysis here by computing the sensitivity kernels of infrasound waveforms in a case where the dimensions are closer to the application cases.

We consider here the explosion at the Hukkakero site on 2017/08/22 11:59:59. The IS37 station of the IMS network, located 320 km away from the Hukkakero site, recorded this explosion. As illustrated by fig. 3 in Vorobeva *et al.* (2023), the record shows three distinct arrivals. However, acoustic geometry predicts only the last two, which are thermospheric arrivals. Vorobeva *et al.* (2023) studied this case
and inverted the recordings for small-scale vertical variations due to gravity waves. We aim here at explaining the non-predicted arrival using the sensitivity kernels.

### 3.5.1 Configuration

We model this event using a domain of 460 km × 150 km size. The spatial step is 100 m. The time step is 0.05 s.

To model this event, we consider an idealized setting. Thus, although the network aperture is 2 km, we use six receivers, located 320–330 km away from the source. The source and the various receivers are assumed to be at an altitude of 50 m. The source time function is modelled by the first derivative of a Gaussian function in time, defined as follows:

$$S(t) = -2\pi^2 f_0^2 (t - t_0) \exp(-\pi^2 f_0^2 (t - t_0)^2). \qquad (27)$$

$f_0 = 0.1$ Hz is the dominant frequency of the source and is chosen so as to limit the calculation time, and $t_0 = \frac{1}{f_0}$ is a time-shift.

We define a prior model: we use the SD-WACCM-X model (Gasperini 2019), also used in the Vorobeva et al. (2023) study.

Observations are also idealized in an effort to validate the sensitivity kernels. Therefore, we do not use the IS37 station records, but we create observations based on the actual records. To do this, two different approaches are adopted.

In the first approach, we model the wave propagation with the prior model. To the waveforms of each receiver, we add a Gaussian derivative (whose equation is proportional to eq. 27). This perturbation is positioned in time to coincide with the first wave train recorded by the IS37 station. Thus, the six stations present a signal around 18 min. after the explosion, with the same apparent velocity as suggested in the IS37 recordings.

For the second approach, the horizontal wind of the SD-WACCM-X model is directly based on the results of Vorobeva et al. (2021). To obtain the set of observations, we simulate the wave propagation for a modified version of this model. We add a perturbation, defined as follows:

$$v_0^1(\mathbf{x}) = v_0^1_{\text{WACCM}}(\mathbf{x}) + 0.5 \left( 1 - \cos \left( 2\pi \frac{(h + x_2 - H)}{\lambda} \right) \right) \delta v_0^1 \text{ if } |x_2 - H| < h, \qquad (28)$$

with $H = 60$ km, the mean altitude of the perturbation, $h = 6$ km, the maximal range of the perturbation, $\delta v_0^1 = 12$ m.s⁻¹, $\lambda = 12$ km, the parameters to define the intensity of the perturbation and $v_0^1_{\text{WACCM}}$, the extracted wind model of SD-WACCM-X on 2017/08/22 12:00:00.

Moreover, for the second approach, a mask is added in the misfit and in the adjoint source definition to select the differences between the ‘target’ observations and the synthetics at a specific time. Using the mask function, we focus either (1) on the first arrival only, (2) on the thermospheric arrivals only, or (3) on all the different arrivals.

Although the effects of gravity in this simulation are negligible (no gravity waves or gravito-acoustic wave interactions), gravity is not ignored in our simulation tool. Our tool should be able to handle both gravity and acoustic waves as well as their interactions. Although we did not validate these features, we ensured that our results are consistent with acoustic simulations performed without gravity terms in the equations. Fig. S1 in the Supplementary Material presents a validation case in the absence of gravity.

### 3.5.2 Results

Figs 3 and 4, respectively, show the sensitivity kernels of infrasound waveforms related to wave speed and wind perturbation for the first approach. Figs 5–7 correspond to the sensitivity kernels of infrasound waveforms to wind perturbation for the second approach, with the different masks applied to the misfit.

First, Figs 3 and 4 clearly identify the area where the first arrival in the waveforms is sensitive to the model. Waveforms are particularly sensitive where the wave is refracted back to the ground. Sensitivity kernels predict that arrival at receivers at around 18 min is due to refraction at 60 km altitude.

Note that the model presents sudden variations in wave speed at 10 km, and the sensitivity kernels underline these variations. Some wave fronts could have been ducted by the tropospheric waveguide and then diffracted into stratospheric waves.

Moreover, we can note that the sensitivity kernels of infrasound related to wind perturbation and wave speed perturbation are similar, even though the perturbation in the model was only added to the wind. Due to the geometry of the problem (receivers are located only to the right of the source), modifying the wind or the wave speed has the same effect on the effective wave speed, and thus on the infrasound behaviour. This tells us that the sensitivity kernels do not inform us about the most relevant way to modify the model in terms of winds or wave speed. They only point in the direction of the misfit improvement. Thus, it is possible to improve the misfit by also playing with the wave speed in that example.

Fig. 5, which only considers the difference between observed and synthetic first arrivals of the first wave trains, presents results similar to Fig. 4. The sensitivity kernel of Fig. 5 has the form of an arc twice wider than that of Fig. 4. The one in Fig. 4 is similar to the sensitivity kernel of arrival times and gives the exact position where the ray is refracted towards the ground. The size of the arc is directly determined by the frequency of the added perturbation, and so by the source frequency, as shown by Hung et al. (2000) for the seismic sensitivity kernels. The arc of the sensitivity kernel in Fig. 5 is larger because the difference between observation and synthetic (Fig. 5e) spans over a longer
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Sensitivity kernel of infrasound waveforms related to wave speed perturbations, when adding a perturbation to synthetics to reproduce IS37 recordings. Infrasound is modelled by Euler's equations in the presence of gravity. (a) Wave speed model. The solid red line shows the model used to get the 'target' observations. The dashed black line shows the prior model. (b) Sensitivity kernel related to wave speed perturbations. The grey line shows the ray path obtained with InfraGa. (c) Integral over $x_1$ of the sensitivity kernel shown in panel (b). (d) Effective wave speed model (same legend as panel a). (e) 'Target' and synthetic pressure time-series acquired at each receiver. Each receiver is identified by its horizontal distance from the source. (f) Difference of time-series ('Target' data - synthetic pressure) presented in the panel (e).</figcaption></figure>

time interval than that in Fig. 4(e). Consequently, the area of the atmosphere affecting this part of the infrasound waveform is also larger. The duration of the difference between observation and synthetic data depends on the evolution of the wave, potential interferences, etc.

Focusing on the last arrivals in the recordings completely changes the sensitivity kernels, as shown in Fig. 6. We can observe the path of the thermospheric rays and particularly the altitude at which they are refracted towards the ground (around 120 and 140 km). The wind model perturbation added in the model does not just refract waves towards the ground, it modifies the waveform of the waves passing through this perturbation, in particular by accelerating the thermospheric waves in this case. As we can see in Fig. 6(d), thermospheric rays are effectively modified. Moreover, the integral over $x_1$ of the sensitivity kernels in Fig. 6(c) presents the sensitivity of infrasound waveforms to the model perturbation. Therefore, waveforms are particularly sensitive where thermosphere rays are refracted towards the ground. From 0 to 140 km, sensitivity is non-zero. This shows that the model can be modified between these altitudes to reduce the discrepancy between synthetics and observations.

Fig. 7 shows the sensitivity kernels when considering all the recordings. These sensitivity kernels are similar to the sum of those shown in Figs 5 and 6. Nevertheless, the part related to thermosphere rays is more pronounced. This shows that the thermospheric arrivals dominate the misfit because of their amplitude.

In this example, the waveform in Fig. 7(d) is simple because it comes from an idealized pulsed signal and the observations come from a relatively smooth model. Consequently, the error on the waveform for the last two wave trains in Fig. 7(e) can be assimilated to a traveltime error. On the other hand, the first wave train is not explained by ray theory, as shown by Fig. 7(b). Hence, we cannot compute the sensitivity of the arrival time of this wave train with respect to the model.

As part of this exercise, we used a network of six sensors. However, according to additional tests, similar results are obtained when a single sensor is used. This is because the disturbance present in the data has physical meaning and accurately reflects the traveltime and the angle of incidence at each station. However, the use of a network of sensors makes sense for deciphering between waves with different incidence angles and for removing noise.

The sensitivity kernels in panel (b) of Figs 3–7 clearly identify the path of each wave packet, whereas the sensitivity kernels in panel (c) of Figs 3–7 highlight only where the rays are refracted. Nevertheless, the 1-D version also measures where and how the vertical profile of the model is sensitive to the infrasound. Thanks to these characteristics, sensitivity kernels are similar to the vertical weights used in the data assimilation framework as defined in Amezcua *et al.* (2020).
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Sensitivity kernel of infrasound waveforms related to wind perturbations, when adding a perturbation to synthetics to reproduce IS37 recordings. Infrasound is modelled by Euler's equations in the presence of gravity. (a) Wind model. The solid red line shows the model used to get the 'target' observations. The dashed black line shows the prior model. (b) Sensitivity kernel related to wind perturbations. The grey line shows the ray path obtained with InfraGa. (c) Integral over $x_1$ of the sensitivity kernel shown in panel (b). (d) Effective wave speed model (same legend as panel a). (e) 'Target' and synthetic pressure time-series acquired at each receiver. Each receiver is identified by its horizontal distance from the source. (f) Difference of time-series ('Target' data - synthetic pressure) presented in the panel (e).</figcaption></figure>

# 4 TOWARDS THE INVERSE PROBLEM USING THE SENSITIVITY KERNELS

## 4.1 Inversion methods

Here we consider the problem (7), which is equivalent to the problem (6). The problem (7) implies that the cost function and its derivative with respect to $m$ are evaluated at $(m, u_m, \lambda_m)$. For simplicity, we denote the cost function as $\tilde{\chi}(m)$ to refer to $\tilde{\chi}(m, u_m, \lambda_m)$ in the following. We explain here the methodology for minimizing the cost function $\tilde{\chi}$ as a function of $m$ using its gradient $\partial_m \tilde{\chi}$, also called the sensitivity kernel.

Since the sensitivity kernel expresses where and how the model should be modified so that the forward modelling corresponds to the observed data, we update iteratively the current model as follows

$$m^{k+1} = m^k + \alpha^k r^k, \qquad (29)$$

where $r^k$ is the direction of improvement of the model, depending on the sensitivity kernel, and $\alpha^k$ is a step length.

We adopt the 'line search' method (Nocedal & Wright 1999) (as opposed to the 'trust-region' method) to improve the model at each iteration. The line search method consists of first identifying the direction of descent and then the step length.

### 4.1.1 Finding the descent direction

Descent direction can be computed in various ways. After test with different descent directions (see Section 2.2 of Supplementary material), we chose to use the Limited-BFGS method, a Quasi-Newton method suitable for nonlinear and large optimization problems. We implemented and used the L-BFGS algorithm 7.4 provided in Nocedal & Wright (1999).

However, the L-BFGS method requires the gradients of the previous iterations (the previous 3 to 5 iterations). Therefore, the L-BFGS algorithm requires an initialization phase. In our case, we initialize the L-BFGS method on the first 5 iterations, where the descent direction is computed using the Steepest Descent direction. Then, starting from the fifth iteration, the descent direction is calculated using the L-BFGS method.

### 4.1.2 Finding the step length

We aim here to find the step length $\alpha^k$, that is, the distance between the current model $m^k$ and the next model $m^{k+1}$.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Sensitivity kernel of infrasound waveforms related to wind perturbations, when adding a perturbation in the wind model, and considering only the first arrival in the observations. Infrasound is modelled by Euler's equations in the presence of gravity. (a) Wind model. The solid red line shows the model used to get the 'target' observations. The dashed black line shows the prior model. (b) Sensitivity kernel related to wind perturbations. The grey line shows the ray path obtained with InfraGa. (c) Integral over $x_1$ of the sensitivity kernel shown in panel (b). (d) Effective wave speed model (same legend as panel a). (e) 'Target' and synthetic pressure time-series acquired at each receiver. Each receiver is identified by its horizontal distance from the source. (f) Difference of time series ('Target' data - synthetic pressure) presented in the panel (e). The hatch part in panels (e) and (f) indicates where the mask on the misfit is applied: only the first wave train arrival is considered here.</figcaption></figure>

The step length $\alpha^k$ must be chosen carefully to avoid too many iterations. In this purpose, $\alpha^k$ is selected using the two criteria of the strong Wolfe conditions (Nocedal & Wright 1999) (i.e. the minimization of the misfit function and its curvature). To select $\alpha^k$, we use Algorithm 3.5 defined in Nocedal & Wright (1999).

As for the descent direction, we apply a different method to choose the step lengths in the first five iterations. For this, we use the backtracking algorithm, as defined in Nocedal & Wright (1999). This algorithm ensures a sufficient decrease in the cost function for the first iterations.

Finally, the inversion algorithm is presented in Algorithm 1 in Appendix C3.

### 4.1.3 Regularization and choice of the optimization algorithm

A regularization term can be added to the cost function $\tilde{\chi}$ to regularize the inverse problem and defined such that:

$$ \psi = \frac{\|m - m^0\|^2}{\|m^0\|^2}. \qquad (30) $$

Assuming that the prior model is fairly close to reality, this term is used to restrict the search space for solutions around the prior model.

The regularization term is similar to the second term of the misfit, defined by Assink (2012) in his part 4.3, where the covariance matrix is proportional to the identity matrix. Other prior information can be used to assess the weight of the different model parameters.

In our case, we choose $\psi = 0$, since the model is sufficiently well constrained thanks to the different receivers.

### 4.1.4 Definition of the test case

Before performing inversions in 2-D cases, we first ran two inversion cases for 1-D acoustic propagation (see Section S2 of Supplementary Material). These initial tests allowed us to choose the most appropriate hyperparameters for the inversion method, such as the descent direction of the L-BFGS method mentioned in Section 4.1.1, but also to validate that the inversion process works on simpler cases. Now, we present here idealized inversions in 2-D cases with realistic atmospheric models.

We aim to model the effective wave speed, defined by $c_0 + v_0^1$, up to 50 km altitude.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Sensitivity kernel of infrasound waveforms related to wind perturbations, when adding a perturbation in the wind model and considering only the end of the observations. Infrasound is modelled by Euler's equations in the presence of gravity. (a) Wind model. The solid grey line shows the model used to get the 'target' observations. The dashed black line shows the prior model. (b) Sensitivity kernel to wind perturbations. The grey line shows the ray path obtained with InfraGa. (c) Integral over $x_1$ of the sensitivity kernel shown in panel (b). (d) Effective wave speed model (same legend as panel a). (e) 'Target' and synthetic pressure time-series acquired at each receiver. Each receiver is identified by its horizontal distance from the source. (f) Difference of time-series ('Target' data - synthetic pressure) presented in the panel (e). The hatch part in panels (e) and (f) indicates where the mask on the misfit is applied: only the last wave train arrivals are considered here.</figcaption></figure>

To do so, we want to retrieve the tropospheric guide formed by wind variations and defined as

$$v_0^1(x) = A_{\text{wind}} \exp \left( - \frac{(x_2 - H)^2}{2\sigma^2} \right), \qquad (31)$$

where $A_{\text{wind}} = 45 \text{ m.s}^{-1}$ is the maximal amplitude of the wind, $H = 10 \text{ km}$ is the altitude, where the wind variation is maximal, and $\sigma^2 = 1.5 \text{ km}^2$. The other quantities (density, pressure) are derived from the reference model 'ToyAtmo.met' of the InfraGA software (Blom & Waxler 2012; Blom 2014). Figs 8(c)–(e) present the target model in pink.

We propose to start the inversion process with the vertical profile, shown in Figs 8(c)–(e) in yellow. The horizontal wind is zero at each altitude. The density quantity is assumed to be known (equal to that of the target model). The wave speed is well known under 10 km of altitude. A similar test is presented in Section S3 (Supplementary Material) where the prior wave speed is equal to the target wave speed.

In addition, a source emits infrasound at 250 m altitude. The source time function is a Gaussian of 0.1 Hz dominant frequency. Fig. 8(a) shows the position of the 34 receivers used. These receivers are divided into four sets, defined as follows:

(i) Ten receivers linearly spaced from (7, 45) to (92, 45) km.

(ii) Ten receivers linearly spaced from (7, 0.2) to (92, 0.2) km.

(iii) Seven receivers linearly spaced from (10, 2) to (10, 44) km.

(iv) Seven receivers linearly spaced from (90, 2) to (90, 44) km.

Fig. 8(b) shows the target recordings with the pink dashed line and the one computed with the prior model in yellow at the stations R01, R04, R07, R11, R14 and R17 located at (10, 2) km, (10, 23) km, (10, 44) km, (90, 2) km, (90, 23) km and (90, 44) km, respectively. The total dimension of the domain is (99, 57.6) km. The spatial step is 100 m, and the time step is 0.05 s.

The inversion algorithm selected to solve the 2-D test case uses the L-BFGS method to compute the descent direction and $(\ln(\rho_0), \frac{v_0^1}{100}, \frac{c_0}{100})$ as parametrization. The scaling factor 100 applied to the wave speed and wind in the previous parametrization has been updated compared to the one presented in Section 3.2 and the 1-D case (see Section S2 of Supplementary material). Empirically, we add a scaling factor so that the relative variations in the logarithm of density and wave speed are more balanced. The same factor is applied to wave speed and wind so as not to introduce a bias between these two variables.
<figure><img src="image_8.png" /></figure>

**Figure 7.** Sensitivity kernel of infrasound waveforms related to wind perturbations. Infrasound is modelled by Euler's equations in the presence of gravity. (a) Wind model. The solid red line shows the model used to get the 'target' observations. The dashed black line shows the prior model. (b) Sensitivity kernel related to wind perturbations. The grey line shows the ray path obtained with InfraGa. (c) Integral over $x_1$ of the sensitivity kernel shown in panel (b). (d) Effective wave speed model (same legend as panel a). (e) 'Target' and synthetic pressure time-series acquired at each receiver. Each receiver is identified by its horizontal distance from the source. (f) Difference of time-series ('Target' data—synthetic pressure) presented in the panel (e). In the current case, no mask is applied on the misfit.

<figure><img src="image_9.png" /></figure>

**Figure 8.** 2-D inversion test case. Panel (a) shows the source and receivers locations. Panel (b) shows the infrasound waveform observed at stations R11, R14, R17, R28, R31 and R34, located at (10, 2) km, (10, 23) km, (10, 44) km, (90, 2) km, (90, 23) km, and (90, 44) km, respectively. Panel (c) presents the evolution of cost function and relative distance between the current model and the target model. Panels (d), (e) and (f) show the wave speed, horizontal wind and the effective wave speed in function of altitude, respectively.
### 4.1.5 2-D inversion results

The model acquired after inversion is presented in Figs 8(d), (e) and (f) and the evolution of the cost function during the inversion procedure is illustrated in Fig. 8(c).

After 44 iterations, the vertical profile of the effective wave speed is largely reconstructed. The profile of the wave speed from 20 to 50 km has been corrected. The position of the wind jet is well identified at 10 km. Moreover, the error in the amplitude of the wind jet is less than 8 m.s⁻¹ and the error in terms of effective wave speed is 0.64 m.s⁻¹. These orders of magnitude of these errors are satisfactory. However, it also means that after 44 iterations, the algorithm has failed to fully distinguish between the contribution of wind speed and sound. However, the final model has some oscillations that are difficult to eliminate because of the peaks of sensitivity at the receiver locations. The oscillation can be explained since it occurs at the location where the gradients of wave speed and of the wind model are opposite and there are not enough receivers in this difficult area to understand the model variations at this specific location. In addition, as being expected the model has not been improved after 45 km of altitude since no receiver illuminates this area.

The cost function, which evaluates the distance between the observations and the synthetic pressure, has decreased to less than 0.1 and no difference is visible between the time-series presented in Fig. 8(b). The evolution of the relative distance of the model solution and the target model seems stable. For wind and wave speed model, the first iteration significantly reduces the distance between the current model and the target model. However, oscillations in the current model are observed. For subsequent iterations of the inversion, the current model is modified but the oscillations persist and do not allow this metric measuring the distance between the current and target models to be further improved.

To improve efficiently these first results and reduce the persistent oscillations in Figs 8(d)–(f), we can perform a multitage approach (Groos et al. 2017; Beraus et al. 2025). The inversion is done in several stages, where the source time function is filtered at different frequency ranges from the lowest to the highest frequencies. Starting the inversion with low frequencies allows to find a first estimate, and finishing the inversion with the higher frequencies allows then to retrieve fine-scale details for the estimated model. Then, the uncertainty on the inverted model is reduced.

For this example, we used a consequent number of receivers, which is far from real conditions. Here, it is a proof of concept of the algorithm at realistic scales and frequencies. Several numerical tests on receiver positions have shown that the location of receivers is important for this test. The vertical lines of the receivers allow to identify the altitude of the disturbance, and the receivers should be positioned to the left and right of the source to distinguish between wind and propagation speed. Other tests should be done to quantify the importance of the receiver locations and identify potential requirements.

## 5 CONCLUSIONS

Here we presented the first sensitivity kernels for infrasound waveforms. These sensitivity kernels, computed using the adjoint method, were first validated and then used on a full-scale case, revealing the sensitivity of infrasound waveforms produced by an explosion at the Hukkakero site in Finland, and recorded by ground sensors located ≈320 km away. For direct rays, the sensitivity is significant around the infrasound ray paths and particularly strong at altitudes of infrasound refraction towards the ground. For acoustic energy packets that cannot be explained by ray theory, the sensitivity kernels allow us to identify regions where scattered waves are refracted to the ground. The kernels have a complex shape due to potential scattering either in the troposphere or the stratosphere. Furthermore, we developed an inversion algorithm based on the sensitivity kernels and the Quasi-Newton method (Limited-BFGS), the latter has provided the best results among those obtained with several inversion methods we have tested in 1-D and 2-D (Steepest Descent, different versions of the nonlinear Conjugate Gradient, Quasi-Newton or BFGS methods). The algorithm has been run in 2-D synthetic cases and is able to recover sound speed and wind profiles.

Further studies should focus on the application of sensitivity kernels to infrasound recorded by stations. In particular, future work will look at the effect of data noise and attenuation effects on sensitivity kernels and inversion results. Additional work is also needed to test the inverse problem on more realistic and/or data assimilation examples, and with a more restricted number of receivers. In this context, the study of a parametrization resulting from the projection of the model onto a space of functions (such as the empirical orthogonal functions) to describe physical quantities could also be considered to reduce the number of parameters to be inverted.

## ACKNOWLEDGMENTS

The authors thank CALMIP for the considerable computing resources allocated to the project P1404, which enabled us to carry out the forward and inverse simulations. We thank the CNES for funding this study through scientific APR projects BIS/Strateole2 and the Région Occitanie (France) for partly funding the S. Gérier’s PhD thesis. We acknowledge International Monitoring System of CTBT Organisation for providing access to their data base through VDEC contract “ISAE-SUPAERO 2021-CSIF-129”. We are grateful to Aerospace Valley for their support. We thank a lot Sven Peter Näsholm and Phil Blom, for their deep and fruitful review of this study.

## SUPPORTING INFORMATION

Supplementary data are available at *GJIRAS* online.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.
The 3-hourly SD-WACCM-X model is available on https://www.earthsystemgrid.org/dataset/ucar.cgd.ccsm4.SD-WACCM-X.v2.1.atm.hist.3hourly_inst.html (Gasperini 2019). The InfraGA infrasound propagation code and ToyAtmo.met model (Blom & Waxler 2012; Blom 2014) are accessible on GitHub https://github.com/LANL-Seismoacoustics/infraGA.

REFERENCES

Amezcua, J., Näsholm, S.P., Blixt, E.M. & Charlton-Perez, A.J., 2020. Assimilation of atmospheric infrasound data to constrain tropospheric and stratospheric winds, *Q. J. R. Meteorol. Soc.*, **146**(731), 2634–2653.

Amezcua, J., Näsholm, S.P. & Rodriguez, I.V., 2024. Using satellite data assimilation techniques to combine infrasound observations and a Full Ray-Tracing Model to constrain stratospheric variables, *Monthly Weather Rev.*, **152**(8), 1883–1902.

Arrowsmith, S.J., Marcillo, O. & Drob, D.P., 2013. A framework for estimating stratospheric wind speeds from unknown sources and application to the 2010 December 25 bolide, *Geophys. J. Int.*, **195**(1), 491–503.

Assink, J., 2012. Infrasound as Upper Atmospheric Monitor, *Ph.D. thesis*, University of Mississippi.

Assink, J., Smets, P., Marcillo, O., Weemstra, C., Lalande, J.M., Waxler, R. & Evers, L., 2019. Advances in Infrasonic Remote Sensing Methods, *Infrasound Monitoring for Atmospheric Studies*, pp. 605–632, eds, Le Pichon, A., Blanc, E. & Hauchecorne, A., Springer International Publishing.

Assink, J.D., Waxler, R., Frazier, W.G. & Lonzaga, J., 2013. The estimation of upper atmospheric wind model updates from infrasound data, *J. Geophys. Res.: Atmospheres*, **118**(19).

Assink, J.D., Pichon, A.L., Blanc, E., Kallel, M. & Khemiri, L., 2014. Evaluation of wind and temperature profiles from ECMWF analysis on two hemispheres using volcanic infrasound, *J. Geophys. Res.: Atmospheres*, **119**(14), 8659–8683.

Aulanier, F., Nicolas, B., Roux, P. & Mars, J.I., 2013. Time-angle sensitivity kernels for sound-speed perturbations in a shallow ocean, *Acoust. Soc. Am. J.*, **134**(1), 88.

Batchelor, G.K., 2000. *An Introduction to Fluid Dynamics*, Cambridge Mathematical Library, Cambridge University Press, Cambridge.

Beraus, S., Köhn, D., Bohlen, T., Burschil, T., Schuster, B., Buness, H. & Gabriel, G., 2025. Seismic crosshole full-waveform inversion of high-frequency SV-waves for glacial sediment characterization, *Geophys. Prospect.*, **73**(5), 1587–1605.

Blanc, E. et al., 2018. Toward an improved representation of middle atmospheric dynamics thanks to the ARISE project, *Surv. Geophys.*, **39**(2), 171–225.

Blom, P., 2014. LANL-Seismoacoustics/infraGA, original-date: 2020-07-17 15:31:47Z.

Blom, P. & Waxler, R., 2012. Impulse propagation in the nocturnal boundary layer: analysis of the geometric component, *J. acoust. Soc. Am.*, **131**(5), 3680–3690.

Blom, P.S. & Marcillo, O.E., 2017. An optimal parametrization framework for infrasonic tomography of the stratospheric winds using non-local sources, *Geophys. J. Int.*, **208**(3), 1557–1566.

Brissaud, Q., Martin, R., Garcia, R.F. & Komatitsch, D., 2016. Finite-difference numerical modelling of gravitoacoustic wave propagation in a windy and attenuating atmosphere, *Geophys. J. Int.*, **206**(1), 308–327.

Brissaud, Q., Martin, R., Garcia, R.F. & Komatitsch, D., 2017. Hybrid Galerkin numerical modelling of elastodynamics and compressible Navier-Stokes couplings: applications to seismo-gravito acoustic waves, *Geophys. J. Int.*, **210**(2), 1047–1069.

Campus, P. & Christie, D.R., 2009. Worldwide observations of infrasonic waves, *Infrasound Monitoring for Atmospheric Studies*, pp. 185–234, eds, Le Pichon, A., Blanc, E. & Hauchecorne, A., Springer Netherlands.

Charlton, A.J., O’neill, A., Lahoz, W.A. & Massacand, A.C., 2004. Sensitivity of tropospheric forecasts to stratospheric initial conditions, *Q. J. R. Meteorol. Soc.*, **130**(600), 1771–1792.

Chunchuzov, I., Kulichkov, S., Perepelkin, V., Popov, O., Firstov, P., Assink, J.D. & Marchetti, E., 2015. Study of the wind velocity-layered structure in the stratosphere, mesosphere, and lower thermosphere by using infrasound probing of the atmosphere, *J. Geophys. Res.: Atmospheres*, **120**(17), 8828–8840.

Chunchuzov, I.P., Kulichkov, S.N. & Firstov, P.P., 2013. On acoustic N-wave reflections from atmospheric layered inhomogeneities, *Izvestiya, Atmospheric and Oceanic Physics*, **49**(3), 258–270.

Cornuelle, B.D. & Skarsoulis, E.K., 2022. Peak-time sensitivity kernels for noise cross-correlation envelopes, *Acoust. Soc. Am. J.*, **151**(4), 2353–2366.

Donn, W.L. & Rind, D., 1971. Natural Infrasound as an Atmospheric Probe, *Geophys. J. Int.*, **26**(1–4), 111–133.

Drob, D., 2019. Meteorology, climatology, and upper atmospheric composition for infrasound propagation modeling, *Infrasound Monitoring for Atmospheric Studies: Challenges in Middle Atmosphere Dynamics and Societal Benefits*, pp. 485–508, eds, Le Pichon, A., Blanc, E. & Hauchecorne, A., Springer International Publishing.

Drob, D.P., Picone, J.M. & Garcés, M., 2003. Global morphology of infrasound propagation, *J. Geophys. Res.: Atmospheres*, **108**.

Drob, D.P., Meier, R.R., Picone, J.M. & Garcés, M.M., 2009. Inversion of infrasound signals for passive atmospheric remote sensing, *Infrasound Monitoring for Atmospheric Studies*, pp. 701–731, eds, Le Pichon, A., Blanc, E. & Hauchecorne, A., Springer Netherlands.

Eyre, J.R., Bell, W., Cotton, J., English, S.J., Forsythe, M., Healy, S.B. & Pavelin, E.G., 2022. Assimilation of satellite data in numerical weather prediction. Part II: recent years, *Q. J. R. Meteorol. Soc.*, **148**(743), 521–556.

Fichtner, A., 2011. Full waveform tomography on continental scales, *Full Seismic Waveform Modelling and Inversion*, pp. 233–265, ed., Fichtner, A., Springer.

Fichtner, A., Bunge, H.P. & Igel, H., 2006. The adjoint method in seismology: I. Theory, *Phys. Earth planet. Inter.*, **157**(1–2), 86–104.

Fricke, J.T., Evers, L.G., Smets, P.S.M., Wapenaar, K. & Simons, D.G., 2014. Infrasonic interferometry applied to microbaroms observed at the Large Aperture Infrasound Array in the Netherlands, *J. Geophys. Res.: Atmospheres*, **119**(16), 9654–9665.

Froment, M., 2023. *Meteorite Impacts on Mars: Seismic Observation, Source Theory and Modeling.*, These de doctorat, Université Paris Cité.

Froment, M. et al., 2024. Inferring the speed of sound and wind in the nighttime martian boundary layer from impact-generated infrasound, *Geophys. Res. Lett.*, **51**(18), e2024GL109726.

Garcés, M., Willis, M., Hetzer, C., Le Pichon, A. & Drob, D., 2004. On using ocean swells for continuous infrasonic measurements of winds and temperature in the lower, middle, and upper atmosphere, *Geophys. Res. Lett.*, **31**(19).

Garcia, R.F., Brissaud, Q., Rolland, L., Martin, R., Komatitsch, D., Spiga, A., Lognonné, P. & Banerdt, B., 2017. Finite-difference modeling of acoustic and gravity wave propagation in Mars atmosphere: application to infrasounds emitted by meteor impacts, *Space Sci. Rev.*, **211**(1–4), 547–570.

Gasperini, F., 2019. 1980-2017 specified dynamics (SD) WACCM-X version v2.1 simulation with nudging of MERRA-2 data, [Dataset]. UCAR/NCAR Climate Data Gateway. https://doi.org/10.26024/5b58-nc53.

Gauthier, O., 1986. Two-dimensional nonlinear inversion of seismic waveforms: numerical results, *Geophysics*, **51**(7), 1387.

Gerber, E.P., Orbe, C. & Polvani, L.M., 2009. Stratospheric influence on the tropospheric circulation revealed by idealized ensemble forecasts, *Geophys. Res. Lett.*, **36**(24).
Gerier, S., 2024. *Étude et caractérisation de la haute atmosphère au moyen des ondes infrasonores acquises par des vols de ballons stratosphériques, et par inversion numérique*, Ph.D. thesis, Thèse de doctorat dirigée par Garcia, Raphael.

Groos, L., Schäfer, M., Forbriger, T. & Bohlen, T., 2017. Application of a complete workflow for 2D elastic full-waveform inversion to recorded shallow-seismic Rayleigh waves, *Geophysics* **82**(2), R109–R117.

Hermand, J.P., 2006. Adjoint-based acoustic inversion for the physical characterization of a shallow water environment, *Acoust. Soc. Am. J.* **119**(6), 3860.

Hung, S.H., Dahlen, F. & Nolet, G., 2000. Fréchet kernels for finite-frequency traveltimeséII. Examples, *Geophys. J. Int.*, **141**(1), 175–203.

Köhn, D., De Nil, D., Kurzmann, A., Przebindowska, A. & Bohlen, T., 2012. On the influence of model parametrization in elastic full waveform tomography, *Geophys. J. Int.*, **191**(1), 325–345.

Komatitsch, D., Xie, Z., Bozdağ, E., Sales de Andrade, E., Peter, D., Liu, Q. & Tromp, J., 2016. Anelastic sensitivity kernels with parsimonious storage for adjoint tomography and full waveform inversion, *Geophys. J. Int.*, **206**(3), 1467–1478.

Lalande, J.M., Sèbe, O., Landès, M., Blanc-Benon, P., Matoza, R.S., Le Pichon, A. & Blanc, E., 2012. Infrasound data inversion for atmospheric sounding, *Geophys. J. Int.*, **190**(1), 687–701.

Laloyaux, P., Kurth, T., Dueben, P.D. & Hall, D., 2022. Deep learning to estimate model biases in an operational NWP assimilation system, *J. Adv. Model. Earth Syst.*, **14**(6), e2022MS003016.

Landau, L. & Lifshitz, E., 1987. *Fluid Mechanics, Vol. 6*, Butterworth-Heinemann.

Le Pichon, A., Blanc, E. & Drob, D., 2005. Probing high-altitude winds using infrasound, *J. Geophys. Res.: Atmospheres*, **110**.

Le Pichon, A., Ceranna, L., Garcés, M., Drob, D. & Millet, C., 2006. On using infrasound from interacting ocean swells for global continuous measurements of winds and temperature in the stratosphere, *J. Geophys. Res.: Atmospheres*, **111**.

Lemke, M. & Sesterhenn, J., 2016. Adjoint-based pressure determination from PIV data in compressible flows—Validation and assessment based on synthetic data, *Eur. J. Mech. - B/Fluids*, **58**, 29–38.

Liu, Q. & Tromp, J., 2006. Finite-frequency kernels based on adjoint methods, *Bull. seism. Soc. Am.* **96**(6), 2383–2397.

Liu, Q. & Tromp, J., 2008. Finite-frequency sensitivity kernels for global seismic wave propagation based upon adjoint methods, *Geophys. J. Int.*, **174**(1), 265–286.

Martin, R. & Coudier-Castaneda, C., 2010. An improved unsplit convolutionally perfectly matched layer absorbing technique for the Navier-Stokes equations using cut-off frequency shift, *Comput. Model. Eng. Sci. (CMES)*, **63**(1), 47–77.

Martin, R., Komatitsch, D., Gedney, S. & Bruthiaux, E., 2010. A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using auxiliary differential equations (ADE-PML), *Comput. Model. Eng. Sci.*, **56**(1), 17–42.

Martire, L., Martin, R., Brissaud, Q. & Garcia, R.F., 2022. SPECFEM2DDG, an open-source software modelling mechanical waves in coupled solid-fluid systems: the linearized Navier-Stokes approach, *Geophys. J. Int.*, **228**(1), 664–697.

Marty, J., 2019. The IMS Infrasound Network: current status and technological developments, *Infrasound Monitoring for Atmospheric Studies: Challenges in Middle Atmosphere Dynamics and Societal Benefits*, pp. 3–62, eds, Le Pichon, A., Blanc, E. & Hauchecorne, A., Springer International Publishing.

Monteiller, V., Chevrot, S., Komatitsch, D. & Wang, Y., 2015. Three-dimensional full waveform inversion of short-period teleseismic wavefields based upon the SEM-DSM hybrid method, *Geophys. J. Int.*, **202**(2), 811–827.

Nocedal, J. & Wright, S.J., 1999. *Numerical Optimization*, Springer.

Peter, D. et al., 2011. Forward and adjoint simulations of seismic wave propagation on fully unstructured hexahedral meshes, *Geophys. J. Int.*, **186**(2), 721–739.

Rieutord, M., 2016. *Fluid Dynamics: An Introduction, Graduate Texts in Physics*, Springer International Publishing.

Schafer, M., Groos, L., Forbriger, T. & Bohlen, T., 2014. Line-source simulation for shallow-seismic data. Part 2: full-waveform inversion-a synthetic 2-D case study, *Geophys. J. Int.* **198**, 1405–1418.

Simmons, A. et al., 2020. Global stratospheric temperature bias and other stratospheric aspects of ERA5 and ERA5.1. https://www.researchgate.net/profile/L-Haimberger/publication/339442429Global_stratospheric temperature_bias_and_other_stratospheric_aspects_of_ERA5_and_ERA51/links/5e529d91a6fddcc2f8f5d6515/Global-stratospheric-temperature-bias-and-other-stratospheric-aspects-of-ERA5-and-ERA51.pdf.

Skarsoulis, E.K. & Cornuelle, B.D., 2004. Travel-time sensitivity kernels in ocean acoustic tomography, *Acoust. Soc. Am. J.*, **116**(1), 227–238.

Spieser, E. & Bailly, C., 2018. Adjoint solution to the linearized Euler equations computed from the direct problem. *CFA’18 LE HAVRE - 14ème Congrès Français d’Acoustique*, pp. 73–79, Le Havre, France. https://doi.org/10.1017/jfm.2020.469.

Spieser, E. & Bailly, C., 2020. Sound propagation using an adjoint-based method, *J. Fluid Mech.*, **900**, A5.

Spieser, E., Bodard, G. & Bailly, C., 2020. Jet noise prediction based on a self-adjoint operator. *e-Forum Acusticum 2020*, pp. 1325–1332, Lyon, France.

Stein, L., Straube, F., Sesterhenn, J., Weinzierl, S. & Lemke, M., 2019. Adjoint-based optimization of sound reinforcement including non-uniform flow, *J. Acoust. Soc. Am.*, **146**(3), 1774–1785.

Tarantola, A., 1984. Inversion of seismic reflection data in the acoustic approximation, *Geophysics*, **49**(8), 1259.

Tarantola, A., 1986. A strategy for nonlinear elastic inversion of seismic reflection data, *Geophysics*, **51**(10), 1893–1903.

Tarantola, A., 1987. *Inverse Problem Theory. Methods for Data Fitting and Model Parameter Estimation*, Elsevier Science.

Tarantola, A., 2005. *Inverse Problem Theory and Methods for Model Parameter Estimation*, Society for Industrial and Applied Mathematics.

Tromp, J., 2019. Seismic wavefield imaging of earth's interior across scales, *Nat. Rev. Earth Environ.*, **1**(1), 40–53.

Tromp, J., Tape, C. & Liu, Q., 2005. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels, *Geophys. J. Int.*, **160**(1), 195–216.

Tromp, J., Komatitsch, D. & Liu, Q., 2008. Spectral-element and adjoint methods in seismology, *Commun. Comput. Phys.*, **3.1–32**

Tromp, J., Tape, C. & Liu, Q., 2014. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels: Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels, *Geophys. J. Int.*, **160**(1), 195–216.

Vera Rodriguez, I., Näsholm, S.P. & Le Pichon, A., 2020. Atmospheric wind and temperature profiles inversion using infrasound: An ensemble model context, *J. Acoust. Soc. Am.*, **148**(5), 2923–2934.

Virieux, J., 1986. P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method, *Geophysics*, **51**(4), 889.

Virieux, J. & Operto, S., 2009. An overview of full-waveform inversion in exploration geophysics, *Geophysics*, **74**(6), WCC1–WCC26.

Vorobeva, E., De Carlo, M., Le Pichon, A., Espy, P.J. & Näsholm, S.P., 2021. Benchmarking microbarom radiation and propagation model against infrasound recordings: a vespagram-based approach, *Ann. Geophys.*, **39**(3), 515–531.

Vorobeva, E., Assink, J., Espy, P.J., Renkwitz, T., Chunchuzov, I. & Näsholm, S.P., 2023. Probing gravity waves in the middle atmosphere using infrasound from explosions, *J. Geophys. Res.: Atmospheres*, **128**(13), e2023JD038725.

Vorobeva, E., Eggen, M.D., Midtfjord, A.D., Benth, F.E., Hupe, P., Brissaud, Q., Orsolini, Y. & Näsholm, S.P., 2024. Estimating stratospheric polar vortex strength using ambient ocean-generated infrasound and stochastics-based machine learning. *Q. J. R. Meteorol. Soc.*, **150**, 2712–2726

Waxler, R., Blom, P., Talmadge, C.L. & Assink, J., 2011. On the waveforms of infrasonic returns from the stratosphere., *J. Acoust. Soc. Am.*, **129**(4,Supplement), 2444.

Xie, Z., Komatitsch, D., Martin, R. & Matzen, R., 2014. Improved forward wave propagation and adjoint-based sensitivity kernel calculations using a numerically stable finite-element PML, *Geophys. J. Int.*, **198**(3), 1714–1747.
APPENDIX A: NOTATIONS

<figure><img src="image_10.png" /><figcaption>**Figure A1.** Illustration of the used numerical scheme: the staggered scheme.</figcaption></figure>

<table><caption>**Table A1.** Notations.</caption><thead><tr><th></th><th>Dimension quantities</th></tr></thead><tbody><tr><td><b>x</b> = (x<sub>1</sub>, x<sub>2</sub>) ∈ ℝ<sup>2</sup></td><td>Spatial coordinate (m) x<sub>1</sub> is the horizontal range and x<sub>2</sub> is the altitude.</td></tr><tr><td><i>t</i> ∈ ℝ</td><td>Time (s)</td></tr><tr><td>Ω</td><td>Spatial domain of the propagation</td></tr><tr><td>Γ</td><td>Boundary of Ω</td></tr><tr><td><i>T</i><sub>max</sub></td><td>Final time of observations</td></tr><tr><td></td><td>Physical quantities</td></tr><tr><td><i>ρ</i>(<b>x</b>, <i>t</i>)</td><td>Density (kg.m<sup>-3</sup>)</td></tr><tr><td><b>v</b>(<b>x</b>, <i>t</i>)</td><td>Velocity (m.s<sup>-1</sup>)</td></tr><tr><td><i>p</i>(<b>x</b>, <i>t</i>)</td><td>Pressure (Pa)</td></tr><tr><td><i>c</i>(<b>x</b>, <i>t</i>)</td><td>Wave speed (m.s<sup>-1</sup>)</td></tr><tr><td><i>T</i>(<b>x</b>, <i>t</i>)</td><td>Temperature (K)</td></tr><tr><td>γ(<b>x</b>)</td><td>Heat capacity ratio</td></tr><tr><td><i>R</i> = 8.3145</td><td>Gas constant (J mol<sup>-1</sup>.K<sup>-1</sup>)</td></tr><tr><td><b>g</b>(<b>x</b>)</td><td>Gravitational acceleration (m.s<sup>-1</sup>)</td></tr><tr><td></td><td>Model</td></tr><tr><td><b>m</b>(<b>x</b>)</td><td>Model with <i>M</i> parameters</td></tr><tr><td><i>ρ</i><sub>0</sub>(<b>x</b>)</td><td>Background density (kg.m<sup>-3</sup>)</td></tr><tr><td><b>v</b><sub>0</sub>(<b>x</b>)</td><td>Wind (m.s<sup>-1</sup>)</td></tr><tr><td><i>p</i><sub>0</sub>(<b>x</b>)</td><td>Background pressure (Pa)</td></tr><tr><td></td><td>Perturbation quantities—State variables</td></tr><tr><td><b>u</b>(<b>x</b>, <i>t</i>)</td><td>State variable</td></tr><tr><td>ρ<sup>*</sup>(<b>x</b>, <i>t</i>)</td><td>Perturbation of density (kg.m<sup>-3</sup>)</td></tr><tr><td><b>v</b><sup>*</sup>(<b>x</b>, <i>t</i>)</td><td>Velocity perturbation (m.s<sup>-1</sup>)</td></tr><tr><td><i>p</i><sup>*</sup>(<b>x</b>, <i>t</i>)</td><td>Perturbation of pressure (Pa)</td></tr><tr><td></td><td>Adjoint variables</td></tr><tr><td>λ(<b>x</b>, <i>T</i><sub>max</sub> − <i>t</i>)</td><td>Adjoint variable</td></tr><tr><td>ρ<sup>†</sup>(<b>x</b>, <i>t</i>)</td><td>First adjoint variable (m.kg<sup>-1</sup>)</td></tr><tr><td><b>v</b><sup>†</sup>(<b>x</b>, <i>t</i>)</td><td>Second adjoint variable (s<sup>2</sup>.kg<sup>-1</sup>)</td></tr><tr><td><i>p</i><sup>†</sup>(<b>x</b>, <i>t</i>)</td><td>Third adjoint variable (s<sup>2</sup>.m<sup>-1</sup>.kg<sup>-1</sup>)</td></tr><tr><td></td><td>Sensitivity kernel</td></tr><tr><td>∂<sub><i>t</i></sub>ρ<sub>0</sub> = K<sub>ρ0</sub></td><td>Sensitivity kernel of infrasound waveform with respect to density (m<sup>3</sup>.kg<sup>-1</sup>)</td></tr><tr><td>∂<sub><i>v</i></sub>ρ<sub>0</sub><sup>*</sup> = K<sub>ρ0</sub><sup>*</sup></td><td>Sensitivity kernel of infrasound waveform with respect to horizontal wind (s.m<sup>-1</sup>)</td></tr><tr><td>∂<sub><i>p</i></sub>ρ<sub>0</sub> = K<sub>p0</sub></td><td>Sensitivity kernel of infrasound waveform with respect to pressure (Pa<sup>-1</sup>)</td></tr></tbody></table>
## B1 Full Navier-Stokes equations with gravity and without dissipation

Density, velocity and energy allow us to describe the fluid under study, which is the atmosphere. Assumptions (H2), (H1) and (H3) allow us to define the Navier-Stokes equations, such as presented in Brissaud et al. (2017)) to model the infrasound propagation.

The last equation of the system of equations is equivalent to an equation of pressure since pressure is related to density, velocity and internal energy (Martire et al. 2022). We favour the equation of pressure since the variation of pressure induced by the infrasound is directly measured.

Moreover, in this approach, we neglect viscous effects (assumption H4) to solve the problem in a simpler case.

Considering a specific event satisfying assumption (H8), we get the following system of equations, also called the Euler equations.

$$ egin{cases} 
\partial_t \rho + \nabla \cdot (\rho \mathbf{v}) = f_{\rho} \\ 
\partial_t (\rho \mathbf{v}) + \nabla p + \nabla (\rho \mathbf{v} \otimes \mathbf{v}) - \rho \mathbf{g} = 0 \\ 
\partial_t p + \mathbf{v} \cdot \nabla p + \gamma p \nabla \cdot \mathbf{v} = f_{\rho} = c^2 f_{\rho} 
\end{cases} 
\quad (B1) $$

with $f_\rho$ the external contribution in pressure, such that $f_\rho = c^2 f_\rho$.

## B2 Numerical implementation of the forward propagation

### B2.1 Numerical grid

To numerically solve the wave propagation problem, we exploit the finite differences method and discretize our domain spatially and temporally. We note $\Delta t$ the time step and $\Delta x_i$ the space step in the $i^{th}$ direction.

In 2-D, each quantity $q$ is written in discretized form as:

$$ q(x_1, x_2, t) \approx q(i \Delta x_1, j \Delta x_2, n \Delta t) = q_{i,j}^n \quad \forall i \in [0, L_1], \forall j \in [0, L_2], \forall n \in [0, N], \quad (B2) $$

where $L_1$ and $L_2$ are the number of space steps in each direction and $N$ is the number of time steps.

Concerning temporal discretization, we use a leap-frog scheme (as performed by Virieux 1986) where density and pressure are known at time $t^n = n \times \Delta t \forall n \in [0, N]$, whereas velocity is known at time $t^{n+1/2} = (n + \frac{1}{2})\Delta t \forall n \in [0, N]$. The notation $[\cdot, \cdot]$ defines a range of integers only.

Concerning 2-D spatial discretization, we use a staggered grid (see Yee 1966 and Virieux 1986) since it efficiently provides a numerically stable high-order discretization (Brissaud et al. 2016). Thus, $\forall (i, j) \in [0, L_1] \times [0, L_2]$, we have on:

(i) density and pressure are defined on $(i \Delta x_1, (j - \frac{1}{2}) \Delta x_2)$;

(ii) the first component of the velocity perturbation and wind is defined on $((i - \frac{1}{2}) \Delta x_1, (j - \frac{1}{2}) \Delta x_2)$;

(iii) the second component of the velocity perturbation $(i \Delta x_1, j \Delta x_2)$.

Fig. A1 illustrates this mesh.

### B2.2 Numerical scheme

Time discretization is carried out based on a 'leap-frog' Euler explicit scheme (see Virieux (1986)). The leap-frog scheme is written as follows

$$ egin{cases} 
\rho^{n+1} = \rho^n + f_1(\rho^n, \mathbf{v}^{n+1/2}, p^n; \mathbf{m}) \\ 
p^{n+1} = p^n + f_2(\rho^n, \mathbf{v}^{n+1/2}, p^n; \mathbf{m}) \\ 
\mathbf{v}^{n+3/2} = \mathbf{v}^{n+1/2} + f_3(\rho^{n+1}, \mathbf{v}^{n+1/2}, p^{n+1}; \mathbf{m}) 
\end{cases} 
\quad \forall n \in [0, N], \quad (B3) $$

with $f_1$, $f_2$ and $f_3$ such that:

$$ \mathcal{L}(\mathbf{u}; \mathbf{m}) = 0 \Leftrightarrow \begin{cases} \partial_t \rho' = f_1(\rho', \mathbf{v}', p'; \mathbf{m}) = f_1(\mathbf{u}; \mathbf{m}) \\ \partial_t p' = f_2(\rho', \mathbf{v}', p'; \mathbf{m}) = f_2(\mathbf{u}; \mathbf{m}) \\ \partial_t \mathbf{v}' = f_3(\rho', \mathbf{v}', p'; \mathbf{m}) = f_3(\mathbf{u}; \mathbf{m}) \end{cases} $$

Evaluating the functions $f_1$, $f_2$ and $f_3$ requires computing the spatial derivatives of $\mathbf{u}$ and $\mathbf{m}$. To this end, we use 4th-order finite-difference schemes.

Let a scalar quantity $q$ be defined at every point $(i, j)$ of the grid. We define two schemes:

(i) The derivative of the scalar $q$ is evaluated at the point $(i - \frac{1}{2}, j)$ for $(i, j) \in [0, L_1] \times [0, L_2]$. In that case, we use the centred scheme, defined as follows:

$$ (\partial_{x_1} q) \left( i - \frac{1}{2}, j \right) = \frac{q_{i-2,j} + 27q_{i-1,j} - 27q_{i,j} - q_{i+1,j}}{24 \Delta x_1}. \quad (B4) $$
(ii) The derivative of the scalar $q$ is computed at a specific point $(i, j)$. We use a non-centred upwind scheme as in Brissaud et al. (2016). To properly treat advection terms of the form $v_0 \cdot \nabla$, we use a non-staggered upwind and downwind scheme depending on the sign of the wind $v_0^1$:

$$ \text{if } v_0^1 \le 0, \quad (\partial_{x_1} q)(i, j) = \frac{1}{6\Delta x_1} (q_{i-2,j} - 6q_{i-1,j} + 3q_{i,j} + 2q_{i+1,j}), \qquad (B5) $$

$$ \text{if } v_0^1 > 0, \quad (\partial_{x_1} q)(i, j) = -\frac{1}{6\Delta x_1} (q_{i+2,j} - 6q_{i+1,j} + 3q_{i,j} + 2q_{i-1,j}). \qquad (B6) $$

Since the vertical wind is zero, the derivatives with respect to $x_2$ of any variable $q$ are computed using the non-staggered downwind (as defined in eq. B5 for the derivative of $q$ with respect to the first component).

In 2-D, $\Delta t$, $\Delta x_1$ and $\Delta x_2$ are chosen to respect the Courant-Friedrichs-Lewy numerical stability condition. Similarly to Martin & Couder-Castaneda (2010), we define the numerical stability condition as follows:

$$ \max_x (c_0 + |v_0^1|) \Delta t \sqrt{\frac{1}{\Delta x_1^2} + \frac{1}{\Delta x_2^2}} \le \text{CFL}_{\text{max}} = 0.3. \qquad (B7) $$

## B2.3 Boundary conditions

### B2.3.1 Ground interface Eq. (3) is a no-slip condition at the ground interface at any time (i.e. the particle velocity is set to zero on the ground interface).

We use the 'Image' boundary conditions defined in Garcia et al. (2017). This condition consists of assuming a mirror domain relative to the ground interface (using ghost points). This domain propagates the symmetrical wave of the wave propagating in our domain. The interference of the two waves at the interface ensures that the velocity is zero at the interface, while null vertical gradients are assumed at the ground interface for physical properties (such as density $\rho_0$, pressure $p_0$ and heat capacity ratio $\gamma$) and perturbations $\rho'$ and $p'$.

Hence, we define:

$$ \begin{aligned} p'(x_1, x_2, t) &= p'(x_1, -x_2, t) & p_0(x_1, x_2, t) &= p_0(x_1, -x_2, t) \\ \rho'(x_1, x_2, t) &= \rho'(x_1, -x_2, t) & \rho_0(x_1, x_2, t) &= \rho_0(x_1, -x_2, t) \\ v'(x_1, x_2, t) &= -v'(x_1, -x_2, t) & v_0^1(x_1, x_2, t) &= -v_0^1(x_1, -x_2, t) \\ g^2(x_1, x_2, t) &= -g^2(x_1, -x_2, t) & \gamma(x_1, x_2, t) &= \gamma(x_1, -x_2, t) \end{aligned} \qquad (B8) $$

### B2.3.2 Atmospheric upper boundaries To simulate an infinite domain, we use PML to absorb the wave in outer buffer areas of 20 grid points size each. We use the ADE-PML as defined by Martin et al. (2010). Each PML thickness is equal to $L_b = 2$ km. Derivatives of density, pressure and velocity perturbations with respect to $x_i$ are divided by the stretching factor $s_{x_i}$, defined such as

$$ s_{x_i} = \kappa_{x_i} + \frac{d_{x_i}}{\alpha_{x_i} + i\omega}, \qquad (B9) $$

where $\kappa_{x_i}$, $d_{x_i}$ and $\alpha_{x_i}$ are computed such as in Martin et al. (2010).

To compute the derivatives of the background quantities (density, pressure and wind) with respect to $x_2$, we adopt the factor $s_{x_2}$ of the PML, so that:

$$ s_{x_2} = \kappa_{x_2} + \frac{d_{x_2}}{\alpha_{x_2}}. \qquad (B10) $$

# APPENDIX C: ADJOINT, SENSITIVITY KERNEL AND INVERSION

## C1 Adjoint state variable and adjoint equations

We derive $\tilde{\chi}$ (see eq. 7) with respect to $u$. In the following, we refer to $E_1$, $E_2$ and $E_3$ the following terms $E_1 = \langle \partial_{\rho'} \tilde{\chi}, \delta \rho' \rangle_D$, $E_2 = \sum_i \langle \partial_{v_i'} \tilde{\chi}, \delta v_i' \rangle_D$ and $E_3 = \langle \partial_{p'} \tilde{\chi}, \delta p' \rangle_D$, respectively.

The optimality conditions require that the derivatives of $\tilde{\chi}$ with respect to $\rho'$, $v'$ and $p'$ are zero. Therefore, $E_1 = E_2 = E_3 = 0$, for all $\delta \rho'$, $\delta v'$ and for all $\delta p'$.
After successively performing integration by parts and the change of variable $t' = T_{max} - t$, we pose $\lambda_1(\mathbf{x}, T_{max} - t') = \rho^\dagger(\mathbf{x}, t')$, $\lambda_2(\mathbf{x}, T_{max} - t') = \mathbf{v}^\dagger(\mathbf{x}, t')$ and $\lambda_3(\mathbf{x}, T_{max} - t') = p^\dagger(\mathbf{x}, t')$ and we obtain:

$$
\begin{align}
E_1 &= \int_{\Omega} \int_{0}^{T_{\text{max}}} \delta\rho' \left[ \partial_{t'} \rho^\dagger - \mathbf{v}_0^1 \partial_{x_1} \rho^\dagger - \mathbf{v}^\dagger \cdot \mathbf{g} \right] dt' d\Omega \nonumber \\
E_2 &= \int_{\Omega} \int_{0}^{T_{\text{max}}} \delta\mathbf{v}' \cdot \left[ \rho_0 \partial_{t'} \mathbf{v}^\dagger - \nabla(\rho^\dagger \rho_0) + \rho^\dagger \nabla \rho_0 - \rho_0 \mathbf{v}_0^1 \partial_{x_1} \mathbf{v}^\dagger \right. \nonumber \\
&\qquad \left. + \rho_0 (\nabla \mathbf{v}_0^1) \mathbf{v}^{\dagger 1} - \nabla(\gamma p_0 p^\dagger) + p^\dagger \nabla p_0 \right] dt' d\Omega \tag{C1} \\
E_3 &= \int_{\Omega} \int_{0}^{T_{\text{max}}} \delta p' \left[ \partial_{t'} p^\dagger - \nabla \cdot \mathbf{v}^\dagger - \mathbf{v}_0^1 \partial_{x_1} p^\dagger + \mathcal{A}^2 \frac{p' - p^{\text{obs}}}{\|p^{\text{obs}}\|^2} \mathbb{1}_{|\mathbf{x}=x_r|} \right] dt' d\Omega \nonumber
\end{align}
$$

with $p'$, $p^{\text{obs}}$ and $\mathcal{A}$ evaluated at $T_{\text{max}} - t'$, and $\rho^\dagger$, $\mathbf{v}^\dagger$ and $p^\dagger$ evaluated at $t'$.

Note that, to cancel the boundary terms, we assumed:

$$ \rho^{\dagger}(\mathbf{x}, 0) = \mathbf{v}^{\dagger}(0) = p^{\dagger}(\mathbf{x}, 0) = 0 \quad \forall \mathbf{x} \in \Omega, \qquad (C2) $$

$$ \mathbf{v}^{\dagger}(\mathbf{x}, t') = (\mathbf{v}^{\dagger}(\mathbf{x}, t')^{1}, 0)^{T}, \quad \forall \mathbf{x} \in \Gamma_{1}. \qquad (C3) $$

In addition, we suppose that adjoint variable is zero on $\Gamma_2 \cup \Gamma_3 \cup \Gamma_4$ since the domain is sufficiently large such that the adjoint field cannot reach the edges during the simulation time.

Consequently, we can identify three equations that define $\rho^\dagger$, $\mathbf{v}^\dagger$ and $p^\dagger$ at any point in space and time. These equations are the adjoint, presented by eq. (9). These equations are associated with the initial conditions, presented by eq. (10), and boundary conditions, expressed in eq. (11).

## C2 Sensitivity kernel

### C2.1 Computations

We derive here the objective function $\tilde{\chi}$ (see eq. 7) with respect to each component of the model $\mathbf{m} = (\rho_0, \mathbf{v}_0^1, p_0)^T$. Let $E_1 = \langle \partial_{\rho_0} \tilde{\chi}, \delta \rho_0 \rangle_M$, $E_2 = \langle \partial_{\mathbf{v}_0^1} \tilde{\chi}, \delta \mathbf{v}_0^1 \rangle_M$ and $E_3 = \langle \partial_{\rho_0} \tilde{\chi}, \delta p_0 \rangle_M$.

After successively performing integration by parts, the change of variable $t' = T_{max} - t$ and substituting $(\lambda_1, \lambda_2, \lambda_3)(\mathbf{x}, T_{max} - t')$ by $(\rho^\dagger, \mathbf{v}^\dagger, p^\dagger)^T(\mathbf{x}, t')$, we have:

$$
\begin{align}
E_1 &= \int_{\Omega} \int_{0}^{T_{\text{max}}} \delta\rho_0 \left[ \rho^\dagger \nabla \cdot \mathbf{v}' - \partial_{x_2} (\rho^\dagger \mathbf{v}'^2) - \mathbf{v}^\dagger \cdot \partial_{t'} \mathbf{v}' + \mathbf{v}^\dagger \cdot (\mathbf{v}_0^1 \partial_{x_1} \mathbf{v}') \right. \nonumber \\
&\qquad \left. + \mathbf{v}^{\dagger 1} (\mathbf{v}'^2 \partial_{x_2} \mathbf{v}_0^1) + p^\dagger \frac{\gamma p_0}{\rho_0^2} f_\rho \right] dt' d\Omega \tag{C4} \\
E_2 &= \int_{\Omega} \int_{0}^{T_{\text{max}}} \delta\mathbf{v}_0^1 \left[ \rho^\dagger \partial_{x_1} \rho' + \mathbf{v}^\dagger \cdot (\rho_0 \partial_{x_1} \mathbf{v}') - \partial_{x_2} (\rho_0 \mathbf{v}'^2 \mathbf{v}^{\dagger 1}) + p^\dagger \partial_{x_1} p' \right] dt' d\Omega \nonumber \\
E_3 &= \int_{\Omega} \int_{0}^{T_{\text{max}}} \delta p_0 \left[ p^\dagger \gamma \nabla \cdot \mathbf{v}' - \partial_{x_2} (p^\dagger \mathbf{v}'^2) - p^\dagger \frac{\gamma}{\rho_0} f_\rho \right] dt' d\Omega \nonumber
\end{align}
$$

with $\rho'$, $\mathbf{v}'$, $p'$ and $f_\rho$ evaluated at $T_{\text{max}} - t'$ and $\rho^\dagger$, $\mathbf{v}^\dagger$ and $p^\dagger$ evaluated at $t'$.

Noting that $E_1 = \int_0^{T_{\text{max}}} \int_{\Omega} \delta\rho_0 K_{\rho_0} \, d\Omega dt'$, $E_2 = \int_0^{T_{\text{max}}} \int_{\Omega} \delta\mathbf{v}_0 K_{\mathbf{v}_0} \, d\Omega dt'$ and $E_3 = \int_0^{T_{\text{max}}} \int_{\Omega} \delta p_0 K_{\rho_0} \, d\Omega dt'$, we have the expression of the sensitivity kernels $K_{\rho_0}$, $K_{\mathbf{v}_0}$ and $K_{\rho_0}$, defined by eq. (14) to (16).

## C3 Optimization algorithm

The inversion algorithm is formulated in algorithm 1.
Require: $\mathbf{m}^0, N_{\max}, \epsilon, m$

LIST_M[k] ← [];

LIST_dF[k] ← [];

for k = 0, ..., 4 do

$F^k \leftarrow \tilde{\chi}(\mathbf{m}^k);$

/* Compute cost function $\chi$ */

$dF^k \leftarrow \partial_m \tilde{\chi}(\mathbf{m}^k);$

/* Compute derivative of $\chi$ */

$dF^k \leftarrow$ Interpolate at source location;

LIST_M[k] ← $\mathbf{m}^k$;

/* Save current model */

LIST_dF[k] ← $dF^k$;

/* Save current gradient */

$r^k \leftarrow -dF^k$;

/* Set the descent direction */

if $k = 1$ then $\alpha_1 \leftarrow 1$ else $\alpha_1 \leftarrow \min(1, \frac{F^k - F^{k-1}}{r^k \cdot dF^k})$ end if

$\alpha \leftarrow \text{backtracking}(\alpha_1);$

$\mathbf{m}^{k+1} \leftarrow \mathbf{m}^k + \alpha r^k;$

end for

$n \leftarrow \text{length}(\mathbf{m}^k);$

$\text{condition} \leftarrow \text{False;}$

while $\text{condition} == \text{False}$ do

$F^k \leftarrow \chi(\mathbf{m}^k);$

/* Compute cost function $\chi$ */

$dF^k \leftarrow \partial_m \chi(\mathbf{m}^k);$

/* Compute derivative of $\chi$ */

$dF^k \leftarrow$ Interpolate at source location;

LIST_M[k] ← $\mathbf{m}^k$;

/* Save current model */

LIST_dF[k] ← $dF^k$;

/* Save current gradient */

Discard model $\mathbf{m}^{k-m}$ in LIST_M and gradient $dF^{k-m}$ in LIST_dF;

$r^k \leftarrow \text{lbfgs}(\text{LIST_M}, \text{LIST_dF}, m);$

/* Descent direction computed with Algorithm 7.4 of Nocedal & Wright

(1999) */

/* Search the step length

*/

$\alpha_1 \leftarrow \min(1, \frac{F^k - F^{k-1}}{r^k \cdot dF^k})$

$\alpha \leftarrow \text{line\_search}(\alpha_1, \alpha_{\max} = 20)$

/* Updating

*/

$\mathbf{m}^{k+1} \leftarrow \mathbf{m}^k + \alpha r^k;$

/* Check the stopping condition

*/

$\text{condition} \leftarrow k > N_{\max}$ or $\|\mathbf{m}^{k+1} - \mathbf{m}^k\| < \epsilon;$

end while

return $\mathbf{m}^{k+1}$;