# Simulation of multiple scattering of seismic waves: energy, displacement and its gradients

Mirko Bracale, Ludovic Margerin, Romain Brossier, Michel Campillo

► To cite this version:

Mirko Bracale, Ludovic Margerin, Romain Brossier, Michel Campillo. Simulation of multiple scattering of seismic waves: energy, displacement and its gradients. Geophysical Journal International, 2025, 244, 10.1093/gji/ggaf504. hal-05527015

HAL Id: hal-05527015

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

Submitted on 25 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
# Simulation of multiple scattering of seismic waves: energy, displacement and its gradients

Mirko Bracale,<fnref n="1" /> Ludovic Margerin,<fnref n="2" /> Romain Brossier<fnref n="1" /> and Michel Campillo<fnref n="1" />

Accepted 2025 November 28. Received 2025 November 24; in original form 2025 May 3

## SUMMARY

The aim of this study is to assess the potential of rotational and strain measurements to provide complementary information on seismic wave scattering, in addition to the conventional seismological observables. We begin by evaluating the accuracy of numerical solutions to the elastic wave equation, solved via the spectral element method, for modelling wave propagation in 3-D complex heterogeneous media. These simulations are benchmarked against predictions from the radiative transfer equation (RTE), which models energy transport in scattering media. The comparison focuses on key scattering parameters: mean free path, diffusion onset and temporal evolution of P/S energy partitioning. Three levels of velocity heterogeneity (10 per cent, 17 per cent and 25 per cent) are tested in both full-space and half-space configurations. The analysis highlights how scattering strength, numerical accuracy and theoretical assumptions, such as those underlying the Born approximation, affect the agreement between the two modelling approaches. This comparison helps define the conditions under which RTE and wave equation-based simulations produce consistent results. Following this assessment, we analyse the energy envelopes of the displacement wavefield and its spatial gradients. The results demonstrate that rotational measurements preserve source-induced polarization longer than other observables. This persistence can provide valuable information for better constraining the source mechanism. Furthermore, analysis of the rotational components can provide complementary constraints on the medium's elastic and scattering properties.

**Key words:** Computational seismology; Theoretical seismology; Wave propagation; Wave scattering and diffraction.

## 1 INTRODUCTION

In recent years, advances in seismic instrumentation have significantly expanded the observational capabilities beyond traditional measurements by enabling the recording of wavefield gradients. Specifically, distributed acoustic sensing (DAS) systems (Parker et al. 2014), also known as distributed dynamic strain sensing, enable dense measurements of normal strain rate along fibre-optic cables, while rotational sensors provide direct observations of the three components of rotational ground motion, corresponding to the curl of the velocity field (Bernauer et al. 2012). The increasing availability of these measurement technologies motivates a systematic investigation into the potential of alternative observables to enhance the characterization of the seismic wavefield. This study focuses on the scattering process, a fundamental aspect of seismic wave propagation in complex heterogeneous media, using numerical simulations of the elastic wave equation. As one of the most

widely used tools in seismology, the solution of the elastic wave equation provides a deterministic description of the seismic wavefield, enabling detailed modelling of wave propagation in heterogeneous environments and allowing for direct comparison with all types of observed seismic data.

To support this analysis, we have first to assess the conditions under which numerical simulations of the elastic wave equation can reliably reproduce theoretical expectations for wave scattering in complex heterogeneous media. This is achieved by comparing simulation results with predictions from radiative transfer theory (RTT) and its mathematical formulation, the radiative transfer equation (RTE) (Weaver 1990; Turner & Weaver 1994a; Ryzhik et al. 1996; Margerin 2005), which provides a statistical description of wave intensity as it propagates through a scattering medium. Applications of the RTE include inverting observed energy envelopes to constrain the medium's attenuation properties (e.g. Hoshiwa 1993; Lacombe et al. 2003; Przybilla et al. 2006; Sens-Schönfelder et al.

<footnote n="1">¹ISTerre, Université Grenoble Alpes, CNRS, Université Savoie Mont Blanc, IRD, Université Gustave Eiffel, 38000 Grenoble, France. E-mail: mirko.bracale@univ-grenoble-alpes.fr</footnote>

<footnote n="2">²Institut de recherche en astrophysique et planétologie (IRAP/OMP), CNRS, CNES, Université Toulouse III - Paul Sabatier, 31400 Toulouse, France</footnote>
2009), characterizing scattering processes and evaluating the time-
dependent sensitivity of coda waves (Obermann et al. 2013, 2016;
van Dinther et al. 2021; Barajas et al. 2022).

Theoretical derivations define the conditions under which the
RTE is expected to approximate statistical properties of wavefields
governed by the elastic wave equation. This typically occurs in the
low-frequency, weak-scattering regime relative to the geometrical
optics limit (Weaver 1990; Ryzhik et al. 1996). Building on this
foundation, several studies over the past decades have compared
RTE and wave equation approaches (e.g. Nakahara & Yoshimoto
2011; Obermann et al. 2016; Emoto & Sato 2018; Celorio et al.
2022). However, the consistency between the two frameworks under
realistic conditions is not yet fully understood, making numerical
validation essential for establishing its limits.

In this study, the elastic wave equation is solved using the spectral
element method (SEM), as implemented in the SEM46 code (Trinh
et al. 2019). This implementation builds upon the foundational
formulations by Komatitsch & Vilotte (1998) and Komatitsch &
Tromp (1999). SEM is widely recognized for its favourable balance
between computational efficiency and accuracy, employing fourth-
order polynomial interpolation within spectral elements. There are
two primary factors that may introduce discrepancies between nu-
merical simulations and actual seismic wavefield behaviour. The
numerical accuracy of SEM has traditionally been assessed in rel-
atively simple configurations, such as homogeneous or stratified
media, typically in terms of the number of grid points per mini-
mum wavelength (Virieux & Operto 2009; Igel 2017; Trinh et al.
2019). However, additional sources of numerical error arise in more
complex settings. Specifically, the cumulative effect of time inte-
gration steps and repeated scattering due to high-contrast hetero-
geneities can significantly impact the fidelity of the simulation. As
the total simulation time and the number of wave-heterogeneity
interactions increase, discretization and interpolation errors may
accumulate, potentially degrading the accuracy of the computed
wavefield in strongly scattering media. Secondly, the RTE is typi-
cally solved using scattering cross-sections derived under the Born
approximation, as described by Weaver (1990) and Ryzhik et al.
(1996). This approximation is considered valid when the square
of the non-dimensional parameter *ka*ε, *k* the wavenumber, *a* the
characteristic size of the scatterers, *ε* the typical perturbation of the
elastic parameters, is much smaller than 1 (Rytov et al. 1989). In 3-
D heterogeneous elastic media the transport approach has not been
comprehensively compared to full wave simulations (see however
Khazaie et al. 2017, for equipartitioning). As a result the condi-
tions of applicability of the theory remain to be clarified. Further
details of the simulation methods can be found in Section 1 of the
Supplementary Material.

We investigate complementary aspects of the propagation in het-
erogeneous media and discuss in detail the relevant physical pa-
rameters used to quantify scattering phenomena. We pay particular
attention to the attenuation of the coherent wave, the onset of the
diffusion regime and the temporal evolution of the *P*-to-*S* and more
exotic rotational energy partitioning ratios. We quantify both theo-
retically and numerically the associated scattering parameters: the
mean free path, the seismic diffusivity and the equipartition time,
respectively. For this aim we considered three different case studies
with increasing strengths of heterogeneities: 10 per cent, 17 per
cent and 25 per cent, in both full-space and half-space media, using
different types of source mechanisms. In the $\epsilon$ = 25 per cent case
study, the 20 s simulations at the maximum analysed frequency of
2.8 Hz span approximately 160 *P*-wave mean free times (mft), al-
lowing us to capture the full development of scattering and energy

