# Unveiling the deep plumbing system of a volcano by a reflection matrix analysis of seismic noise

Elsa Giraudat, Arnaud Burtin, Arthur Le Ber, Mathias Fink, Jean-Christophe

Komorowski, Alexandre Aubry

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

To cite this version:

Elsa Giraudat, Arnaud Burtin, Arthur Le Ber, Mathias Fink, Jean-Christophe Komorowski, et al. Unveiling the deep plumbing system of a volcano by a reflection matrix analysis of seismic noise. 2023. [hal-04268602]

HAL Id: hal-04268602

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

Preprint submitted on 2 Nov 2023

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Unveiling the deep plumbing system of a volcano by a reflection matrix analysis of seismic noise

Elsa Giraudat,<fnref n="1" /> Arnaud Burtin,<fnref n="2" />, Arthur Le Ber,<fnref n="1" /> Mathias Fink,<fnref n="1" />, Jean-Christophe Komorowski,<fnref n="2" /> Alexandre Aubry<fnref n="1" />

In geophysics, volcanoes are particularly difficult to image because of the multi-scale heterogeneities of fluids and rocks that compose them and their complex non-linear dynamics. By exploiting seismic noise recorded by a sparse array of geophones, we are able to reveal the magmatic and hydrothermal plumbing system of La Soufrière volcano in Guadeloupe. Spatio-temporal cross-correlation of seismic noise actually provides the impulse responses between virtual geophones located inside the volcano. The resulting reflection matrix can be exploited to numerically perform an auto-focus of seismic waves on any reflector of the underground. An unprecedented view on the volcano's inner structure is obtained at a half-wavelength resolution. This innovative observable provides fundamental information for the conceptual modeling and high-resolution monitoring of volcanoes.

In everyday life, a multitude of sensors surround us to monitor our environment. In wave physics, those sensors can be active and work together to control the wave-field at will whether

<footnote n="1">¹Institut Langevin, ESPCI Paris, PSL University, CNRS, Paris, France</footnote>

<footnote n="2">²Institut de Physique du Globe de Paris, Université Paris Cité, CNRS, Paris, France</footnote>

<footnote>*To whom correspondence should be addressed; E-mail: alexandre.aubry@espci.fr.</footnote>
it be for focusing (1) or communication (2) purposes. For imaging, the problem is often ill-
posed because of the medium complexity and/or the sensor array sparsity. This is particularly
the case in seismology, where the topography of the site under investigation can be so irregular
that it is illusory to deploy a large and dense network of geophones.

This paper addresses the issue of seismic imaging in complex areas such as volcanoes or
fault zones based on data recorded by a sparse array of seismometers. The goal is to provide high
spatial resolution and in-depth imaging of such critical areas that are of paramount importance
for Earth sciences. To that aim, we will build on a matrix imaging approach imported from
other fields than geophysics, such as medical ultrasonics (3, 4) and optical microscopy (5, 6)
that were designed for scales ranging from a few centimeters for ultrasonic waves to a few
hundreds of nanometers for light. In contrast with concurrent seismic methods such as full
waveform inversion (7), the strength of matrix imaging lies in the fact that: (i) it does not
rely on a sophisticated wave velocity model whose knowledge is often limited and uncertain in
geophysics; (ii) it is robust with respect to data quality which is a frequent issue in seismology.

Matrix imaging relies on the array response matrix that contains the set of impulse responses
between each seismometer. Although a geophone is purely passive, cross-correlation of seis-
mic noise received at two stations is known to converge toward the Green's function between
receiving stations (8, 9), as if one of them had been used as source, thus paving the way to
passive matrix imaging (10–12). As surface waves dominate ambient noise, most past studies
on the topic aimed at extracting surface wave properties from ambient noise correlations (13).
However, they also contain the contribution of body waves reflected by deep structures (14) and
fluid reservoirs (10).

As a proof-of-concept, we here exploit seismic noise recorded by a sparse geophone network
deployed at the surface of the La Soufrière volcano of Guadeloupe (15, 16). The covariance ma-
trix of this seismic noise provides the reflection matrix that contains all the available information
on the underground reflectivity. A numerical focusing process, often referred to as redatuming
in seismology (17), can then be applied to provide a confocal image of the subsoil reflectiv-
ity (10). This image is directly proportional to the axial fluctuations of the acoustic impedance
associated with length scales typically of the order or smaller than the wavelength. It is therefore
an extremely relevant observable for highlighting the presence of fluid-rock interfaces. How-
ever, the quality of the confocal image is drastically degraded by: (i) the mismatch between the
wave velocity model and its real distribution that gives rise to a foggy image; (ii) the sparsity
and finite size of the geophone network that limit its resolution. The former problem is solved
by compensating wave distortions from the Earth surface, thereby revealing an helical conduit
in the upper part of the volcano. The sparsity issue is addressed by an iterative phase rever-
sal process driven from the k-space (18, 19) that resolves the deep reflectors with a transverse
resolution of the order of a half-wavelength, thereby breaking the free space diffraction-limit
usually limited by the array aperture. The inner structure of the volcano is revealed up to a depth
of 10 km. It shows sub-horizontal bodies linked by thinner sub-vertical structures that match
the current state-of-the-art conceptual and data-derived view of transcrustal magmatic systems.
Such a mush-based model applies to numerous volcanic systems and has indeed been theorized
for La Soufrière of Guadeloupe (20–22).

**Canonical Reflection Matrix**

Figure 1A shows the virtual network of 76 geophones whose distribution has been dictated by
the topography of the volcano. It spans over a lateral extension $d_\| = 1300$ m and a vertical range
$d_z = 500$ m. The impulse response $R(i, j, t)$ between each pair of stations $(i, j)$ is estimated
by cross-correlation of ambient seismic noise (9) (section S1). The set of impulse responses is
stored in a time-dependent response matrix $R(t)$.

This canonical reflection matrix is powerful since it enables a post-processing projection of
<figure><img src="image_3.png" /><figcaption>Figure 1: **Passive imaging of La Soufrière volcano.** (A) Map of the 76 geophones installed above La Soufrière (Guadeloupe, France). Both permanent stations (red) and temporary nodal array (green) are used. (B) Covariance matrix of seismic noise acquired during 2 months is post-processed to obtain the impulse responses between a set of virtual geophones identified by their position r<sub>in/out</sub> = (ρ<sub>in/out</sub>, z) and mapping the inside of the volcano. (C) 1D-velocity model (23) used for the seismic data redatuming process.</figcaption></figure>

seismic data into different mathematical bases. The reflection matrix can be investigated into
the plane wave basis (or k-space) or any plane in the real space that sits between the Earth
surface (u) and the expected focal plane (ρ) at a given time-of-flight t. To project the seismic
data in these latter bases, a wave velocity model is nevertheless required.

**Wave Velocity Model**

As we consider only the vertical component of the impulse responses, collected echoes are
assumed to be mainly associated with P-waves (14). Thus, we adopt in the following a ho-
mogeneous P-wave velocity model. More precisely, for each depth, we define a homogeneous
velocity distribution whose value is calculated on the basis of the four-layer large scale velocity
model (23) (Fig. 1C). This value ranges from c₀ = 2000 m.s⁻¹ at shallow depth to c₀ = 4300
m.s⁻¹ at depth z = 10 km below the surface. The detailed evolution of the wave velocity model
$c_0(z)$ with respect to depth is given in Figure 1C. The assumed background velocity model is
rough but is, as we will see, sufficient to image the volcano by leveraging the matrix approach.

Confocal Redatuming

In a first step, the velocity model is used to back-propagate in-depth the recorded echoes gathered in the canonical reflection matrix $\mathbf{R}$ in order to retrieve local reflectivity information at each depth of interest. Back-propagation is commonly achieved by applying appropriate time delays at emission and at reception to migrate echoes in post-processing. Such focusing operations are frequently used in imaging and are in particular known as redatuming in seismic exploration (24). The matrix formalism offers a convenient framework to easily perform such beamforming in post-processing, especially in the frequency domain where these operations are described using simple matrix products (10–12, 25) (section S2).

The result is a focused reflection matrix $\mathbf{R}_{\rho\rho}(z) = [R(\rho_{\text{out}}, \rho_{\text{in}}, z)]$ at each depth $z$ that contains the inter-element impulse responses between a set of virtual sources at $\mathbf{r}_{\text{in}} = (\rho_{\text{in}}, z)$ and virtual receivers at $\mathbf{r}_{\text{out}} = (\rho_{\text{out}}, z)$ mapping the inner structure of the volcano (Fig. 1B). Its diagonal elements are associated with coincident *input* and *output* focusing points ($\rho_{\text{in}} = \rho_{\text{out}}$, see section S3). After compensation of wave attenuation with depth (section S4), a 3D confocal image of the volcano is obtained (Fig. 2A) with horizontal cross-sections shown for different depths in Fig. 2B : (i) $z = 1.6$ km *i.e* where the most abundant seismicity occurs at La Soufrière; (ii) $z = 6.9$ km *i.e* at the level of the magma reservoir whose depth range is expected between 5.6 and 8.5 km (26); (iii) $z = 9.1$ km *i.e* beyond the magma reservoir.