redistribution. The strong velocity fluctuations considered in this
study are not only commonly observed in rock samples but have
also been inferred from seismological observations. For a compre-
hensive overview of the range of velocity fluctuations observed in
the Earth, we refer the reader to the review by Sato (2019). High
levels of heterogeneity have also been inferred from modelling seis-
mogram envelopes on other planetary bodies, such as the Moon (e.g.
Dainty et al. 1977; Gillet et al. 2017; Onodera et al. 2021) and Mars
(van Driel et al. 2021; Menina et al. 2023; Onodera et al. 2023).
Incorporating such conditions into our analysis is important not
only for testing the limits of theoretical and numerical models, but
also for ensuring that the results remain relevant and robust across
a wide range of geophysical applications.

The relevant theoretical scattering parameters are computed using established formulations. The scattering mean free path is computed using scattering cross-sections derived following the approach of Sato (1984) which is based on the Born approximation (see also Weaver 1990; Ryzhik et al. 1996). The diffusivity of the medium is obtained from the asymptotic (long lapse-time) solution of the transport equation developed by Weaver (1990). These results have been concisely and pedagogically summarized for seismologists by Turner (1998). The global equipartition time, identified as the coupling time between *P* and *S* waves by Turner & Weaver (1994b) is again computed with the Born approximation. Finally, the ensemble average energy envelopes of *P* and *S* waves are numerically evaluated by solving the full RTE using Monte Carlo simulations (Margerin et al. 2000). To ensure conservation of energy, the exact same scattering cross-sections based on the Born approximation are used in the RTE and the mean free path computations (Weaver 1990; Ryzhik et al. 1996).

A significant challenge of our study lies in the fundamental dif-
ferences between the approaches used in scattering theory and those
employed to describe wave propagation via the wave equation.
While intensity can, in principle, be derived from ground motion
obtained via wave equation simulations, it cannot be directly com-
pared with predictions from scattering theory, as it constitutes a
local measure. Consequently, special methodologies are required to
compare the outputs from these two approaches called ‘ensemble
average’ and ‘spatial average’. The first method involves running
multiple simulations of the wave equation with the same scatter-
ing characteristics of the medium but different configurations of
heterogeneities (referred to as realizations). The average of the re-
sults across these different realizations is then used to extract the
statistical characteristics of the wavefield. While this approach is
more accurate, as it accounts for different media that share the same
statistical properties as described by theory, it is computationally
expensive, particularly for 3-D simulations, as extensively discussed
in this work.

In this study, we employed the second method, the spatial av-
erage, which requires only a single realization. This approach is
valid provided that the conditions of complete isotropy and uniform
statistical properties are met. Receivers must be symmetrically dis-
tributed around the source, and averages are calculated at specific
distances. This approach assumes that wave paths are independent,
meaning that waves arriving at different receivers primarily interact
with distinct heterogeneities. This condition is not fully satisfied
but can be considered valid if the receivers are sufficiently spaced.
In this case, the spatial average is assumed to be equivalent to the
ensemble average, following the ergodic assumption (e.g. Shapiro
& Kneib 1993). The only issue arises because heterogeneities lo-
cated near the source may generate a correlation between the signals
observed at all positions (Shapiro 1999). In this work, we focus on
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) Source-receiver configuration. A total of 2521 receivers were uniformly distributed over a spherical surface with an average angular spacing of approximately 5°, based on discretization of the polar angles $θ$ and $φ$ in spherical coordinates. (b) Seismograms of displacement along the radial direction at a distance of 2 km from the source, for increasing strength of heterogeneity, single realization.</figcaption></figure>

the transport regime in which the scattering mean free path $l$ is much greater than the wavelength $λ$. In this case, the contribution of long-range correlations of the order of $λ/l$ can be considered negligible.

We would like to draw the attention of the reader to the fact that while RTE is well established in the full space case, the half-space and waveguide geometries have been the topic of intense research in the last 20 yr (Tregoures & van Tiggelen 2002; Zeng 2006; Maeda et al. 2008; Margerin et al. 2019; de Hoop et al. 2022; Xu et al. 2022). From a practical perspective the key difference between the different approaches lies in the formulation of the coupling between surface and body waves. Tregoures & van Tiggelen (2002) and de Hoop et al. (2022) employ an eigenmode decomposition of the wavefield while Maeda et al. (2008), Margerin et al. (2019) and Xu et al. (2022) use a more straightforward ray theoretical representation of the body waves. In the latter formulations, our numerical approach offers a unique opportunity to estimate the coupling terms between body and surface waves based on the wave equation.

This article is organized into two main parts. In the first part, we perform a reliability assessment, comparing results from SEM simulations with those derived from scattering theory in both full-space and half-space media. It should be noted that we have not considered the possibility of a rough topography, which would also be locally the source of significant conversions between propagation modes (e.g. Sánchez-Sesma & Campillo 1991; Rodgers et al. 2010; Kumagai et al. 2011; Imperatori & Mai 2015; Takemura et al. 2015).

In the second part, we revisit these cases to investigate the energy evolution of the three components of the displacement wavefield and their spatial gradients, namely normal strain and rotation. These quantities are particularly interesting because they are now directly accessible through rotational sensors and DAS, and because such measurements have been shown to be highly sensitive to small heterogeneities (Capdeville & Sladen 2024).

## 2 SIMULATIONS IN A FULL SPACE MEDIUM

In this section, we assess the reliability of the numerical simulations in accurately reproducing the dynamics of seismic wave scattering in a full-space medium. We consider three levels of heterogeneity, $ϵ = 10$ per cent, $ϵ = 17$ per cent and $ϵ = 25$ per cent, although 10

per cent and 25 per cent are the cases we discuss in the most detail in this study. These values define the range within which the validity and robustness of our modelling approach were tested. In all cases, the medium contains heterogeneities characterized by a von Kármán correlation with a correlation length of 320 m. The full simulation parameters are listed in the Supplementary Material, Table S1. For the given parameters, we observe at central frequency of 2 Hz a $P$ and $S$-mft of respectively 1.30 and 0.79 s in the 10 per cent case and 0.25 and 0.15s in the 25 per cent case. In Fig. 1(a), we show the source-receiver configuration used for three sets of receivers on concentric spheres. 2521 receivers were placed on 36 concentric spheres, with a minimum radius of 2 km, increasing by 0.5 km. The source, characterized by an isotropic moment tensor, was positioned at the centre. On these spheres, the receivers' axes are rotated to obtain a radial component directed towards the source and transverse components tangent to the spheres. By averaging observations from these symmetrically distributed observation points, we provide an estimation of the average wave behaviour, minimizing the influence of random fluctuations.

Fig. 1(b) shows the seismogram of radial displacements recorded at 2 km distance from the source in three cases of increasing heterogeneity standard deviation, $ϵ$. It is possible to observe how, with increasing $ϵ$, scattering attenuates ballistic arrivals and transfers the coherent energy to the coda.

### 2.1 Mean free path computation

The scattering mean free path represents the average distance a wave travels between successive interactions with heterogeneities. It is a fundamental parameter in the description of wave propagation in complex media, as it governs the rate at which the wavefield transitions between different scattering regimes. In this section, we estimate the scattering mean free path by analysing the decay of the frequency components of the mean $P$ wave field. These numerical estimates are then compared with theoretical predictions obtained following the approach of Sato (1984). In Figs 2(a) and 3(a), the coherent field is shown at distances ranging from 2 to 14 km from the source in the two cases with $ϵ = 10$ per cent and $ϵ = 25$ per cent. It is evident that as one moves away from the source, the interaction with heterogeneities increases, and in some cases, even spatial averaging is unable to remove the fluctuations. In the first case study, the coherent field remains well preserved even at a distance of 14 km
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Simulation with $\epsilon = 10$ per cent. (a) Coherent field, obtained by averaging the radial component of the displacement at increasing distance from the centre. (b) Comparison between the mean free path computed with SEM and RTE. A good agreement is found between about 1.7 and 3.0 Hz.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Simulation with $\epsilon = 25$ per cent. (a) Coherent field. (b) Comparison between the mean free path computed with SEM and RTE. The frequency dependence is well captured between about 0.5 and 2.5 Hz.</figcaption></figure>

from the source. Conversely, in the second case study, because of the strong scattering even at 6 km distance, the coherent part of the signal can be extracted, but the fluctuations are not adequately suppressed. At greater distances, these fluctuations significantly affect the signal and its spectrum. To quantitatively analyse this phenomenon, we examined the mean field at a reference distance and computed the ratio between the standard deviation of the coda and that of the ballistic arrivals. If the spatial averaging is effective in suppressing fluctuations, the standard deviation of the coda should be close to zero, resulting in a ratio approaching zero. Conversely, if the averaging fails to eliminate these fluctuations, the ratio will be significantly higher. By normalizing with respect to the standard deviation of the ballistic arrivals, we obtain a dimensionless metric that enables consistent comparison across different source-receiver distances and varying wavefield amplitudes. In the simulation with $\epsilon = 10$ per cent, at distances of 4 and 8 km, the ratios obtained were 0.010 and 0.012, respectively. In the second simulation, where $\epsilon = 25$ per cent, the ratios obtained were 0.176 and 0.424 at the same distances. This indicates that, despite spatial averaging, a substantial amount of energy remains in the averaged coda in the second case. As a result, comparison with theoretical predictions becomes difficult due to the limited effectiveness of the averaging procedure. The mean free paths obtained from the numerical simulations were then compared with the theoretical results given by Sato (1984), Sato et al. (2012) and the computation procedure is described in the Supplementary Material Section 2.

The results, shown in Fig. 2(b), where $\epsilon = 10$ per cent, indicate that our findings align with theoretical expectations, at least within the frequency range of 1.7 to 3 Hz. These boundaries are set by the simulation configuration. The lower limit is due to the fact that when the scattering mean free path becomes comparable to the medium size, wave attenuation is too weak to be measured effectively. The upper limit arises from the issue that if the scattering mean free path is too small relative to receiver spacing, signal decay occurs too rapidly to be measured.

This analysis allows us to assert that frequencies within this range are reliable. However, it is important to note that frequencies outside of this range may also be well represented. The accuracy of this measurement could be influenced solely by the method used. In the following analyses, good agreement between SEM and RTE has been found in the range of frequency between 0.6 and 2.8 Hz.

In the case $\epsilon = 25$ per cent, Fig. 3 shows that the computed mean free path agrees with expectations only within a limited frequency range, specifically between about 0.6 and 2.5 Hz. The difference from the $\epsilon = 10$ per cent case arises because, at low frequencies, when the mean free path becomes comparable to the size of the medium, it is no longer possible to correctly estimate the decay. At high frequencies, the discrepancy is due to the fact that the averaging process does not fully eliminate the fluctuations, since much stronger scattering is produced (see Fig. 3a). Furthermore, the limits of the Born approximation are reached: for example, above 3.5 Hz, the mean free path becomes smaller than
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Comparison between energies derived from SEM simulations and predictions of the diffusion approximation. The displacement energy $E_u = u_x^2 + u_y^2 + u_z^2$, rotation energy $E_\omega = \omega_x^2 + \omega_y^2 + \omega_z^2$ and strain energy $E_\epsilon = \epsilon_{xx}^2 + \epsilon_{yy}^2 + \epsilon_{zz}^2$ are shown for two levels of heterogeneity. Computations correspond to source to receiver distance of 6 km. The diffusion regime starts at about 12 s in the $\epsilon = 10$ per cent simulation (a) and at about 7 seconds in the $\epsilon = 25$ per cent simulation (b).</figcaption></figure>

the correlation length. This issue is discussed at the end of Section 2.3.

## 2.2 Transition from multiple scattering to diffusion

As lapse time increases, seismic waves interact with an increasing number of heterogeneities, gradually transitioning from the single and multiple-scattering regimes to the diffusive regime. The rate at which this transition occurs is governed by the scattering mean free time, defined as the ratio between the scattering mean free path and the average phase velocity. After multiple scattering events, the wavefield reaches the diffusive regime, which is characterized by a weakly anisotropic energy flux and an energy density that evolves according to the diffusion equation. The solution of this equation for a 3-D medium is the following:

$$E(\mathbf{r}, t) = \frac{S_0}{(4\pi Dt)^{3/2}} e^{-\frac{t^2}{4Dt}}, \quad (1)$$

where $S_0$ is the source term and $D$ is the diffusivity. Close to the source and for large $t$, the decay of the energy with time is proportional to the factor $t^{-3/2}$. In this paragraph, we examine whether the observed decay in the simulations is consistent with this theoretical expectation. A fundamental parameter of this equation is the characteristic time or Thouless time, defined as $\tau_D = r^2/D$, which quantifies the rise time of the energy envelope. For the computation of this parameter the diffusivity was derived by fitting the observed energy decay to the diffusion equation solution. The comparison between the diffusivity constant, in the 25 per cent case, obtained in the SEM simulations (3.4 km$^2$ s$^{-1}$) and that obtained from the RTE (3.2 km$^2$ s$^{-1}$) was found to be consistent. This is illustrated by the fit shown in the Fig. S2 of the Supplementary Material.

Fig. 4 shows the energy decay observed at a distance of 6 km for the two cases: (a) $\epsilon = 10$ per cent, corresponding to a Thouless time $\tau_D = 2.2$ s, and (b) $\epsilon = 25$ per cent, corresponding to $\tau_D = 10.6$ s.

The energy has been averaged considering the receivers located at that distance. In the first case, Fig. 4(a), a good fit between energy decay and diffusion equation solution becomes visible starting at around 12 s, 8 P-mfts, from P waves arrival, while in the second case, Fig. 4(b), a good fit is achieved already at 7 s, 13 P-mfts. This difference arises from variations in scattering anisotropy, mean free time and consequently diffusivity, all of which depend on the heterogeneity strength $\epsilon$, which differs significantly between the two

scenarios. As $\epsilon$ increases, the mean free time decreases, leading to a faster transition between scattering regimes. The theoretical values of the mean free times are shown in Fig. S3 of the Supplementary Material. Note that the time of transition between regimes is not dependent of the source mechanism (Morioka et al. 2017).

## 2.3 Energy partition ratio

The interaction of seismic waves with heterogeneities leads to energy conversion between longitudinal ($E_P$) and shear ($E_S$) components. The rate of this conversion is also governed by the scattering mean free times. The system ultimately converges toward an asymptotic constant energy ratio, which depends on the mean elastic properties of the medium but not on their fluctuations. This state, known as equipartition, is a fundamental concept in multiple scattering theory, stipulating that at sufficiently long lapse times, the wave energy becomes evenly distributed among all propagation modes supported by the medium. We refer the reader to Weaver (1990), Ryzhik et al. (1996) and de Hoop et al. (2022) for the theoretical framework and to Khazaie et al. (2017) for a detailed numerical investigation in 3-D using SEM simulations. In this paragraph, we compare the energy partition ratios obtained from SEM simulations with those predicted by radiative transfer theory, in order to assess whether the scattering behaviour in both frameworks is consistent over long lapse times. In SEM simulations, longitudinal and shear energy were computed using the normal strains ($\epsilon_{ij}$) and rotations ($\omega_i$), which provide access to the respective contributions of P and S-wave energy in the wavefield.

This calculation is simplified by the fact that the S field has a null divergence while the P field is irrotational. These properties allow us to define the energies in the following way:

$$E_P = (\lambda + \mu)(\nabla \cdot \mathbf{u}(t))^2 = (\lambda + \mu)(\epsilon_{xx} + \epsilon_{yy} + \epsilon_{zz})^2 \quad (2)$$

$$E_S = \mu(\nabla \times \mathbf{u}(t))^2 = \mu(2\omega_x^2 + 2\omega_y^2 + 2\omega_z^2), \quad (3)$$

where $(\mathbf{x}, \mathbf{y}, \mathbf{z})$ are a set of orthogonal reference axes. Once $E_P$ and $E_S$ have been obtained, it is then possible to analyse the energy partition as the ratio of the two. Because of the strong fluctuations, the energy was smoothed, before computing the ratio, using a window of 4 s length.