Whether it be on the transverse or the vertical view (see Supplementary Movies 1, 2 and 3), some scattering structures seem to emerge at different locations in Figs. 2A,B but the overall structure appears to be fully blurred, suggesting a high level of aberrations. Such a raw confocal image is indeed very sensitive to aberrations and its interpretation should be extremely cautious.
<figure><img src="image_4.png" /><figcaption>Figure 2: **Confocal redatuming** (A) Vertical slice of the 3D confocal image along the South-North direction. This image is shown after depth compensation of seismic wave attenuation (section S4). (B) Horizontal slices and (C) associated RPSFs at depths *z* = 1.6 km, 6.9 km and 9.1 km below the summit. The spatial extension *δρ*<sub>u</sub> (Eq. 1) of the theoretical diffraction-limited focal spot is denoted as a white circle.</figcaption></figure>

**Focusing Quality**

The focusing quality can actually be assessed by considering the off-diagonal elements of
**R**<sub>ρρ</sub>(z) (ρ<sub>in</sub> ≠ ρ<sub>out</sub>) that provide an estimator of the point spread function in reflection (RPSF, see sections S5 and S6) as a function of the relative position Δρ = ρ<sub>out</sub> − ρ<sub>in</sub> (11, 12). Figure 2C displays the evolution of RPSF for different depths and highlights a significant spreading of the back-scattered energy over off-diagonal elements of **R**<sub>ρρ</sub>(z). This is a direct manifestation of the gap between the wave velocity model and its real distribution in the volcano. In absence of aberration, all the back-scattered energy would be contained in a diffracted-limited confocal spot (white circle in Figure 2C) whose size is governed by the angle θ<sub>u</sub> = tan<sup>-1</sup>(d<sub>||</sub>/2z) under
which the geophone array is seen by the focusing point:

$$ \delta\rho_u = \lambda / (2 \sin \theta_u). \qquad (1) $$

In Fig. 2C, the focused wave-field spans over a much larger area than this ideal focal spot and significant side lobes appear around the main central lobe, indicating that images suffer from a high level of aberration.

## Overcoming Aberrations

To isolate and compensate for these aberration effects, we build upon a physical phenomenon referred to as the memory effect in wave physics (27). Waves produced by nearby points inside a complex medium can generate highly correlated, but tilted, reflected wave-fields (Fig. 3A). To exploit this *tilt-tilt* memory effect, our strategy is thus the following (28) (section S7): (i) project the reflection matrix between the focused basis ($\mathbf{r}$) and the Earth surface basis ($\mathbf{u}$) (Fig. 3A); (ii) highlight the angular correlations of the reflected wave-field by building a dual-basis matrix (the distortion matrix $\mathbf{D}$) that connects any input focal point in the medium with the distortion exhibited at the Earth surface by the corresponding reflected wavefront (Fig. 3B); (iii) take advantage of the angular correlations between those wave distortions to accurately estimate the aberration phase transmittance in the Earth surface basis through an iterative phase reversal algorithm (Fig. 3C); (iv) phase conjugate the resulting transmittance to tailor adaptive focusing laws that shall compensate for the volcano's heterogeneities (Fig. 3D).

Vertical and horizontal cross-sections of the resulting confocal image are displayed in Figs. 3E and F, respectively. The comparison with the initial image demonstrates the benefit of the correction process, especially at shallow depths ($z < 4$ km) where the twisted conduit of the volcano is revealed. The comparison of the original and the corrected RPSFs (Figs. 2C and 3G) confirms that the focusing quality is significantly improved in this depth range: Whereas the original RPSF (top panel in Fig. 2B) spreads far beyond the diffraction-limited focal spot, the
<figure><img src="image_5.png" /><figcaption>**Figure 3: Exploiting memory effect for overcoming aberrations.** (A) Schematic view of back-scattered wave-fronts in the Earth surface basis (u) generated by virtual sources (r<sub>in</sub>) at a given depth z. (B) The extraction of wave-distortions amounts to an angular de-scan of each input focal spot. (C) Iterative phase reversal applied to the D-matrix at each depth provides an aberration transmittance whose phase is here shown at depths z = 1.6 km, 6.9 km and 9.1 km and whose modulus is encoded as a transparency mask. Black dots indicate the lateral position of geophones. (D) Phase conjugation and tilt of such aberration phase laws enable an adaptive focusing process on each point of the subsoil. (E) Vertical slice of the 3-D confocal image obtained by means of these optimized focusing laws. This image is shown after a depth compensation of seismic wave attenuation (section S4). (F) Corresponding horizontal slices at depths z = 1.6 km, 6.9 km and 9.1 km. (G) Modulus and (H) phase of the resulting RPSFs at the same depths. The spatial extension δρ<sub>u</sub> (Eq. 1) of the theoretical diffraction-limited focal spot is denoted as a white circle.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Overcoming diffraction by operating the focusing process from the k-space. (A) The distorted wave-field in the k-space amounts to a lateral de-scan of each reflected echoes. Under the Franhofer approximation, the far-field projection of each focal spot is limited by the de-scanned geophone network aperture. (B) In the present case, Fresnel diffraction gives rise to a modulation of each focal spot by a parabolic phase law. The diffraction pattern of each reflector then corresponds to Fresnel rings that cover the whole diffraction disk of radius $k_0$. For sake of clarity, only the contribution of a central reflector is displayed. (C) Iterative phase reversal applied to D extracts those diffraction patterns here shown at depths $z = 1.6$ km, 6.9 km and 9.1 km. (D) Vertical slice of the resulting image (same view as in Figs. 2A and 3D). This image is shown after depth compensation of seismic wave attenuation (section S4). (E) Corresponding horizontal slices and (F) RPSFs at the same depths as panel C. The spatial extension $\delta\rho_u$ (Eq. 1) of the aperture-limited focal spot is denoted as a white circle.</figcaption></figure>

transverse extension of the corrected RPSF is drastically reduced. However, the gain in image and focusing quality is more modest at larger depths (Figs. 3E,F). The RPSFs still exhibit secondary lobes, a manifestation of residual aberrations (Fig. 3G). Moreover, the spatial extension of the central lobe is limited by the geophone network aperture (Eq. 1). As a consequence, the deep plumbing system of the volcano, in particular the deepest regions of the transcrustal magmatic system and its magma storage zones beyond 5 km depth, cannot be resolved.

## Beating Diffraction

Strikingly, an analysis of wave distortions from the k-space will allow us to break this fundamental limit. In the plane wave basis, each distorted wave-field corresponds to the diffracted
patterns of each laterally de-scanned output focal spot (4). In a far-field approximation, the
contribution of each scatterer would emerge onto limited parts of the k-space because of the
finite size of the geophone array (see Fig. 4A and section S8). However, the focal spots also
exhibit a parabolic phase law scaling as $\exp(jk_0|\Delta\rho|^2/z)$ and resulting from the curvature of
focused wave-fronts (see Fig. 3H and section S5). Projected in the k-space, the associated trans-
fer function is thus a superposition of Fresnel rings associated with each reflector. The support
of those Fresnel rings is not limited by the geophone network aperture ($k_0 \sin \theta_u$) but covers the
whole diffraction disk of radius $k_0$ (Fig. 4B). Iterative phase reversal applied to the D-matrix
expressed in the k-space leads to a focusing law (Fig. 4C) that realigns the phase of each spatial
frequency component such that the focal spot size reduces to the diffraction limit $\delta\rho_0 \sim \lambda/2$ (sec-
tion S9). It leads to a new confocal image whose several cross-sections are displayed in Figs. 4D
and E. In particular, a complex multi-lens melt reservoir is revealed by Fig 4D beyond a depth
of 5 km. The comparison with the previous image (Figs. 3E and F) highlights the spectacular
gain in terms of contrast and resolution provided by a k-space analysis of the D-matrix. This
observation is confirmed by the new RPSFs displayed in Fig. 4F. Compared to their previous
version (Fig. 3G), the diffuse background has been suppressed by a compensation of residual
wave distortions (28) exhibiting a *shift-shift* memory effect (27). More importantly, the RPSF
extension is now of the order of $\lambda/2$ over the whole considered depth range, thereby beating the
usual aperture-limited resolution (Eq. 1) displayed by conventional imaging methods.

One necessary condition for this striking performance is the sparsity of the volcano reflec-
tivity with only a few reflectors emerging at each depth (Fig. 4E). As the signature of each
reflector is independent, we are able to focus simultaneously on each scatterer provided that
they are not too numerous. More precisely, the contrast of the confocal image will typically
scale as the ratio between the number of independent geophones and the number of reflectors
lying at each depth.
<figure><img src="image_7.png" /><figcaption>Figure 5: **Three-dimensional view of the first 10 km of the hydrothermal and magmatic system of La Soufrière.** (A,B) Isosurface plots of the three-dimensional image of the volcano viewed from East and North, respectively. The isosurface is fixed to be -15 dB. (C,D) Corresponding zooms on the first 3 km depth. The isosurface is fixed to be -10 dB. This image is shown after a depth compensation of seismic wave attenuation (section S4).</figcaption></figure>