In this section, we show the evolution of the energy partition ratio observed at a distance of 6 km from the source, in two cases
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Energy partition comparison, SEM computation and RTE solved with the Born approximation for a source-to-receiver distance of 6 km. (a) $\epsilon = 10$ per cent (b) $\epsilon = 25$ per cent.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Source-receiver configuration. 1500 receivers have been deployed on each circle, centred at the source location.</figcaption></figure>

(Fig. 5). To ensure comparable frequency content between the two simulation methods, the SEM signal was filtered between 0.6 and 2.2 Hz.

In both case studies, Figs 5(a) and (b), we solved the RTE and compared the results with the SEM computations. The main difference between the two simulations is that while in the SEM simulations the source function is a Ricker wavelet, in the Monte Carlo simulations the source is a single pulse. To align with this difference, we used the method of Celorio *et al.* (2022) and conducted an additional simulation in a completely homogeneous medium. The solution of the RTE was subsequently convolved with the observed source function at a 6 km distance to simulate the same source conditions as in the SEM simulations. This allowed us to analyse the entire evolution of the partition ratio, starting from the ballistic arrivals. In both cases we observe an excellent match between the results obtained from SEM simulations and the energy ratio predicted by RTE. In the first case, where $\epsilon = 10$ per cent, we observe the dynamic of energy partition up to a certain time lapse, which is insufficient to achieve complete equilibrium, while in the second case, where $\epsilon = 25$ per cent, we can observe the transition up to the asymptotic limit, $0.5 \cdot v_s^3/v_p^3 \sim 0.1$, which occurs at about 30 P-mfts.

A note here concerns the applicability of the RTE coupled with the Born approximation in the 25 per cent case. Theoretically the Born approximation is expected to be valid in the regime $(ka\epsilon)^2 \ll 1$ with $k$ the wavenumber and $a$ the correlation length of fluctuations (Rytov *et al.* 1989). In the case $\epsilon = 25$ per cent, we certainly reach the theoretical limit of applicability $(k^2a^2\epsilon^2 \approx 0.6)$. In practice, the application threshold appears to be flexible. In this study, we demonstrate that the Born approximation appears to be valid up to 25 per cent velocity perturbations as far as energy partition ratios are concerned.

We suspect that the range of validity of the Born approximation might be broader for certain observables than for others. Indeed, it is well known that the Born approximation tends to overpredict the scattering amplitude around the forward direction, which eventually leads to its breakdown in the high-frequency (geometrical optics) limit. This breakdown is apparent in Fig. 3(b) where the mean free path calculated theoretically above approximately 3.5 Hz becomes smaller than the correlation length, which is obviously non-physical. Nevertheless, it is worth pointing out that in the computation of certain transport quantities such as the diffusivity or the equipartition time, a weighting factor $(1 - \cos\Theta)$, with $\Theta$ the scattering angle, appears in the integrand (Turner 1998). This factor de facto eliminates the offending directions. For these observables specifically, we may expect that transport theory (with Born approximation) is valid for perturbations larger than 10 per cent as suggested by our study. Finally, we note that in the case of scalar and electromagnetic waves, modified RTEs have been derived that allow the treatment of resonant discrete scatterers, thereby extending the validity of transport theory beyond the Born approximation (Barabanenkov *et al.* 1995; Barabanenkov & Barabanenkov 1997).

## 3 SIMULATION IN A HALF-SPACE MEDIUM

In this section, we repeat the computation of the mean free path, together with the evaluation of the transition from the multiple scattering to the diffusion regime and of the equipartition ratio, this time assuming an isotropic source and accounting for the presence of the free surface. The media analysed here are characterized by the same scattering properties as in the full-space case. The source-receivers configuration shown in Fig. 6 was used. 1500 receivers were arranged along concentric circles starting from a radius of 2 km, increasing by 1 km. The source was positioned at the centre of these circles but at a depth of 300 m.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Comparison of two simulations in a heterogeneous half-space medium with $\epsilon = 10$ per cent (a) and $\epsilon = 25$ per cent (b). Radial field observed increasing source-to-receiver distance. In the first case (a), *P* and Rayleigh waves are clearly visible, while in the second case (b), the average is not efficient enough to remove the fluctuations, and therefore the presence of coherent body and surface waves does not appear clearly. The signals have been filtered between 0.9 and 2.8 Hz.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Mean free paths of the half-space, for *P* (*P*<sub>HS</sub>) and Rayleigh waves (RW). Comparison with *P* mean free path observed in the full space using SEM (*P*<sub>FS</sub>) and RTE (*P*<sub>FS</sub>(RTE)).</figcaption></figure>

The parameters used in these simulations are once again those listed in the Supplementary Material, Table S1.

## 3.1 Mean free path

The mean field was estimated as in the full space medium, by averaging the radial component at a defined distance from the source. Fig. 7 presents the mean field obtained from the two simulations with $\epsilon = 10$ per cent (a) and $\epsilon = 25$ per cent (b).

In the first case, the main components of the wavefield are clearly visible; ballistic arrivals, that is, *P* waves, and, from 4 km onward, Rayleigh waves can be distinctly observed. In the second case, however, strong fluctuations obscure these components, resulting in a less distinct average compared to the full-space simulation, as shown in Fig. 3(a). This occurs because, in the full-space medium, the receivers were arranged on a spherical surface and more widely spaced. In contrast, the reduced distance between receivers in this case results in nearby receivers capturing very similar signals, making it difficult for the averaging process to fully eliminate the fluctuations.

The mean free path was calculated for both *P* waves and Rayleigh waves. The results are presented in Fig. 8 with a comparison to previous results from the full-space medium for *P* waves. The comparison between *P* mean free paths in full space (*P*<sub>FS</sub>) and half-space (*P*<sub>HS</sub>) highlights the significant influence of the free surface: *P*<sub>FS</sub> is approximately 2 times larger than *P*<sub>HS</sub>, up to 2.4 Hz, suggesting substantial modal coupling, likely toward *S* or surface waves, introduced by the surface. Furthermore, the intermediate value of the Rayleigh wave mean free path (RW), falling between the two *P*-wave estimates, indicates that Rayleigh waves actively participate in the scattering process, though they experience weaker scattering compared to body waves in the half-space.

We also would like to remark that Margerin *et al.* (2019) and Xu *et al.* (2022) have suggested that in the vicinity of the free surface the mean free path of body waves may actually have to be considered as a depth-dependent quantity (even if the variance of the medium fluctuations is actually independent of the depth). This must be so because as the body waves get scattered closer to the free surface, the coupling with the Rayleigh wave becomes increasingly probable. If we neglect any perturbation of body wave to body wave coupling, the additional scattering channel offered by the Rayleigh wave causes extra-attenuation of the body wave compared to the full-space case. This physical scenario is supported by the numerical results shown in Fig. 8.

## 3.2 Transition from multiple scattering to diffusion

The decay of energy, calculated as described in paragraph 2.2, shows that while the dynamics of the observables differ from the full space medium case due to the influence of surface waves, the transition to the diffusive regime occurs later. Fig. 9 shows the analysis of energy dynamics at a source to receiver distance of 2 km. In both case studies, the decay follows a $t^{-3/2}$ trend, characteristic of 3-D diffusion and consistent with the predictions of Margerin *et al.* (2019).

## 3.3 Energy partition ratio