## Unveiling the plumbing system of La Soufrière

Figure 5 shows two perpendicular views of La Soufrière down to a depth of 10.5 km below the summit (see also Supplementary Movies 1, 2 and 3). Based on an analysis of the P-wave reflected wave-field, it displays the iso-surfaces of the confocal image obtained at the end of the matrix imaging process. As outlined above, the superior part of the volcano, from a depth of 5 km up to the surface, exhibits the clear signature of a tortuous conduit that finds its way through the hostrock forming the upper part of the volcano. On the contrary, its deep structure, between *ca.* 5 and 8.5 km depth, induces a more diffuse scattering that is compatible with the existence of a vertical succession of several sub-horizontal and irregular globular coalescing structures. Those elements are superimposed over a distance of a few kilometers and linked together by narrow sub-vertical diffuse structures. The sub-horizontal structures extend
laterally over a distance of about 8 km. The presence of the superimposed magma storage
zones is also highlighted by the depth-dependence of unnormalized scattering signal displayed
in Fig. S4. The magma storage system ($z = 5-8.5$ km) exhibits a weaker reflectivity probably
due to the presence of extended magma volumes. The enhancement of the confocal signal above
the outer carapace of the magma storage zone ($z = 3.5-5$ km) may be induced by gases and/or
liquid and/or supercritical hydrothermal and magmatic fluids that are present in the pores of
the hostrock along special zones of elevated porosity-permeability. The increase of reflectivity
observed at the bottom of the magma storage system ($z = 8.5-10$ km) is probably due to a strong
back-reflection at the eruptible melt / hostrock interface. The 3D-view of the internal structure
of La Soufrière volcano displayed in Figure 5 thus constitutes a remarkable advance beyond
the current state-of-the-art because it confirms, for the first time, with great detail and striking
similarity the typical structure of transcrustal magmatic systems below volcanoes that has been
predicted by previous conceptual and petrological models (29–31).

Transcrustal magmatic systems consist of vertically-arranged piles of lenses of magmatic
mushes more or less ductile (intricate network of crystals and intersticial melt fraction), erupt-
ible melt, and magmatic fluids that extend laterally. This model of a magmatic plumbing system
has been described at many other volcanoes (32,33). The internal image of the volcano revealed
by Fig. 5 strikingly matches the complex structure described by recent studies on La Soufrière
of Guadeloupe (21, 22, 34, 35).

Last but not least, the seismic confocal image of La Soufrière shows that the main magmatic
plumbing systems extends from about 5 km below the surface to a depth of at about 8.5 km,
values in agreement with those determined by independent petrological studies (26) who
showed that, for the last magmatic eruption of La Soufrière in 1530 CE, the top of the magma
storage zone was located between 5.6 and 7 km and the base could not exceed 8.5 km.
Discussion

Inspired by pioneering works in optical microscopy (6, 18), ultrasound (4, 19) and passive seis-
mology (10–12), a novel matrix approach to volcano imaging is proposed in this paper. Ex-
ploiting Green’s functions retrieval from seismic noise (8, 9), it takes advantage of bulk seismic
waves reflected by heterogeneities in depth to map the internal structure of La Soufrière of
Guadeloupe. This volcano is currently in a state of significant long-lasting unrest that could
result from dynamic changes in the magmatic plumbing system. Such changes can propagate
to the shallow hydrothermal system and could herald eruptive unrest in the future.

The seismic matrix image of the magmatic plumbing system reveals, for the first time,
high-resolution features of the magma storage zone, its geometry and dimensions, its complex
layered structure, its relative connectedness with other regions of the multi-layer transcrustal
magmatic system, and the size and geometry of the upper final eruptive conduit.

The impedance contrast in this complex image also offers the potential, upon further anal-
ysis, to distinguish zones of mush from those of eruptible melt, their relative volume, their
position in the system. Hence, it can lead to the estimation of parameters such as pressure,
temperature, volatile saturation, density contrast, and the connectivity to the surface in evolving
magmatic systems, parameters that drive volcanic eruptions.

The strength of this new imaging method lies in its robustness with respect to sparsity of the
geophone array and inaccuracy of the wave velocity model. In the future, it will be combined
with time-lapse ability resulting from reiteration surveys at active unresting volcanoes and can
be coupled with multiparameter data analysis from other classic monitoring networks. Matrix
imaging can therefore become a revolutionary game changer in the way scientists understand
and model volcanic systems and how volcano observatories monitor their evolving dynamics
to forecast their potential for hazardous eruptive activity that threatens the lives of 800 million
people living within 100 km from a dangerous volcano (36).
References

1. A. P. Mosk, A. Lagendijk, G. Lerosey, M. Fink, *Nat. Photonics* **6**, 283 (2012).

2. A. L. Moustakas, H. U. Baranger, L. Balents, A. M. Sengupta, S. H. Simon, *Science* **287**, 287 (2000).

3. G. Chau, M. Jakovljevic, R. Lavarello, J. Dahl, *Ultrason. Imaging* **41**, 3 (2019).

4. W. Lambert, L. A. Cobus, T. Frappart, M. Fink, A. Aubry, *Proc. Nat. Acad. Sci. USA* **117**, 14645 (2020).

5. S. Kang, et al., *Nat. Commun.* **8**, 2157 (2017).

6. A. Badon, et al., *Sci. Adv.* **6**, eaay7170 (2020).

7. J. Virieux, S. Operto, *Geophysics* **74**, 1 (2009).

8. R. L. Weaver, O. I. Lobkis, *Phys. Rev. Lett.* **87**, 134301 (2001).

9. M. Campillo, A. Paul, *Science* **299**, 547 (2003).

10. T. Blondel, J. Chaput, A. Derode, M. Campillo, A. Aubry, *J. Geophys. Res.: Solid Earth* **123**, 10936 (2018).

11. R. Touma, T. Blondel, A. Derode, M. Campillo, A. Aubry, *Geophys. J. Int.* **226**, 780 (2021).

12. R. Touma, A. Le Ber, M. Campillo, A. Aubry, *J. Geophys. Res.: Solid Earth* **128**, e2023JB026704 (2023).

13. N. M. Shapiro, M. Campillo, L. Stehly, M. H. Ritzwoller, *Science* **307**, 1615 (2005).

14. P. Roux, K. G. Sabra, P. Gerstoft, W. A. Kuperman, M. C. Fehler, *Geophys. Res. Lett.* **32**, L19303 (2005).
15. Institut De Physique Du Globe De Paris (IPGP), Data collection of the seismological and volcanological observatory of Guadeloupe, https://doi.org/10.18715/guadeloupe.ovsg (2021).

16. A. Burtin, et al., Dense seismic monitoring of La Soufrière de Guadeloupe hydrothermal system, https://www.fdsn.org/networks/detail/ZK_2017 (2018).

17. G. T. Schuster, et al., Seismic Interferometry (Society of Exploration Geophysicists, 2008), pp. 331–448.

18. U. Najar, et al., arXiv: 2303.06119 (2023).

19. F. Bureau, et al., Nat. Commun. 14, 6793 (2023).

20. R. Moretti, et al., J. Volcanol. Geotherm. Res. 393, 1067–69 (2020).

21. A. Metcalfe, et al., Front. Earth Sci. 9 (2021).

22. A. Metcalfe, S. Moune, J.-C. Komorowski, R. Moretti, J. Volcanol. Geotherm. Res. 424, 107–488 (2022).

23. J. Dorel, S. Eschenbrenner, M. Feuillard, Pure Appl. Geophys. 117, 1050 (1979).

24. A. J. Berkhout, C. P. A. Wapenaar, J. Acoust. Soc. Am. 93, 2017 (1993).

25. W. Lambert, L. A. Cobus, M. Couade, M. Fink, A. Aubry, Phys. Rev. X 10, 021048 (2020).

26. M. Pichavant, S. Poussineau, P. Lesne, C. Solaro, J. L. Bourdier, J. Petrol. 59, 257 (2018).

27. G. Osnabrugge, R. Horstmeyer, I. N. Papadopoulos, B. Judkewitz, I. M. Vellekoop, Optica 4, 886 (2017).
28. W. Lambert, L. A. Cobus, J. Robin, M. Fink, A. Aubry, *IEEE Trans. Med. Imag.* **41**, 3921 (2022).

29. K. V. Cashman, R. S. J. Sparks, J. D. Blundy, *Science* **355** (2017).

30. M. Cassidy, M. Manga, K. Cashman, O. Bachmann, *Nat. Commun.* **9** (2018).

31. M. Edmonds, K. V. Cashman, M. Holness, M. Jackson, *Philos. Trans. Royal Soc. A* **377**, 20180298 (2019).

32. M. Edmonds, S. Kohn, E. Hauri, M. Humphreys, M. Cassidy, *Lithos* **252-253**, 216 (2016).

33. I. Koulakov, et al., *Sci. Rep.* **8** (2018).

34. V. Combier, et al., *J. Geophys. Res.: Solid Earth* **120**, 5379 (2015).

35. H. Winslow, et al., *Sci. Rep.* **12**, 9734 (2022).

36. P. Papale, W. Marzocchi, *Science* **363**, 1275 (2019).

37. E. Giraudat, et al., Unveiling the deep plumbing system of a volcano by a reflection matrix analysis of seismic noise [data, code]. Zenodo (2023).

38. S. A. Shapiro, G. Kneib, *Geophys. J. Int.* **114**, 373 (1993).

39. K. Aki, *J. Geophys. Res.* **74**, 615 (1969).

40. G. Bensen, et al., *Geophys. J. Int.* **169**, 1239 (2007).

41. C. Prada, M. Fink, *Wave Motion* **20**, 151 (1994).

42. W. Lambert, J. Robin, L. A. Cobus, M. Fink, A. Aubry, *IEEE Trans. Med. Imag.* **41**, 3907 (2022).
## Acknowledgments

The authors wish to thank R. Touma, M. Campillo and A. Derode for initial discussions on the project; the colleagues at the OVSG-IPGP for field assistance in installing and running the node network, SISMOB network of RESIF (now EPOS-FRANCE); the Parc National de Guadeloupe for permission to install the node in the field.

## Funding

The authors are grateful for the funding provided by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement no. 819261, REMINISCENCE project). This work has also been supported by the AO-IPGP 2017 project "Dense seismic monitoring of the hydrothermal system of La Soufrière de Guadeloupe", the AO-Tellus-INSU 2017 action ALEAS (coord.: A. Burtin), the project "Vers la Plateforme Régionale de Surveillance Tellurique du futur" - (PREST) co-funded by INTERREG Caraïbes V for the European Regional Development Fund, and the European Union's Horizon 2020 research and innovation programme (grant agreement, no. 731070, EUROVOLC project). The authors also thank IPGP for general funding to the Observatoires Volcanologiques et Sismologiques (OVS), the INSU-CNRS for funding provided by the Service National d'Observation en Volcanologie (SNOV), and the Ministère pour la Transition Ecologique et Solidaire (MTES) as well as the IdEx project "Université Paris Cité" (ANR-18-IDEX-0001) for financial support.

## Author Contributions

A.A. and A.B. initiated the project. A.B. designed and conducted the collection of seismic data. A.B. performed the cross-correlation of seismic data. E.G., A.L.B, and A.A. developed the post-processing tools. E.G. and A.A. performed the theoretical analysis. E.G. prepared the
figures. A.B. and J.-C. K. provided the geophysical and volcanological interpretation for the obtained image. E.G., J.-C. K. and A.A. prepared the manuscript. E.G., A.B., A.L.B., J.-C. K., M.F., and A.A. discussed the results and contributed to finalizing the manuscript.

**Competing Interests**

The authors declare no competing interests.

**Data and Materials Availability**