In this section, we compare the energy partition ratio $E_s/E_p$ observed in the SEM simulations of the two case study with the one predicted asymptotically at long lapse-time by equipartition theory (Hennino *et al.* 2001).
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Decay rate of the late coda in a half-space medium. Comparisons of the displacement energy $E_u = u_x^2 + u_y^2 + u_z^2$, rotation energy $E_\omega = \omega_x^2 + \omega_y^2 + \omega_z^2$ and strain energy $E_\epsilon = \epsilon_{xx}^2 + \epsilon_{yy}^2 + \epsilon_{zz}^2$ in the two cases. Envelopes have been computed for a source-to-receiver distance of 2 km. The transition to the diffusive regime starts at about 9 s in the $\epsilon = 10$ per cent simulation (a) and at about 6 s in the $\epsilon = 25$ per cent simulation (b). We show the decay rate, proportional to $t^{-3/2}$, which is consistent with predictions of Margerin *et al.* (2019).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Energy partition ratio of the two case with $\epsilon = 10$ per cent (a) and $\epsilon = 25$ per cent (b) at the free surface for a source-to-receiver distance of 2 km. The asymptotic limit was obtained by Hennino *et al.* (2001).</figcaption></figure>

In the first case, Fig. 10(a), where $\epsilon = 10$ per cent, we observe the arrival of *P* waves at about 2.5 s, the arrival of the *S* and Rayleigh wave train at about 5 s and then we see the ratio rise to about 0.2 before starting to decrease. In this case, the asymptotic limit is not reached, but the temporal evolution suggests a slow convergence toward the expected limit.

By contrast, in Fig. 10(b), after the arrival of the *S* and Rayleigh waves, the ratio stabilizes around the expected value, approximately 0.14 for the medium analysed (Hennino *et al.* 2001). For further examples and detailed discussions, the reader may consult Khazaie *et al.* (2017).

# 4 ENERGY ENVELOPE ANALYSIS

In this section, we present an analysis of the evolution of energy of motion in the three directions for the three observables: displacement $u(t)$, normal strain $\epsilon(t)$ and rotation $\omega(t)$. The aim is to investigate whether they may provide different information about the wavefield. It is natural to study the direction of displacement and, therefore, to use the term ‘polarization’ to describe the distribution of amplitude or energy between the three components of displacement. Similarly to displacement, we can consider the

three components of the rotation vector and analyse the distribution of signal energy using the term ‘polarization’. In the following, we extract the three normal components of the deformation tensor (by analogy with DAS measurements) and again use the term ‘polarization’ to describe the distribution of gradient intensities in the three corresponding specific directions. We examine two case studies: first, the full-space medium, and secondly, the half-space medium. In contrast with previous numerical work focusing on the capability of RTE to correctly model a wide variety of scattering media (e.g. Przybilla *et al.* 2006; Celorio *et al.* 2022), our interest is in the potential of wavefield gradients as new observables to characterize the scattering mechanisms in heterogeneous media. For an application to seismic data, the reader is referred to Gaebler *et al.* (2015) who employ Monte-Carlo solutions of RTE to model rotational motions in the teleseismic *P*-coda detected by a ringlaser.

## 4.1 Full space medium

The energy of individual components of the observables was averaged at specific distances from the source. In this section we analyse separately the two cases with increasing strength of heterogeneity
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Envelope simulation based on SEM with $\epsilon$ = 10 per cent in full space medium. Comparison of the energy evolution of the radial and transverse components of the displacement (a), normal strain (b) and rotation (c). The vertical black line indicates the theoretical arrival of the *P* wave, while the red line marks the arrival of the *S* wave.</figcaption></figure>

$\epsilon$ = 10 per cent and 25 per cent, the case with $\epsilon$ = 17 per cent is shown in the Supplementary Material, Fig. S5.

In the case of least scattering, as shown in Fig. 11(a) for the displacement, a strong ballistic arrival is visible in the radial component, followed by a train of *S* waves visible in the other components. Even though the source is an isotropic moment tensor, heterogeneous near the source can generate conversion from *P* wave to *S* wave. If these waves are generated close to the source, they appear to propagate as if they were emitted directly from the source itself. The observation of *S* waves arrival in a heterogeneous medium, given an isotropic *P*-wave source, is consistent with the works of Hirakawa *et al.* (2016); Margerin (2017) and appears in all the cases analysed. To confirm this hypothesis, it is possible to note that this contribution is absent in the divergence of the field, Fig. 11(b). The normal strain components show features similar to the displacements.

As the scattering strength increases ($\epsilon$ = 25 per cent), the transition to multiple scattering occurs more rapidly. In Figs 12(a) and (b) for $u(t)$ and $\epsilon(t)$, we observe that the ballistic arrivals are less pronounced and the contribution of *S* waves is spread over longer time windows. At a distance of 8 km, for instance, a perfect example of a diffusive regime is shown. In this case, the direct waves are absent and the envelope begins with a smooth increase of the energy, characteristic of the diffusive regime.

Even before the onset of the diffusive regime, repeated scattering interactions in the heterogeneous medium progressively redistribute energy among components. This reduces the initial differences imposed by the source. In complex media with numerous heterogeneous, this redistribution eventually leads the three components

toward comparable energy levels. The analysis of the different observables, however, reveals substantial differences when comparing displacement-strain with rotations.

In the first two observables, displacement and strain, shown in Figs 11(a) and (b), the convergence to the same value of the energy of the three components occurs a few seconds after the arrival of *S* waves. For example, for the displacements at a distance of 6 km (Fig. 11a), it takes place approximately 2 s after the *S*-wave arrival, corresponding to about 3 *S*-mft at a central frequency of 2 Hz. The same behaviour is observed for strain in Fig. 11(b). However, a notable difference appears in the case of rotations (Fig. 11c), where, at the same distance, the convergence occurs around 8 s after the *S*-wave arrival, approximately 10 *S*-mfts.

The important aspect of rotational measurements is that they are exclusively sensitive to *S* waves. High amplitude *P* waves scattered in the vicinity of the isotropic source are converted into *S* waves resulting in a clear arrival of energy after the *S* traveltime in the transverse components. Again, we remark that this observation aligns with previous findings by Pitarka *et al.* (2015) and Hirakawa *et al.* (2016), who demonstrated that such *S*-wave phases are closely tied to velocity perturbations in the vicinity of the source. At finite lapse-times, the flux of energy is strongly directive outwards, away from the source (Paul *et al.* 2005). This results in a predominance of polarization of *S* waves in the plane perpendicular to the flux, that leads for the rotations to a weaker amplitude for the radial component. This effect is clearly visible with rotations, whereas it is impossible to discern in displacements due to the presence of *P* waves that are mainly polarized in the radial direction. This leads
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Envelope simulations based on SEM with $\epsilon = 25$ per cent in full space medium. Comparison of the energy evolution of the radial and transverse components of the displacement (a), normal strain (b) and rotation (c). The vertical black line indicates the theoretical arrival of the *P* wave, while the red line marks the arrival of the *S* wave.</figcaption></figure>

to a weaker radial component and a significantly delayed depolarization in the rotational signal. However, as scattering increases, the depolarization time for rotational components decreases, as shown in Fig. 12(c). In conclusion, rotational measurements can clearly reveal the *S*-wave behaviour and the changes in polarization with increasing scattering, which are difficult to identify from displacement records.

## 4.2 Half-space medium

The presence of the free surface significantly alters the dynamics of the wavefield. The existence of Rayleigh waves propagating at the free surface affects the partition of energy between *P* and *S* in the diffusion limit (Hennino *et al.* 2001). Although the evolution of the partition of energy between body waves and surface waves is difficult to analyse, in particular because of the finite depth of penetration of Rayleigh waves (see e.g. Margerin *et al.* 2019), we expect that the distribution of the wave-field energy over the different components of the observables will be very different from the case of full space. A first-order argument is that the Rayleigh wave displacements have an elliptical polarization in the plane defined by the vertical and the direction of propagation, so there is no contribution for rotations around a vertical axis. Furthermore, observations at the free surface imply relationships between the components of the deformations in order to satisfy the local stress boundary condition. These relationships between the components of the deformations are reflected in the amplitude distributions for the deformations in the different directions, but are independent of the propagation processes in

the medium responsible for the movements at the surface. We can therefore expect to observe overlapping effects of different natures when comparing the wave field in a full space and on a free surface.

In this configuration, the directions of observation have been rotated towards the source's epicentre, and not the source location, as usually done with actual seismic records. This means that considering the widely used cylindrical RTZ coordinate system, *z* corresponds to *Z*, $\Phi$ corresponds to *T* and *r* to *R*. We consider the case of an isotropic source.

In Fig. 13, we show the evolution of the three radial and transverse components of displacement (a), strain (b) and rotation (c) in the simulation with $\epsilon = 10$ per cent.

For the displacements, we observe *P* waves polarized along the radial direction, followed by the arrivals of *S* waves and Rayleigh waves, the latter being clearly identifiable. At distances greater than 5 km from the source and after the arrival time of the direct Rayleigh wave, the envelopes of the *r* and $\phi$ components show lower amplitudes compared to the *z* component. This effect arises because the predominant scattered Rayleigh waves contribute to the horizontal components depending on their propagation direction, whereas Rayleigh waves from all directions contribute to the vertical component.

For the normal strain in Fig. 13(b), two peaks related to body waves and Rayleigh waves are still observed. However, an offset can be observed between the strain energy for the vertical component and for the other components. This difference is due to the free surface conditions. Specifically, at free surface, the three components of the stress, along the three axes $i = x, y, z$, the stress
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Envelope simulations based on SEM with $\epsilon = 10$ per cent in a half-space medium and isotropic source. Comparison of the energy evolution of the radial and transverse components of the displacement (a), normal strain (b) and rotation (c). The vertical black line marks the theoretical P-wave arrival, the red line indicates the S-wave arrival and the blue line corresponds to the Rayleigh waves.</figcaption></figure>

component $\tau_{iz}$ must be zero. For the vertical component we can write (Aki & Richards 2002):

$$ \tau_{zz} = \lambda \left( \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} + \frac{\partial u_z}{\partial z} \right) + 2\mu \frac{\partial u_z}{\partial z} = \lambda \nabla \cdot u + 2\mu \frac{\partial u_z}{\partial z} = 0 \quad (4) $$

This means that:

$$ \nabla \cdot u = - \frac{2\mu}{\lambda} \frac{\partial u_z}{\partial z}. \quad (5) $$

For a Poisson medium, $2\mu/\lambda = 2$. Which means that the amplitude of the vertical strain component is approximately half that of the signal's divergence, and four times less when considering the energy. In the simulations, we find full agreement with the theory. Apart from this effect, the energy dynamics of the strain are similar to those observed in the displacement.

In Fig. 13(c) for rotations, we observe a large peak on the $\Phi$ component, due to the arrival of Rayleigh waves, as expected, and a global dynamics similar to that observed in the full space medium case of Figs 11 and 12. The key difference observed here is that the radial component exhibits significantly less energy compared to the other two components, up until the arrival of the Rayleigh waves. Once the Rayleigh waves arrive, within a few seconds, the horizontal components, that is, radial and $\phi$ components, converge and exhibit similar energy levels, dominating the wavefield. This effect is due to the specific polarization of the Rayleigh waves that, as already mentioned, does not contribute to rotation around the vertical axis.

In the case of $\epsilon = 25$ per cent, we observe a dynamics similar to that shown in the previous case across all observables, with the difference that, due to the lower conversion times, the convergence of the energy of the three components is accelerated. For the displacement, Fig. 14(a), we observe a very quick convergence of the components, a small increase in energy is notable in the vertical component $z$, attributable to the scattered Rayleigh waves. This pattern does not appear in the strain field (b), for which the main characteristics is the lower energy of the vertical component as imposed by the free surface condition.

As in this case with fluctuations of 25 per cent the convergence towards the diffusive regime seems well established, we can compare our results with the values derived from equipartition. Using the computational approach of Margerin (2009), we have calculated the distribution of energy of the rotational wavefield on the vertical and horizontal axes at equipartition. The theory predicts that the ratio $(\nabla \times \mathbf{u})_h^2/(\nabla \times \mathbf{u})_z^2$, (the subscript $h$ designates the sum of the two horizontal components; see Supplementary Material, Section 3, for the full formula), should tend to 5.3 at the surface of a Poisson solid. This result is in good agreement with the numeric shown in the Fig. S4 of Supplementary Material. The analysis of the case with a source corresponding to a dipole $M_{zz} \neq 0$, presented in the Supplementary Material (Fig. S6), shows dynamics consistent with the isotropic moment tensor case, except that the radial rotation component exhibits higher energy, while the equilibration times remain unchanged. For all observables, the two source mechanisms show the remarkably similar dynamics of the envelopes rapidly after the direct arrivals.
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Envelope simulations based on SEM with $ϵ = 25$ per cent in a half-space medium and isotropic source. Comparison of the energy evolution of the radial and transverse components of the displacement (a), normal strain (b) and rotation (c). The vertical black line marks the theoretical P-wave arrival, the red line indicates the S-wave arrival and the blue line corresponds to the Rayleigh waves.</figcaption></figure>

## 5 CONCLUSIONS

In this study, we explore how gradient measurements can complement conventional seismological observables, providing additional perspectives on seismic wave scattering.

Initially, we assess two simulation approaches, the wave equation solved with SEM and the RTE solved with Monte Carlo under the Born approximation, for 10–25 per cent heterogeneity in both full-space and half-space media. The comparison is based on the analysis of scattering mean free path, the transition from multiple scattering to the diffusive regime inferred from energy decay and P/S energy partition. Over 0.6–2.8 Hz, SEM and RTE exhibit excellent consistency. Under these conditions, SEM remains accurate slightly below the nominal frequency limit. Although the mesh-based minimum is 3.1 Hz, the simulations indicate a practical lower bound ≈ 0.3 Hz lower, corresponding to about 1.4 elements per minimum wavelength. The comparison with RTE suggests that the validity range of the Born approximation is broader than commonly assumed, at least for the observables considered, it remains accurate for velocity fluctuations up to $ϵ = 25$ per cent. This arises from the observable-dependent sensitivity to forward scattering: although the Born scheme overemphasizes forward energy at high frequency (geometrical-optics regime), metrics such as diffusivity and equipartition time, which incorporate angular weighting, effectively attenuate these contributions and thus preserve the approximation.

The second part of the study examines the temporal evolution of the energy of displacement, strain and rotation generated by an isotropic moment–tensor source. Although the source does not radiate S-waves directly, near-source P–S conversion produces a ‘ghost’

S source (as in Hirakawa et al. 2016), predominantly polarized along the transverse axis; at the same time, because the scattered-energy flux is anisotropic (Paul et al. 2005), radial rotations remain weak. Analysis of the energy envelopes at a central frequency of 2 Hz ($k_p a = 1.3$) and a source–receiver distance of 1.5 P-wave mean free paths shows that the diffusive regime is reached about 7 P-mft (≈12 S-mft) after the P arrival, with P–S energy equipartition occurring at ≈30 P-mft (≈50 S-mft). Comparing full-space and half-space simulations, we find that the free surface strongly influences the dynamics through Rayleigh waves after their arrival; prior to the Rayleigh arrival, the envelope evolution matches that of the full-space case. In both configurations, the two horizontal displacement components converge rapidly following the direct S waves, and the strain components display a broadly similar time evolution, except for a systematically smaller vertical (z) component due to the free-surface condition at the observation point. A similar overall pattern is obtained for a dipole source.

The rotation components provide additional insights into the wavefield dynamics and the medium properties. The joint evolution of the three rotational components offers complementary constraints on the source mechanism, which can be particularly useful for complex or volcanic sources. The diffusion onset time, inferred from the convergence of the rotational components, allows estimation of the medium’s diffusivity and hence the diffusive mean free time (or path). Finally, in the equipartition regime, the asymptotic ratio between horizontal and vertical rotational energies is governed by the elastic properties of the medium, providing further information about its behaviour.
## ACKNOWLEDGMENTS

This work is supported by the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie (Grant agreement No 955515 for MB and MC) and by The European Research Council under the European Union Horizon 2020 research and innovation program (Grant agreements no 742335, F-IMAGE for MB, MC and LM).

The authors thank the two anonymous reviewers and the editor for their helpful comments and suggestions, which significantly improved the manuscript.