Seismic data used in this manuscript has been deposited at the Data collection of the seismo-logical and volcanological observatory of Guadeloupe (15, 16). The seismic noise correlation data generated in this study and codes used to post-process it are available at Zenodo (37) (https://zenodo.org/record/10066910).

**Supplementary materials**

Supplementary Text

Figs. S1-S2-S3-S4

Reference (38–42)
# Supplementary Information

This document provides further information on: (i) the seismic data and noise correlation processing; (ii) the building of the focused reflection matrix; (iii) the formation of the confocal image; (iv) its depth gain compensation for display purpose; (v) the transmit point spread function; (vi) the reflection point spread function; (vii) the iterative phase reversal algorithm for compensation of wave distortions from the Earth surface basis; (viii) iterative phase reversal driven from k-space and its comparison with a singular value decomposition approach; (ix) the spatial resolution of the final image; (x) the depth evolution of the maximum confocal signal.

## S1 Seismic data and noise correlation processing

The seismic data used in this study consists of a temporary nodal array of 65 geophones (16) and 6 permanent stations (15) operated by the OVSG-IPGP (Volcanologic and Seismologic Observatory of Guadeloupe). The geophone sensors were Zland 3C Gen2 (Fairfieldnodal) with a natural frequency of 5 Hz, recording at 500 samples per second and along 3 orthogonal directions (Vertical, North and East). The 6 OVSG seismic stations are 3 components broadband sensors, all having a flat response in the [1-50] Hz frequency band. The seismic records are sampled at 100 Hz. For this study only vertical components are used. The temporary nodal array was deployed from mid-November 2017 to mid-January 2018 during 2 sessions in order to download seismic data and recharge batteries. Since we moved the location of 5 geophones between both acquisition sessions, we ended with a virtual network of 76 sites, for which we could apply the computation of seismic noise correlation functions (NCFs).

The procedure to compute the seismic NCF mainly follows the stages detailed in Ref. (40). Here, we summarize each step that we apply on seismic recordings whether it was a temporary geophone or a permanent seismic sensor. (1) We detrend each hourly vertical seismic record and removed the mean. (2) We remove the instrument response to homogenize the seismic signals and we applied a band-pass filter between 1 Hz and 20 Hz. (3) We resample the seismic record to a unique sample frequency of 100 Hz. (4) We apply a spectral and temporal normalization by proceeding to a spectral whitening followed by a 1-bit normalization to only keep the sign of the seismic signal. (5) We end with the computation of the NCF by cross-correlating hourly seismic records at each stations pair for time delays ranging from -30 to +30 seconds. To increase the signal to noise ratio of a NCF, we apply some quality checks and a waveform summation by first averaging the 24 hourly NCFs in a daily one, for which we discard hourly segments that were not coherent with the raw daily average (correlation coefficient threshold of 0.5). The average over each daily NCF estimated during the 2 months of nodal array deployment finally provides an estimation of the impulse response $R(g_j, g_i, t)$ between each couple of geophones $i$ and $j$ whose positions are identified by vector $g$. The set of the estimated 2850 vertical impulse
responses forms the canonical reflection matrix $\mathbf{R}_{\text{gg}}(t) = [R(\mathbf{g}_{\text{out}}, \mathbf{g}_{\text{in}}, t)]$ that is used to image the inner structure of La Soufrière volcano in the accompanying paper.

## S2 Broadband focused reflection matrix

To that aim, a temporal Fourier transform is first applied to $\mathbf{R}_{\text{gg}}(t)$ to get the set of monochromatic canonical reflection matrices $\overline{\mathbf{R}}_{\text{gg}}(f)$ over the desired frequency bandwidth [10–20] Hz. The monochromatic matrices are then propagated at emission and at reception towards a focal plane at depth $z$ using the corresponding free-space Green propagator $\mathbf{G}_0(z, f)$:

$$G_0(\boldsymbol{\rho}, \mathbf{g}, z, f) = \frac{e^{-i2\pi f \sqrt{\|\boldsymbol{\rho}-\mathbf{g}\|_2^2 + \|z-g_z\|_2^2}}}{4\pi \sqrt{\|\boldsymbol{\rho}-\mathbf{g}\|_2^2 + \|z-g_z\|_2^2}} \quad (S1)$$

$\mathbf{G}_0(z, f)$ describes the causal 3-D propagation of waves between any geophone $g = (g_x, g_y, g_z) = (\mathbf{g}_{\|}, g_z)$ and any focusing point $\boldsymbol{\rho} = (x, y)$ in the focused basis at depth $z$ in a supposed homogeneous medium with a wave velocity $c_0(z)$. The evolution of the wave velocity $c_0(z)$ with respect to depth is provided in Fig. 1C of the accompanying paper.

Within the framework of matrix imaging, the projection of $\overline{\mathbf{R}}_{\text{gg}}(f)$ at each depth $z$ is described by the following matrix product:

$$\overline{\mathbf{R}}_{\rho\rho}(z, f) = \mathbf{G}_0^*(z, f) \times \overline{\mathbf{R}}_{\text{gg}}(f) \times \mathbf{G}_0^\dagger(z, f) \quad (S2)$$

or in terms of matrix coefficients:

$$\overline{R}(\boldsymbol{\rho}_{\text{out}}, \boldsymbol{\rho}_{\text{in}}, z, f) = \sum_{\mathbf{g}_{\text{out}}} G_0^*(\boldsymbol{\rho}_{\text{out}}, \mathbf{g}_{\text{out}}, z, f) \sum_{\mathbf{g}_{\text{in}}} \overline{R}(\mathbf{g}_{\text{out}}, \mathbf{g}_{\text{in}}, f) G_0^*(\boldsymbol{\rho}_{\text{in}}, \mathbf{g}_{\text{in}}, z, f) \quad (S3)$$

where the symbols $*$, $\dagger$ and $\times$ stand for phase conjugate, transpose conjugate and matrix product respectively. It leads to the set of monochromatic focused reflection matrices $\overline{\mathbf{R}}_{\rho\rho}(z, f) = [\overline{R}(\boldsymbol{\rho}_{\text{out}}, \boldsymbol{\rho}_{\text{in}}, z, f)]$. Physically, each coefficient of $\overline{\mathbf{R}}_{\rho\rho}(z, f)$ contains the inter-element impulse response between a virtual source located at $\mathbf{r}_{\text{in}} = (\boldsymbol{\rho}_{\text{in}}, z)$ and a virtual detector at $\mathbf{r}_{\text{out}} = (\boldsymbol{\rho}_{\text{out}}, z)$ (see Fig. 1B of the accompanying paper). In the single scattering regime an in absence of reverberations, the axial dimension of this focal spot, $\delta\bar{z}_0$, is only limited by diffraction, such that

$$\delta\bar{z}_0 \sim 2\lambda_z / \sin^2 \theta_u \quad (S4)$$

with $\lambda_z = c(z)/f$, the wavelength at depth $z$.

In order to enhance this axial resolution, a broadband focused reflection matrix $\mathbf{R}_{\rho\rho}$ can be derived at each depth by coherently summing the monochromatic matrices over the frequency bandwidth:

$$\mathbf{R}_{\rho\rho}(z) = \int_{f^{-}}^{f^{+}} df \overline{\mathbf{R}}_{\rho\rho}(z, f) \quad (S5)$$
with $f^{\pm} = f_0 \pm \Delta f/2$, $f_0 = 15$ Hz, and $\Delta f = 10$ Hz. The operation amounts to a ballistic time gating of singly-scattered echoes at times $t \sim 2z/c_0(z)$. Thanks to this operation, the axial dimension of virtual geophones is greatly reduced and only limited by the frequency bandwidth $\Delta f$:

$$ \delta z_0 \sim c_0(z)/\Delta f. \qquad (S6) $$

In the single scattering regime, the coefficients of $R_{\rho\rho}(z)$ can be theoretically expressed as follows (25):

$$ R(\rho_{\text{out}}, \rho_{\text{in}}, z, f) = \int d\rho H(\rho, \rho_{\text{out}}, z)\gamma(\rho, z)H(\rho, \rho_{\text{in}}, z) \qquad (S7) $$

where $\gamma(\rho, z)$ is the medium reflectivity at depth $z$. $H(\rho, \rho_{\text{in/out}}, z)$ corresponds to the point-spread-function (PSF), that is to say the spatial amplitude distribution of the focal spot around the focusing point $\mathbf{r}_{\text{in/out}}$. Its support defines the characteristic size of each virtual source at $\mathbf{r}_{\text{in}}$ and detector at $\mathbf{r}_{\text{out}}$.

## S3 Confocal Imaging

A confocal image of the medium can be easily retrieved from the focused reflection matrix at each depth $z$ by considering the diagonal elements which verify $\rho_c = \rho_{\text{in}} = \rho_{\text{out}}$:

$$ \mathcal{I}(\rho_c, z) = R(\rho_c, \rho_c, z) \qquad (S8) $$

Injecting Eq. S7 into the last equation leads to the following expression for the confocal image:

$$ \mathcal{I}(\rho_c, z) = \int d\rho H^2(\rho, \rho_c, z)\gamma(\rho, z) \qquad (S9) $$

Each line of the confocal image results from the sample reflectivity $\gamma$ and the confocal PSF $H^2$. The raw confocal image is displayed in Fig. 2 but note that a time gain compensation is priorly applied to compensate for geometrical spreading as well as scattering and absorption losses that drastically decreases the backscattered energy as a function of depth.

## S4 Depth gain compensation of the confocal image

The raw confocal image exhibits a strong amplitude drop with depth (Fig. S1). This attenuation is due to the decay of energy experienced by seismic waves while they propagate. Without compensation, this attenuation strongly degrades the contrast of the confocal image at large depths.

The depth attenuation of the confocal signal can be caused by several factors such as geometrical spreading, scattering and absorption (intrinsic or anelastic attenuation) (38, 39). In the present case, the geometrical spreading of waves is compensated, at least partially, by the focusing process performed both at input and output of the reflection matrix. The attenuation of the
<figure><img src="image_8.png" /><figcaption>Figure S1: **Scattering intensity vs. depth.** The logarithm of the mean confocal intensity, $\langle |I'(\rho, z)|^2 \rangle_{\rho}$, is plotted as a function of depth. The depth decay of intensity is fitted with an exponential curve over four depth ranges summarized in Tab. S1.</figcaption></figure>

confocal image is thus mainly due to scattering and absorption. In a statistically homogeneous disordered medium, the mean intensity, $\langle |I(\rho, z)|^2 \rangle$, shall scale as $\exp(-2z/\ell_{ext})$. $\ell_{ext}$ is the extinction length that combines the scattering and absorption losses as follows: $\ell_{ext}^{-1} = \ell_s^{-1} + \ell_a^{-1}$, with $\ell_s$, the scattering mean free path and $\ell_a$, the absorption length.

To retrieve such an exponential decay, the random-like fluctuations of the confocal image due to lateral reflectivity variations should be priorly smoothed out by averaging. The resulting mean confocal intensity, $\langle |I(\rho, z)|^2 \rangle_{\rho}$, is displayed in log-scale as a function of effective depth $z$ in Fig. S1. It highlights four depth ranges with distinct decay rates. For each depth range, the decrease of the mean confocal intensity is fitted by an exponential curve whose decay provides an estimation of $\ell_{ext}$ reported in Tab. S1.

<table><caption>Table S1: **Extinction length.** Estimation of the extinction length $\ell_{ext}$ from the depth decay of the confocal intensity displayed in Fig. S1.</caption><thead><tr><td>Depth range</td><td>Extinction length</td></tr></thead><tbody><tr><td>0 - 0.5 km</td><td>1665 m</td></tr><tr><td>0.5 - 1.9 km</td><td>310 m</td></tr><tr><td>1.9 - 3.8 km</td><td>765 m</td></tr><tr><td>3.8 - 10.5 km</td><td>3070 m</td></tr></tbody></table>

Beyond the near-field (0 - 0.5 km) that is probably be polluted by surface waves, the confocal intensity decays rapidly with an exponential fit (red dashed line in Fig. S1) yielding a relatively small value of $\ell_{ext} \sim 310$ m, a manifestation of the strong heterogeneities encountered by the seismic waves in the upper part of the volcano. From depth $z = 1.9$ km to 3.8 km, this decay slows down. As a consequence, the fitting curve (green dashed line in Fig. S1) is associated with a larger value for the extinction length: $\ell_{ext} \sim 760$ m. Beyond $z = 3.8$ km (purple dashed line
in Fig. S1), the mean confocal intensity decays even more slowly and the value of $\ell_{ext}$ increases drastically: $\ell_{ext} \sim 3075$ m. In this region (magma storage zone), there is more back-reflection at interfaces between rocks and fluids but less scattering in average (see Fig. 4D and Fig. 5). This point will be discussed in more details in Sec. S10

The overall fitting curve, $\exp[-\beta(z)]$, displayed in Fig. S1, can be used to normalize at each depth the confocal images shown in the manuscript, such that:

$$ \mathcal{I}_N(\rho, z) = \exp[\beta(z)/2] \mathcal{I}(\rho, z), \qquad (S10) $$

with $\mathcal{I}_N(\rho, z)$, the normalized confocal image displayed in the accompanying manuscript.

Despite our effort to improve the image contrast by compensating attenuation, the raw confocal image displayed in Fig. 2A appears to be fully blurred, suggesting a high level of aberrations. To understand the effects of diffraction and aberrations on this confocal image, the theoretical expression of the transmit PSF is now investigated.

## S5 Transmit point-spread function

The PSF $H(\rho, \rho_{in/out}, z)$ can be expressed using the real Green's function $G(\rho, g_i, z)$ between the geophones and the focused basis (25):

$$ H(\rho, \rho_{in/out}, z) = \sum_{g_i} G(\rho, g_i, z, f_0) G_0^*(\rho_{in/out}, g_i, z, f_0) \qquad (S11) $$

This discrete equation can be rewritten under a continuous form as a function of a coordinate $\mathbf{u}$ describing the Earth surface at a depth origin $z = 0$ defined by the average elevation of the seismic stations $g_i$:

$$ H(\rho, \rho_{in/out}, z) = \int d\mathbf{u}\bar{O}(\mathbf{u})G(\rho, \mathbf{u}, z, f_0)G_0^*(\rho_{in/out}, \mathbf{u}, z, f_0) \qquad (S12) $$

with $\bar{O}(\mathbf{u}) = N_g^{-1} \sum_i \delta(\mathbf{u} - g_{i,i})$, the distribution of geophones, $N_g$, their number and $\delta$ the Dirac distribution.

In the absence of aberrations, i.e if the wave velocity model is valid ($G \equiv G_0$), the expression of the PSF becomes under the Fresnel approximation:

$$ H_0(\rho, \rho_{in/out}, z) = \exp\left[i\frac{k_0}{2z}(|\rho|^2 - |\rho_{in/out}|^2)\right]O\left(\frac{\rho - \rho_{in/out}}{\lambda z}\right) \qquad (S13) $$

In absence of aberration, the reference PSF $H_0$ is thus the product between: (i) a geometrical spreading term; (ii) a parabolic phase law that accounts for the curvature of the focused wave-front; (iii) the Fourier transform $O$ of the geophone network aperture $\bar{O}$, such that $O(\frac{\rho-\rho_{in/out}}{\lambda z}) = \int d\mathbf{u}\bar{O}(\mathbf{u}) \exp[i\frac{k_0}{z}\mathbf{u}(\rho - \rho_{in/out})]$. The transverse dimension of the focal spot, $\delta\rho_u$, is then only limited by diffraction:

$$ \delta\rho_u(z) \sim \lambda z / [2 \sin[\theta_u(z)]] \qquad (S14) $$
with $\theta_u$, the mean angle under which the geophone network is seen by the focusing point $\mathbf{r}_{in/out}$.

In the presence of aberrations, i.e., if the velocity model is inaccurate, there is a mismatch between the true Green's matrix $\mathbf{G}$ and its model $\mathbf{G}_0$. If aberrations are moderate ($\delta\rho < \sqrt{\lambda z}$), they can be accounted for, at each depth $z$, by a phase screen $\bar{A}(\mathbf{u}, z)$ at the Earth surface, such that

$$G(\boldsymbol{\rho}, \mathbf{u}, z, f_0) = \bar{A}(\mathbf{u}, z)G_0(\boldsymbol{\rho}, \mathbf{u}, z, f_0) \quad (S15)$$

Equation S12 then simplifies into:

$$H(\boldsymbol{\rho}, \boldsymbol{\rho}_{in/out}, z) = \int d\mathbf{u}\bar{F}(\mathbf{u}, z)G_0(\boldsymbol{\rho}, \mathbf{u}, z, f_0)G_0^*(\boldsymbol{\rho}_{in/out}, \mathbf{u}, z, f_0) \quad (S16)$$

with $\bar{F}(\mathbf{u}, z) = \bar{O}(\mathbf{u})\bar{A}(\mathbf{u}, z)$, the overall transmittance that combines the array aperture $\bar{O}$ and the aberration phase screen $\bar{A}$. Under the Fresnel approximation, the previous equation becomes:

$$H(\boldsymbol{\rho}, \boldsymbol{\rho}_{in/out}, z) = \exp\left[i\frac{k_0}{2z}(|\boldsymbol{\rho}|^2 - |\boldsymbol{\rho}_{in/out}|^2)\right] \underbrace{\int d\mathbf{u}\bar{F}(\mathbf{u}, z) \exp\left[i\frac{k_0}{z}\mathbf{u} \cdot (\boldsymbol{\rho} - \boldsymbol{\rho}_{in/out})\right]}_{F\left(\frac{\boldsymbol{\rho} - \boldsymbol{\rho}_{in/out}}{\lambda z}\right)} \quad (S17)$$

The PSF $H$ is thus the product of a parabolic phase law that results from the curvature of focused wave-fronts and a focusing function $F$, that results from the convolution between the network PSF $O$ that accounts for diffraction and the aberration PSF $A$ defined as the Fourier transform of the aberration transmittance $\bar{A}$:

$$F\left(\frac{\boldsymbol{\rho} - \boldsymbol{\rho}_{in/out}}{\lambda z}\right) = O \otimes A\left(\frac{\boldsymbol{\rho} - \boldsymbol{\rho}_{in/out}}{\lambda z}, z\right) \quad (S18)$$

where the symbol $\otimes$ stands for the convolution product.

## S6 Reflection point spread function

Interestingly, the focused reflection matrix can provide a local assessment of the focusing quality. Lambert et al. (25) and Touma et al. (11) showed that the amplitude distribution along each antidiagonal of $R_{\rho\rho}(z)$ provides a key quantity that we will refer to as the reflection point-spread function (RPSF):

$$RPSF(\Delta\boldsymbol{\rho}, \boldsymbol{\rho}_c, z) = R(\boldsymbol{\rho}_c - \Delta\boldsymbol{\rho}, \boldsymbol{\rho}_c + \Delta\boldsymbol{\rho}) \quad (S19)$$

Along an antidiagonal of $R_{\rho\rho}(z)$, all couple of points on a given antidiagonal share the same midpoint $\boldsymbol{\rho}_c = (\boldsymbol{\rho}_{out} + \boldsymbol{\rho}_{in})/2$ but with a varying relative position $\Delta\boldsymbol{\rho} = (\boldsymbol{\rho}_{out} - \boldsymbol{\rho}_{in})/2$. This RSPF is a direct indicator of the local focusing quality. For a sparse scattering medium like a volcano, the RPSF at a scatterer position $\boldsymbol{\rho}_s$ scales as:

$$RSPF(\Delta\boldsymbol{\rho}, \boldsymbol{\rho}_s, z) = \exp\left(i\frac{k_0}{z}|\Delta\boldsymbol{\rho}|^2\right) F\left(\frac{\Delta\boldsymbol{\rho}}{\lambda z}\right) F\left(-\frac{\Delta\boldsymbol{\rho}}{\lambda z}\right). \quad (S20)$$
Therefore, the energy spreading in the vicinity of each scatterer position shall enable one to probe the spatial extension of the PSF. As the scatterer position is a priori unknown, the RPSF is, in practice, probed by considering the antidiagonal whose common mid-point exhibits the maximum confocal signal.

## S7 Iterative phase reversal from the Earth surface basis

As highlighted in the previous section, the focused basis is the proper framework for imaging and quantification of focusing quality. However, a dual basis is a better framework to analyse and compensate for aberrations. In the accompanying paper, aberrations are unscrambled by projecting the reflection matrix between the Earth surface (u) and focused basis (r).

The broadband focused reflection matrix can be projected in a geophone basis (here at output) by the following matrix product:

$$R_{u\rho}(z) = G_{0}^{\top}(z, f_{0}) \times R_{\rho\rho}(z) \qquad (S21)$$

where the symbol $\top$ stands for matrix transpose. An angular de-scan of the input focusing points as sketched in Fig. 3B can be performed by a Hadamard product between $R_{u\rho}$ and its ideal counterpart $G_{0}(z)$:

$$D_{u\rho}(z) = G_{0}^{\dagger}(z, f_{0}) \circ R_{u,\rho}(z) \qquad (S22)$$

Each column of the resulting distortion matrix $D_{u\rho}(z)$ maps the phase-distortions with the respect to the ideal wave-front that would be obtained for a point-like source at $r_{in}$. Injecting Eqs. S7 and Eq. S1 leads to the following expression of the D-matrix coefficients under the Fresnel approximation:

$$D(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_{\text{in}}, z) = \underbrace{\overline{F}(\mathbf{u}_{\text{out}}, z)}_{\text{transmittance}} \int d\boldsymbol{\rho} \underbrace{\gamma(\boldsymbol{\rho}, z) \exp\left(i \frac{k_0}{z} |\boldsymbol{\rho}|^2\right)}_{\text{apparent reflectivity}} \underbrace{F\left(\frac{\boldsymbol{\rho} - \boldsymbol{\rho}_{\text{in}}}{\lambda z}\right)}_{\text{transmit PSF}} \underbrace{\exp\left(i \frac{k_0}{z} \mathbf{u}_{\text{out}} \cdot (\boldsymbol{\rho} - \boldsymbol{\rho}_{\text{in}})\right)}_{\text{angular de-scan}} \qquad (S23)$$

In previous works, the aberration transmittance was extracted through a singular value decomposition (SVD) of $D_{u\rho}$ or, equivalently, an eigenvalue decomposition (EVD) of the correlation matrix $C_{uu} = D_{u\rho}D_{u\rho}^{\dagger} (11,28)$. This result can be understood if we assume a point-like transmit PSF in the previous equation. In that case, we have

$$D(\mathbf{u}_{\text{out}}, \boldsymbol{\rho}_{\text{in}}, z) = \overline{F}(\mathbf{u}_{\text{out}}) \gamma(\boldsymbol{\rho}_{\text{in}}, z) \exp\left(ik_0 |\boldsymbol{\rho}_{\text{in}}|^2 / z\right).$$

The D-matrix is then of rank 1 and its first singular vector $U_1$ directly provides the aberration transmittance at depth $z$.

Physically, the first eigenvector $U_1$ is the result of a virtual iterative time reversal (ITR) experiment on a guide star whose reflectivity corresponds to the transmit PSF intensity (28).
Mathematically, the time reversal invariant can be found by solving the following iterative relation (41):

$$ \lambda_1^{(n+1)} U_1^{(n+1)} = C_{uu} \times U_1^{(n)} \quad (S24) $$

with $U_1 = \lim_{n \to \infty} U_1^{(n)}$ and $\lambda_1 = \lim_{n \to \infty} \lambda_1^{(n)}$, the first eigenvalue of $C_{uu}$. ITR converges towards the wave-front that maximizes the energy backscattered by the virtual guide star. If this guide star is point like, $U_1^{(n)}$ thus converges towards the aberration transmittance. However, in reality, the transmit PSF is of course not point-like and its blurring biases the estimation of the aberration transmittance with a time reversal invariant that concentrates on the central part of the geophone array and vanishes on its edge (28).

To circumvent this problem, a related approach consists in an iterative phase reversal (IPR) process (18, 19) that forces a uniform amplitude for the phase reversal invariant $W$, such that

$$ W^{(n+1)} = \exp(i\arg\{C \times W^{(n)}\}) \quad (S25) $$

with $W^{(0)} = [1 \cdots 1]^T$ chosen arbitrarily as a unit wave-front. The resulting wave-front, $W = \lim_{n \to \infty} W^{(n)}$, provides a satisfying estimation of the aberration transmittance over the whole array aperture, as shown in Fig. 3D of the accompanying paper.

Mathematically, the success of IPR can be explained by expressing the correlation matrix $C_{uu}$. In the speckle regime (42) or in a sparse medium made of a few point-like reflectors, its coefficients can be expressed as follows:

$$ C(\mathbf{u}_{\text{out}}, \mathbf{u}'_{\text{out}}, z) = \overline{F}(\mathbf{u}_{\text{out}}, z) \overline{F}^*(\mathbf{u}'_{\text{out}}, z) [\overline{F} \otimes \mathbf{u}_{\text{out}} \otimes \overline{F}] (\mathbf{u}_{\text{out}} - \mathbf{u}'_{\text{out}}, z) \quad (S26) $$

The correlation term of the right hand side results from the Fourier transform of the input focal spot intensity distribution $|F(\rho - \rho_{in}/\lambda z)|^2$. This formulation is reminiscent of the Van Cittert Zernike theorem for an aberrating layer, which links the spatial correlation of a wavefield to the Fourier transform of the intensity distribution from the virtual guide stars (here the input focal spots). In other words, the support of the coherence function scales as the inverse of the input focal spot size. Injecting Eq. S26 into Eq. S25 leads to the following equation:

$$ W(\mathbf{u}_{\text{out}}, z) = \exp\left(j\arg\left\{\overline{F}(\mathbf{u}_{\text{out}}, z) \sum_{\mathbf{u}'_{\text{out}}} W(\mathbf{u}'_{\text{out}}, z) \overline{F}^*(\mathbf{u}'_{\text{out}}, z) [\overline{F} \otimes \mathbf{u}_{\text{out}} \otimes \overline{F}] (\mathbf{u}_{\text{out}} - \mathbf{u}'_{\text{out}}, z)\right\}\right). \quad (S27) $$

For a real autocorrelation function $\overline{F} \otimes \overline{F}$, or equivalently, a symmetric input PSF $|F|^2$, the solution of the previous equation is

$$ W(\mathbf{u}_{\text{out}}, z) = \overline{F}(\mathbf{u}_{\text{out}}, z). \quad (S28) $$

If the previous condition is not fulfilled, the estimation of the aberration transmittance suffers from a bias that can be reduced by iterating the aberration correction process, thereby gradually reducing the size of the virtual guide star and flattening the autocorrelation function $\overline{F} \otimes \overline{F}$.
To do so, the phase conjugate of the estimator W is first used as a focusing law to compensate (partially) for wave distortions. An updated focused reflection matrix is obtained through the following relation:

$$R_{\rho\rho}(z) = G_{0}^{\dagger}(z, f_{0}) \times [G_{0}(z, f_{0}) \circ W^{*}(z) \circ D_{u\rho}(z)] \quad (S29)$$

The whole process is then iterated to improve the estimation of the aberration transmittance by alternating aberration correction at input and output (4). In practice, two iterations of the aberration correction process were enough to converge in the present case.

At the end of the process, a novel confocal image is obtained by considering the diagonal elements of the updated focused reflection matrix [Figs. 3 E and F]. The fine compensation of wave distortions is highlighted by the RSPF (Eq. S20) deduced from the updated focused R-matrix (Fig. 3G). As expected theoretically, compensation of aberrations in the geophone basis enables the recovery of a resolution only limited by the geophone aperture. As explained in the accompanying paper, this is nevertheless not sufficient to have a contrasted image of volcano in depth.

## S8 Iterative phase reversal driven from the k-space

To go beyond, we will now tackle the parabolic phase law exhibited by the transmit PSF (Eq. S17) that was not addressed in the geophone basis and highlighted by the phase of the RPSF in Fig. 3H. To that aim, the R-matrix can be investigated from the k-space.

As previously done for the geophone basis, the first step consists in a projection of the focused reflection matrix in the plane wave basis:

$$R_{k\rho}(z) = T_{0} \times R_{\rho\rho}(z), \quad (S30)$$

with $T_0$, the Fourier transform operator

$$T_0(k_{\parallel}, \rho) = \exp(-i k_{\parallel} \cdot \rho). \quad (S31)$$

Then, the D-matrix is built by isolating the difference between each reflected wave-field and the reference wave-field that would be obtained for a point-like guide star at $r_{in}$:

$$D_{k\rho}(z) = R_{k\rho}(z) \circ T_{0}^{*}(z) \quad (S32)$$

To derive an expression for the D-matrix coefficients in the k-space, one can inject Eqs. S7, S13, S30 and S31 into the last equation. Assuming that aberrations have been fully compensated
in the geophone basis ($F \equiv O$ and $\bar{F} \equiv \bar{O}$), this expression writes as follows:

$$ \mathbf{D}(\mathbf{k}_{\text{out}}, \boldsymbol{\rho}_{\text{in}}, z) = \int d\boldsymbol{\rho} \overbrace{K(\mathbf{k}_{\text{out}}, \boldsymbol{\rho}, z)}^{\text{transmittance}} \overbrace{\gamma(\boldsymbol{\rho}, z) \exp\left(-i \frac{k_0}{2z} |\boldsymbol{\rho}|^2\right)}^{\text{reflectivity}} \underbrace{O\left(\frac{\boldsymbol{\rho} - \boldsymbol{\rho}_{\text{in}}}{\lambda z}\right) \exp\left(i \frac{k_0}{2z} |\boldsymbol{\rho}_{\text{in}}|^2\right)}_{\text{input PSF}} \underbrace{\exp[i\mathbf{k}_{\text{out}} \cdot (\boldsymbol{\rho} - \boldsymbol{\rho}_{\text{in}})]}_{\text{lateral de-scan}}, \quad (\text{S33}) $$

with

$$ K(\mathbf{k}_{\text{out}}, \boldsymbol{\rho}, z) = \mathcal{H}_{|\mathbf{k}_{\text{out}}|<k_0} \left\{ \bar{O}(\boldsymbol{\rho} + \lambda z \mathbf{k}_{\text{out}}) \mathbf{k}_{\text{out}} \otimes \exp(i\pi\lambda z |\mathbf{k}_{\text{out}}|^2) \right\}, \quad (\text{S34}) $$

the aperture transmittance projected in the k-space. The symbol $\mathcal{H}_{|\mathbf{k}_{\text{out}}|<k_0}$ accounts for the low pass filter operated by diffraction in the spatial frequency domain: Only spatial frequencies $\mathbf{k}_{\text{out}}$ whose magnitudes are smaller than $|k_0|$ can propagate into the Earth; higher spatial frequency components are evanescent and cannot probe in depth the Earth beyond a wavelength.

If we compare Eq. S33 with the D-matrix in the geophone basis (Eq. S23), the main difference lies in the fact that the aperture transmittance is no longer isoplanatic from the k-space. Indeed, the angular component of the wave-field recorded by the geophone network and induced by one scatterer in the field-of-view strongly depends on its position, as shown by the term $\bar{O}(\boldsymbol{\rho} + \lambda z \mathbf{k}_{\text{out}})$ in Eq. S34 and highlighted by Fig. 4A. While this would be an issue in the speckle regime (random reflectivity), this property can become an asset in a sparse scattering medium made of a few reflectors: $\gamma(\boldsymbol{\rho}, z) = \sum_s \gamma_s \delta(\boldsymbol{\rho} - \boldsymbol{\rho}_s(z))$. Under this assumption, Eq. S33 becomes:

$$ \mathbf{D}(\mathbf{k}_{\text{out}}, \boldsymbol{\rho}_{\text{in}}, z) = \sum_{s} K(\mathbf{k}_{\text{out}}, \boldsymbol{\rho}_{s}, z) \gamma_{s} \exp\left(-i \frac{k_{0}}{2z} |\boldsymbol{\rho}_{s}(z)|^{2}\right) O\left(\frac{\boldsymbol{\rho}_{s}(z) - \boldsymbol{\rho}_{\text{in}}}{\lambda z}\right) \exp\left(i \frac{k_{0}}{2z} |\boldsymbol{\rho}_{\text{in}}|^{2}\right) \exp[i\mathbf{k}_{\text{out}} \cdot (\boldsymbol{\rho}_{s}(z) - \boldsymbol{\rho}_{\text{in}})]. \quad (\text{S35}) $$

The coefficients of the associated correlation matrix, $C_{kk}(z) = D_{k,\rho}(z)D_{k\rho}^{\dagger}(z)$, are given by:

$$ C(\mathbf{k}_{\text{out}}, \mathbf{k}'_{\text{out}}, z) = \sum_{s} |\gamma_s|^2 K(\mathbf{k}_{\text{out}}, \boldsymbol{\rho}_s, z) K^*(\mathbf{k}'_{\text{out}}, \boldsymbol{\rho}_s, z) \left[ \overline{O} \otimes^{\mathbf{k}_{\text{out}} \otimes \overline{O}} \right] (\mathbf{k}_{\text{out}} - \mathbf{k}'_{\text{out}}) \quad (\text{S36}) $$

For analytical tractability, we will consider, in first approximation, the correlation term $\overline{O} \otimes^{\mathbf{k}_{\text{out}} \otimes \overline{O}}$ as constant:

$$ C(\mathbf{k}_{\text{out}}, \mathbf{k}'_{\text{out}}, z) \approx \sum_{s} |\gamma_s|^2 K(\mathbf{k}_{\text{out}}, \boldsymbol{\rho}_s) K^*(\mathbf{k}'_{\text{out}}, \boldsymbol{\rho}_s) \quad (\text{S37}) $$

Under this assumption and provided that the scatterers belong to different resolution cells ($K^\dagger K \propto I$, with I the identity matrix), Eq. S36 has the form of an eigenvalue decomposition of $C_{kk}$:

$$ C(\mathbf{k}_{\text{out}}, \mathbf{k}'_{\text{out}}) = \sum_{s} \sigma_{s}^{2} U_{s}(\mathbf{k}_{\text{out}}) U_{s}^{*}(\mathbf{k}'_{\text{out}}) \quad (\text{S38}) $$
<figure><img src="image_9.png" /><figcaption>Figure S2: Comparison between iterative time reversal (left) and iterative phase reversal (right) processes for aberration phase law extraction - Illustration at depth *z* = 6.9 km. (A) Singular value histogram of the distortion matrix D<sub>kr</sub>. (B) Phase and (C) modulus of the first three singular vectors U<sub>s</sub> = [U<sub>s</sub>(k<sub>out</sub>)]. (D) Sum of the modulus of the first three eigenvectors. (E) Final image resulting from the combination of each image displayed in panel F. (F) Image obtained using the phase conjugate of U<sub>s</sub> as the focusing law. (G) Phase φ of the IPR invariant W and (H) associated angular spectrum obtained by considering the modulus of C<sub>kk</sub> × W. (I) Confocal image using the phase conjugate of W as the focusing law.</figcaption></figure>

with $\sigma_s$, the singular values of $D_{k,\rho}$, and $U_s = [U_s(k_{out})]$, the eigenvectors of $C_{kk}$, or equivalently the output singular vectors of $D_{k,\rho}(z)$. The identification of Eqs. S37 and S38 shows a one-to-one association between each eigenstate of $C_{kk}$ and each scatterer. The singular values $\sigma_s$ are proportional to each scatterer reflectivity $\gamma_s$, while each eigenvector $U_s$ provides the far-field transmittance $K(k_{out}, \rho_s)$ of each scatterer.

Figure S2 confirms this conjecture by showing the result of the ITR processing applied to $D_{k\rho}(z)$ at depth $z = 6.9$ km. As seen before, this process is mathematically equivalent to the SVD of $D_{k\rho}(z)$. Figure S2A shows the singular value spectrum of $D_{k\rho}(z)$ dominated by three singular values eigenvalues. Each corresponding eigenvector $U_s$ covers a distinct part of the k-space (Fig. S2D), as predicted by the term $\bar{O}(\rho_s + \lambda z k_{out})$ in Eq. S34. The phase of these eigenvectors (Fig. S2B) shows the Fresnel rings corresponding to the parabolic phase term $\exp(i\pi\lambda z|k_{out}|^2)$ in Eq. S34.

Note that this observation also enables to revisit the results of Ref. (11) that showed a similar
feature in the fault area of San Jacinto. As shown in that previous paper, the phase conjugate of each eigenstate can provide the focusing law to image each scatterer (Figure S2F). A compound image can be built by combining the result provided by each eigenstate. However, this approach only provides an image of the main structures of the volcano. As we will see further, it actually fails in highlighting smaller reflectors.

The IPR process that we previously introduced above in the geophone basis can provide a much more complete view of the inner volcano structure. By forcing a transmittance estimation with the same weight over the whole k-space, the resulting focusing law $W$ (Fig. S2G) can address simultaneously all scatterers in the field-of-view. The angular spectrum addressed by this focusing law can be estimated by considering the modulus of the vector $C_{kk} \times W$ (Fig. S2H). The comparison with the angular spectrum covered by the three first eigenvectors of $C_{k,k}(z)$ (Fig. S2D) shows the benefit of IPR for tailoring a focusing law operating over the whole k-space. Also applied at input of the R-matrix, the focusing law derived by IPR leads to the confocal image displayed in Fig. S2I. Compared with the ITR process that only focuses on the three main reflectors at the considered depth (Fig. S2E), the IPR algorithm provides a full-field image of the subsoil highlighting six main structures at the same depth (Fig. S2I).

One could argue that this difference comes from the fact that we did not consider enough eigenstates in the ITR process. However, the higher-order eigenstates cannot be used for imaging purposes (11). Figure S3C illustrates this fact by showing the fourth eigenstate $U_4$ of $D_{k\rho}$. Its support emerges in the same angular range as $U_2$. The link between the second and fourth eigenstates is confirmed by Fig. S3D that displays the phase difference between $U_4$ and $U_2$. $U_4$ is an higher-order eigenstate associated with the same reflector as $U_2$. Each reflector gives actually rise to a set of eigenmodes induced by the autocorrelation term in Eq. S36. Only the fundamental modes corresponding to the highest singular values can be considered for imaging. Higher-order eigenmodes as $U_4$ in Fig. S3 correspond to smaller singular values and give rise to distorted PSFs (28). They cannot be used for imaging but they pollute the singular value spectrum of $D_{k\rho}$. Hence they can prevent from imaging scatterers of smaller reflectivity within the framework of an ITR process. The proposed IPR process allows to circumvent this limit by finding a phase reversal invariant over the whole angular spectrum (Fig. S2G).

## S9 Overcoming diffraction to attain a half-wavelength resolution

As highlighted by the final confocal image (Fig.4E and F) and corresponding RPSF (Fig.4G and F), the IPR algorithm driven from the k-space leads to a resolution of the order of $\lambda/2$ much thinner than the usual diffraction limit dictated by the geophone aperture: $\delta\rho_u \sim \lambda z/d_0$. Mathematically, this can be understood by the convolution product between the geophone aperture and Fresnel rings exhibited by the transmittance $K(k_{out}, \rho_s)$ in the k-space (Eq. S34). These Fresnel rings originate from the parabolic phase law exhibited by the focal spots in real space (Eq. S20). Encoded in the secondary lobes of the PSF, this Fresnel phase law exhibits spatial
<figure><img src="image_10.png" /><figcaption>Figure S3: **Limit of iterative time reversal processing** - Illustration at depth *z* = 6.9 km. (A) Singular value histogram of the distortion matrix D<sub>kr</sub>. (B) Modulus and (C) phase of the second and fourth singular vectors, U<sub>2</sub> and U<sub>4</sub>. (D) Phase of the Hadamard product between U<sub>2</sub> and U<sub>4</sub>*.</figcaption></figure>

frequency components from 0 to *k*₀. When properly realigned in phase, those frequency com-
ponents lead to a corrected PSF whose extension spans over λ/2 instead of the usual aperture
limited resolution Δρ*u* ~ λz/d₀.

S10 Depth evolution of the maximum confocal signal

Figure S4 shows the depth evolution of the maximum confocal signal at each depth without normalization. This curve is plotted at the end of the matrix imaging process. The compensation of aberrations and diffraction operated by matrix imaging allows a finer analysis than the preliminary study provided in Sec. S4. It shows different behaviors in each main part of the volcano. While a strong attenuation is observed in the superior part of the volcano (z =0-3.5 km, see also Sec. S4), the deeper part of the volcano (Fig. S4B) exhibits fluctuations around a relatively constant reflectivity. The magma storage zone (z=5-8.5 km) shows a weaker reflectivity probably due the presence of extended magma volumes. This region is surrounded by two areas of larger reflectivity around 4 and 9 km. This larger confocal signal is probably a manifestation of the important impedance mismatch existing at the boundaries of the magma storage zone (fluid-rock interface). Above the outer carapace of the magma storage system, there are lenses of supercritical acid fluids/brines, and then closer to the surface these supercritical fluids become zones with gases and/or liquid hydrothermal fluids that are present in the pores of the host-rocks, along special zones of porosity-permeability (20). On the one hand, this porous region may account for the enhancement of the confocal signal observed between z = 3.5 and 5 km in Fig. S4. On the other hand, the increase of reflectivity at the bottom of the magma storage system (z = 8.5 – 10 km) is a priori due to the back-reflection echo induced by
<figure><img src="image_11.png" /><figcaption>Figure S4: **Depth evolution of the maximum confocal intensity at each depth:** (A) between $z = 1$ and 10.5 km; (B) between $z = 3$ and 10.5 km.</figcaption></figure>

the interface between eruptive melt of the magma storage system and the deeper host-rock.