**Authors' contribution:**

*Conceptualization*—MC, LM, MB; *Data curation*—MB; *Formal analysis*—MB, MC, LM; *Funding acquisition*—MC; *Investigation*—MB, LM, MC, RM; *Methodology*—MB, MC, LM; *Project administration*—MC; *Resources*—LM; *Software*—RB, LM; *Supervision*—MC, LM; *Validation*—LM, MC, RB; *Visualization*—MB; *Writing—original draft*—MB; *Writing—review and editing*—MC, RB, LM, MB.

## SUPPORTING INFORMATION

Supplementary data are available at *GJIRAS* online.

**Figure S1.** Autocorrelation of the medium (a) and its spectrum (b), fitted with a Von Kármán function, characterized by $k = 2$ and a correlation length of 320 m.

**Figure S2.** Comparison between the diffusion constant obtained from the radiative transfer equation ($D_{RTE}$) and that derived from the elastic wave equation ($D_{WE}$) by fitting the diffusion equation. Case study: full space with $\epsilon = 25$ per cent, evaluated at 14 km from the source.

**Figure S3.** Mean free times for the two case studies: $\epsilon = 10$ per cent (a) and $\epsilon = 25$ per cent (b). $t_p$ and $t_s$ represent the average times that P or S seismic waves interact with a heterogeneity before being polarized into S or P modes, respectively. $t_{ps}$ and $t_{sp}$ are the average times after which, through interaction with a heterogeneity, P waves convert to S waves and S waves convert to P waves, respectively. In the case of $\epsilon = 25$ per cent, it can be observed that $t_p v_p$, above 3 Hz, would become smaller than the correlation length. This is physically impossible and arises from the fact that the mean free time calculations were performed using the Born approximation, which becomes unreliable under strong scattering conditions at high frequencies.

**Figure S4.** Case study in a half-space medium with $\epsilon = 25$ per cent and an isotropic moment tensor. The figure displays the energy partition ratio of the kinetic energy of the three rotational components, defined as the sum of the horizontal components divided by the vertical component.

**Figure S5.** Case study in full space medium with $\epsilon = 17$ per cent and isotropic moment tensor. Comparison of the energy evolution of the radial and transverse components of the displacement (a), linear strain (b) and rotation (c). The vertical black line marks the theoretical P-wave arrival and the red line indicates the S-wave arrival.

**Figure S6.** Case study in half-space medium with $\epsilon = 10$ per cent and $M_{zz} = 1$. Comparison of the energy evolution of the radial and transverse components of the displacement (a), linear strain (b) and rotation (c). The vertical black line marks the theoretical P-wave arrival, the red line indicates the S-wave arrival and the blue line corresponds to the Rayleigh waves.

**Table S1.** Settings of the simulations of the elastic wave equation.

## DATA AVAILABILITY

Seismograms generated by wave equation simulations, along with the results of the RTE for the full-space media, are available at: Bracale *et al.* (2025). Seismograms corresponding to the half-space media are available at: Bracale *et al.* (2025).

## REFERENCES

Aki, K. & Richards, P.G., 2002. *Quantitative Seismology: Theory and Methods*, 2nd edn., University Science Books.

Barabanenkov, Y.N. & Barabanenkov, M.Y., 1997. Radiative transfer theory with time delay for the effect of pulse entrapping in a resonant random medium: general transfer equation and point-like scatterer model, *Waves Random Media*, 7(4), 607.

Barabanenkov, Y.N., Zurk, L. & Barabanenkov, M.Y., 1995. Poynting's theorem and electromagnetic wave multiple scattering in dense media near resonance: modified radiative transfer equation, *J. Electromagnet. Waves Appl.* 9(11-12), 1393-1420.

Barajas, A., Margerin, L. & Campillo, M., 2022. Coupled body and surface wave sensitivity kernels for coda-wave interferometry in a three-dimensional scalar scattering medium, *Geophys. J. Int.*, **230**(2), 1013-1029.

Bernauer, F., Wassermann, J. & Igel, H., 2012. Rotational sensors—a comparison of different sensor types, *J. Seismol.*, **16**(4), 595-602.

Bracale, M., Campillo, M. & Brossier, R., 2025. Numerical simulations in full space heterogeneous media [Dataset]

Bracale, M., Campillo, M. & Brossier, R., 2025. Numerical simulations in half space heterogeneous media [Dataset]

Capdeville, Y. & Sladen, A., 2024. Das sensitivity to heterogeneity scales much smaller than the minimum wavelength, *Seismica*, **3**(1).

Celorio, M., Chaljub, E., Margerin, L. & Stehly, L., 2022. Propagation of 2-D SH waves in random media: Insights from *ab initio* numerical simulations and transport theory, *Front. Earth Sci.* **10**.

Dainty, A., Toksoz, M. *et al.* 1977. Elastic wave propagation in a highly scattering medium—a diffusion approach, *J. Geophys.*, **43**(1), 375-388.

de Hoop, M.V., Garnier, J. & Sølna, K., 2022. System of radiative transfer equations for coupled surface and body waves, *Z. angew. Math. Phys.*, **73**(5), 177.

Emoto, K. & Sato, H., 2018. Statistical characteristics of scattered waves in three-dimensional random media: comparison of the finite difference simulation and statistical methods, *Geophys. J. Int.*, **215**(1), 585-599.

Gaebler, P.J., Sens-Schönfelder, C. & Korn, M., 2015. The influence of crustal scattering on translational and rotational motions in regional and teleseismic coda waves, *Geophys. J. Int.*, **201**(1), 355-371.

Gillet, K., Margerin, L., Calvet, M. & Monnereau, M., 2017. Scattering attenuation profile of the moon: implications for shallow moonquakes and the structure of the megaregolith, *Phys. Earth planet. Inter.*, **262**, 28-40.

Hennino, R., Trégourès, N., Shapiro, N.M., Margerin, L., Campillo, M., van Tiggelen, B.A. & Weaver, R.L., 2001. Observation of equipartition of seismic waves, *Phys. Rev. Lett.*, **86**, 3447-3450.

Hirakawa, E., Pitarka, A. & Mellors, R., 2016. Generation of shear motion from an isotropic explosion source by scattering in heterogeneous media, *Bull. seism. Soc. Am.* **106**(5), 2313-2319.

Hoshiya, M., 1993. Separation of scattering attenuation and intrinsic absorption in japan using the multiple lapse time window analysis of full seismogram envelope, *J. geophys. Res.: Solid Earth*, **98**(B9), 15 809-15 824.

Igel, H., 2017. *Computational Seismology: A Practical Introduction*, Oxford Univ. Press.

Imperatori, W. & Mai, P., 2015. The role of topography and lateral velocity heterogeneities on near-source scattering and ground-motion variability, *Geophys. J. Int.*, **202**(3), 2163-2181.

Khazaie, S., Cottereau, R. & Clouteau, D., 2017. Numerical observation of the equipartition regime in a 3D random elastic medium, and discussion of the limiting parameters, *Comput. Geosci.*, **102**, 56-67.
Komatitsch, D. & Tromp, J., 1999. Introduction to the spectral element method for three-dimensional seismic wave propagation, *Geophys. J. Int.*, **139**(3), 806–822.

Komatitsch, D. & Vilotte, J.P., 1998. The spectral element method: An efficient tool to simulate the seismic response of 2d and 3d geological structures, *Bull. seism. Soc. Am.* **88**(2), 368–392.

Kumagai, H., Saito, T., O'Brien, G. & Yamashina, T., 2011. Characterization of scattered seismic wavefields simulated in heterogeneous media with topography, *J. geophys. Res.: Solid Earth*, **116**(B3).

Lacombe, C., Campillo, M., Paul, A. & Margerin, L., 2003. Separation of intrinsic absorption and scattering attenuation from lg coda decay in Central France using acoustic radiative transfer theory, *Geophys. J. Int.*, **154**(2), 417–425.

Maeda, T., Sato, H. & Nishimura, T., 2008. Synthesis of coda wave envelopes in randomly inhomogeneous elastic media in a half-space: single scattering model including Rayleigh waves, *Geophys. J. Int.*, **172**(1), 130–154.

Margerin, L., 2005. Introduction to radiative transfer of seismic waves, in *Seismic Earth: Array Analysis of Broadband Seismograms, Volume 157*, pp. 229–252, eds Levander, A. & Nolet, G., American Geophysical Union (AGU).

Margerin, L., 2009. Generalized eigenfunctions of layered elastic media and application to diffuse fields, *J. acoust. Soc. Am.*, **125**(1), 164–174.

Margerin, L., 2017. Computation of Green's function of 3-D radiative transport equations for non-isotropic scattering of *P* and unpolarized *S* waves, *Pure appl. Geophys.*, **174**(11), 4057–4075.

Margerin, L., Campillo, M. & Van Tiggelen, B., 2000. Monte Carlo simulation of multiple scattering of elastic waves, *J. geophys. Res.: Solid Earth*, **105**(B4), 7873–7892.

Margerin, L., Bajaras, A. & Campillo, M., 2019. A scalar radiative transfer model including the coupling between surface and body waves, *Geophys. J. Int.*, **219**(2), 1092–1108.

Menina, S. et al., 2023. Stratification of heterogeneity in the lithosphere of mars from envelope modeling of event s1222a and near impacts: Interpretation and implications for very-high-frequency events, *Geophys. Res. Lett.*, **50**(7), e2023GL103202.

Morioka, H., Kumagai, H. & Maeda, T., 2017. Theoretical basis of the amplitude source location method for volcano-seismic signals, *J. geophys. Res.: Solid Earth*, **122**(8), 6538–6551.

Nakahara, H. & Yoshimoto, K., 2011. Radiative transfer of elastic waves in two-dimensional isotropic scattering media: semi-analytical approach for isotropic source radiation, *Earth Planets Space*, **63**(6), 459–468.

Obermann, A., Planès, T., Larose, E., Sens-Schönfelder, C. & Campillo, M., 2013. Depth sensitivity of seismic coda waves to velocity perturbations in an elastic heterogeneous medium, *Geophys. J. Int.*, **194**(1), 372–382.

Obermann, A., Planès, T., Hadziioannou, C. & Campillo, M., 2016. Lapse-time-dependent coda-wave depth sensitivity to local velocity perturbations in 3-D heterogeneous elastic media, *Geophys. J. Int.*, **207**(1), 59–66.

Onodera, K., Kawamura, T., Tanaka, S., Ishihara, Y. & Maeda, T., 2021. Numerical simulation of lunar seismic wave propagation: investigation of subsurface scattering properties near Apollo 12 landing site, *J. geophys. Res.: Planets*, **126**(3), e2020JE006406.

Onodera, K., Maeda, T., Nishida, K., Kawamura, T., Margerin, L., Menina, S., Lognonné, P. & Banerdt, W.B., 2023. Seismic scattering and absorption properties of mars estimated through coda analysis on a long-period surface wave of s1222a marsquake, *Geophys. Res. Lett.*, **50**(13), e2022GL102716.

Parker, T., Shatalin, S. & Farhadiroshan, M., 2014. Distributed acoustic sensing—a new tool for seismic applications, *First Break*, **32**(2), 61–67.

Paul, A., Campillo, M., Margerin, L., Larose, E. & Derode, A., 2005. Empirical synthesis of time-asymmetrical green functions from the correlation of coda waves, *J. geophys. Res.: Solid Earth*, **110**(B8).

Pitarka, A. et al., 2015. Analysis of ground motion from an underground chemical explosion, *Bull. seism. Soc. Am.*, **105**(5), 2390–2410.

Przybilla, J., Korn, M. & Wegler, U., 2006. Radiative transfer of elastic waves versus finite difference simulations in two-dimensional random media, *J. geophys. Res.: Solid Earth*, **111**(B4).

Rodgers, A.J., Petersson, N.A. & Sjogreen, B., 2010. Simulation of topographic effects on seismic waves from shallow explosions near the North Korean nuclear test site with emphasis on shear wave generation, *J. geophys. Res.: Solid Earth*, **115**(B11).

Rytov, S.M., Kravtsov, Y.A. & Tatarskii, V., 1989. *Principles of Statistical Radiophysics: Wave Propagation through Random Media*, Springer.

Ryzhik, L., Papanicolaou, G. & Keller, J.B., 1996. Transport equations for elastic and other waves in random media, *Wave Motion*, **24**(4), 327–370.

Sánchez-Sesma, F.J. & Campillo, M., 1991. Diffraction of *P*, *SV*, and Rayleigh waves by topographic features: a boundary integral formulation, *Bull. seism. Soc. Am.*, **81**(6), 2234–2253.

Sato, H., 1984. Attenuation and envelope formation of three-component seismograms of small local earthquakes in randomly inhomogeneous lithosphere, *J. geophys. Res.: Solid Earth*, **89**(B2), 1221–1241.

Sato, H., 2019. Power spectra of random heterogeneities in the solid earth, *Solid Earth*, **10**(1), 275–292.

Sato, H., Fehler, M.C. & Maeda, T., 2012. *Seismic Wave Propagation and Scattering in the Heterogeneous Earth*, 2nd edn., Springer.

Sens-Schönfelder, C., Margerin, L. & Campillo, M., 2009. Laterally heterogeneous scattering explains lg blockage in the pyrenees, *J. geophys. Res.: Solid Earth*, **114**(B7), e2008JB006107.

Shapiro, B., 1999. New type of intensity correlation in random media, *Phys. Rev. Lett.*, **83**, 4733–4735.

Shapiro, S.A. & Kneib, G., 1993. Seismic attenuation by scattering: Theory and numerical results, *Geophys. J. Int.*, **114**(2), 373–391.

Takemura, S., Furumura, T. & Maeda, T., 2015. Scattering of high-frequency seismic waves caused by irregular surface topography and small-scale velocity inhomogeneity, *Geophys. J. Int.*, **201**(1), 459–474.

Tregoures, N.P. & van Tiggelen, B.A., 2002. Quasi-two-dimensional transfer of elastic waves, *Phys. Rev. E*, **66**(3), 036601.

Trinh, P.T., Brossier, R., Métivier, L., Tavard, L. & Virieux, J., 2019. Efficient time-domain 3D elastic and viscoelastic full-waveform inversion using a spectral-element method on flexible Cartesian-based mesh, *Geophysics*, **84**(1), R75–R89.

Turner, J.A., 1998. Scattering and diffusion of seismic waves, *Bull. seism. Soc. Am.*, **88**(1), 276–283.

Turner, J.A. & Weaver, R.L., 1994a. Radiative transfer and multiple scattering of diffuse ultrasound in polycrystalline media, *J. acoust. Soc. Am.*, **96**(6), 3675–3683.

Turner, J.A. & Weaver, R.L., 1994b. Radiative transfer of ultrasound, *J. acoust. Soc. Am.*, **96**(6), 3654–3674.

van Dinther, C., Margerin, L. & Campillo, M., 2021. Implications of laterally varying scattering properties for subsurface monitoring with coda wave sensitivity kernels: Application to volcanic and fault zone setting, *J. geophys. Res.: Solid Earth*, **126**(12), e2021JB022554.

van Driel, M. et al., 2021. High-frequency seismic events on mars observed by insight, *J. geophys. Res.: Planets*, **126**(2), e2020JE006670.

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

Weaver, R., 1990. Diffusivity of ultrasound in polycrystals, *J. Mech. Phys. Solids*, **38**(1), 55–86.

Xu, Z., Margerin, L. & Mikesell, T.D., 2022. Monte Carlo simulations of coupled body-and Rayleigh-wave multiple scattering in elastic media, *Geophys. J. Int.*, **228**(2), 1213–1236.

Zeng, Y., 2006. Scattered surface wave energy in the seismic coda, *Pure appl. Geophys.*, **163**, 533–548.