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

Elsa Giraudat,1 Arnaud Burtin,2, Arthur Le Ber,1 Mathias Fink,1,
Jean-Christophe Komorowski,2 Alexandre Aubry1
1Institut Langevin, ESPCI Paris, PSL University, CNRS, Paris, France
2Institut de Physique du Globe de Paris, Université Paris Cité, CNRS, Paris,
France

∗To whom correspondence should be addressed; E-mail<EMAIL_ADDRESS>
> 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 it be for focusing (?) or communication
(?) 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 (?, ?) and
optical microscopy (?, ?) 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 (?),
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 seismic noise received at two stations is known
to converge toward the Green’s function between receiving stations (?, ?), as
if one of them had been used as source, thus paving the way to passive matrix
imaging (?, ?, ?). As surface waves dominate ambient noise, most past studies
on the topic aimed at extracting surface wave properties from ambient noise
correlations (?). However, they also contain the contribution of body waves
reflected by deep structures (?) and fluid reservoirs (?).

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 (?, ?). The covariance matrix 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 (?), can then be applied to provide a confocal image
of the subsoil reflectivity (?). 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. However, 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 reversal process driven from the
$\mathbf{k}-$space (?, ?) 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 (?, ?, ?).

### 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.

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
$\mathbf{r}_{\textrm{in/out}}=(\bm{\rho}_{\textrm{in/out}},z)$ and mapping the
inside of the volcano. (C) 1D-velocity model (?) used for the seismic data
redatuming process.

The impulse response $R(i,j,t)$ between each pair of stations $(i,j)$ is
estimated by cross-correlation of ambient seismic noise (?) (section S1). The
set of impulse responses is stored in a time-dependent response matrix
$\mathbf{R}(t)$.

This canonical reflection matrix is powerful since it enables a post-
processing projection of seismic data into different mathematical bases. The
reflection matrix can be investigated into the plane wave basis (or
$\mathbf{k}$-space) or any plane in the real space that sits between the Earth
surface ($\mathbf{u}$) and the expected focal plane ($\bm{\rho}$) 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 (?). Thus, we adopt in
the following a homogeneous 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 (?) (Fig. 1C). This
value ranges from $c_{0}=2000$ m.s-1 at shallow depth to $c_{0}=4300$ m.s-1 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 (?). 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 (?, ?, ?, ?) (section S2).

The result is a focused reflection matrix
$\mathbf{R}_{\bm{\rho\rho}}(z)=[R(\bm{\rho}_{\textrm{out}},\bm{\rho}_{\textrm{in}},z)]$
at each depth $z$ that contains the inter-element impulse responses between a
set of virtual sources at
$\mathbf{r}_{\textrm{in}}=(\bm{\rho}_{\textrm{in}},z)$ and virtual receivers
at $\mathbf{r}_{\textrm{out}}=(\bm{\rho}_{\textrm{out}},z)$ mapping the inner
structure of the volcano (Fig. 1B). Its diagonal elements are associated with
coincident input and output focusing points
($\bm{\rho}_{\textrm{in}}=\bm{\rho}_{\textrm{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 (?);
(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 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 $\delta\rho_{u}$ (Eq. 1) of the theoretical diffraction-
limited focal spot is denoted as a white circle.

### Focusing Quality

The focusing quality can actually be assessed by considering the off-diagonal
elements of $\mathbf{R}_{\bm{\rho\rho}}(z)$
($\bm{\rho}_{\textrm{in}}\neq\bm{\rho}_{\textrm{out}}$) that provide an
estimator of the point spread function in reflection (RPSF, see sections S5
and S6) as a function of the relative position
$\Delta\bm{\rho}=\bm{\rho}_{\textrm{out}}-\bm{\rho}_{\textrm{in}}$ (?, ?).
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 $\mathbf{R}_{\bm{\rho\rho}}(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 $\theta_{u}=\tan^{-1}\left(d_{||}/2z\right)$ under which
the geophone array is seen by the focusing point:

$\delta\rho_{u}=\lambda/\left(2\sin\theta_{u}\right).$ (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 (?).
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 (?) (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).

Figure 3: Exploiting memory effect for overcoming aberrations. (A) Schematic
view of back-scattered wave-fronts in the Earth surface basis ($\mathbf{u}$)
generated by virtual sources ($\mathbf{r}_{\textrm{in}}$) 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
$\mathbf{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 $\delta\rho_{u}$ (Eq. 1) of
the theoretical diffraction-limited focal spot is denoted as a white circle.

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 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

Figure 4: Overcoming diffraction by operating the focusing process from the
$\mathbf{k}-$space. (A) The distorted wave-field in the $\mathbf{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 $\mathbf{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.

Strikingly, an analysis of wave distortions from the $\mathbf{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 (?). In a far-field approximation, the
contribution of each scatterer would emerge onto limited parts of the
$\mathbf{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\bm{\rho}|^{2}/z)$ and resulting from the
curvature of focused wave-fronts (see Fig. 3H and section S5). Projected in
the $\mathbf{k}$-space, the associated transfer 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 $\mathbf{D}-$matrix
expressed in the $\mathbf{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$
(section 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 $\mathbf{k}-$space analysis of the
$\mathbf{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
(?) exhibiting a shift-shift memory effect (?). 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 reflectivity 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.

### Unveiling the plumbing system of La Soufrière

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).

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 strorage system ($z=$5-8.5 km) exhibits a weaker reflectivity probably
due 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 (?, ?, ?).

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), eruptible melt, and magmatic fluids that extend
laterally. This model of a magmatic plumbing system has been described at many
other volcanoes (?, ?). 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 (?, ?, ?, ?).

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 this those determined by
independent petrological studies (?) 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 (?, ?), ultrasound (?, ?)
and passive seismology (?, ?, ?), a novel matrix approach to volcano imaging
is proposed in this paper. Exploiting Green’s functions retrieval from seismic
noise (?, ?), 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 analysis, 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 (?).

## 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, 106769 (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, 107488 (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 seismological and volcanological observatory of Guadeloupe (?, ?). The
seismic noise correlation data generated in this study and codes used to post-
process it are available at Zenodo (?) (https://zenodo.org/record/10066910).

### Supplementary materials

Supplementary Text
Figs. S1-S2-S3-S4
Reference (?, ?, ?, ?, ?)

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 $\mathbf{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 (?) and 6 permanent stations (?) 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 applied the computation of seismic noise correlation functions
(NCFs).

The procedure to compute the seismic NCF mainly follows the stages detailed in
Ref. (?). 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(\mathbf{g}_{j},\mathbf{g}_{i},t)$ between each couple of geophones $i$ and
$j$ whose positions are identified by vector $\mathbf{g}$. The set of the
estimated 2850 vertical impulse responses forms the canonical reflection
matrix
$\mathbf{R}_{\mathbf{g}\mathbf{g}}(t)=[R(\mathbf{g}_{out},\mathbf{g}_{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}_{\mathbf{g}\mathbf{g}}(t)$ to get the set of monochromatic
canonical reflection matrices
$\overline{\mathbf{R}}_{\mathbf{g}\mathbf{g}}(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}(\bm{\rho},\mathbf{g},z,f)=\frac{e^{-i2\pi
f\sqrt{{\|\bm{\rho}-\mathbf{g}_{||}\|}^{2}+{\|z-g_{z}\|}^{2}}/c_{0}(z)}}{4\pi\sqrt{{\|\bm{\rho}-\mathbf{g_{||}}\|}^{2}+{\|z-{g_{z}}\|}^{2}}}.$
(S1)

$\mathbf{G}_{0}(z,f)$ describes the causal 3-D propagation of waves between
any geophone $\mathbf{g}=(g_{x},g_{y},g_{z})=(\mathbf{g_{||}},g_{z})$ and any
focusing point $\bm{\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}}_{\mathbf{g}\mathbf{g}}(f)$ at each depth $z$ is
described by the following matrix product:

$\overline{\mathbf{R}}_{\bm{\rho}\bm{\rho}}(z,f)=\mathbf{G}^{*}_{0}(z,f)\times\overline{\mathbf{R}}_{\mathbf{g}\mathbf{g}}(f)\times\mathbf{G}_{0}^{{\dagger}}(z,f)$
(S2)

or in terms of matrix coefficients :

$\overline{R}(\bm{\rho}_{out},\bm{\rho}_{in},z,f)=\sum_{\mathbf{g}_{out}}G_{0}^{*}(\bm{\rho}_{out},\mathbf{g}_{out},z,f)~{}\sum_{\mathbf{g}_{in}}\overline{R}(\mathbf{g}_{out},\mathbf{g}_{in},f)~{}G_{0}^{*}(\bm{\rho}_{in},\mathbf{g}_{in},z,f)$
(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}}_{\bm{\rho}\bm{\rho}}(z,f)=[\overline{R}(\bm{\rho}_{out},\bm{\rho}_{in},z,f)]$.
Physically, each coefficient of
$\overline{\mathbf{R}}_{\bm{\rho}\bm{\rho}}(z,f)$ contains the inter-element
impulse response between a virtual source located at
$\mathbf{r}_{\textrm{in}}=(\bm{\rho}_{\textrm{in}},z)$ and a virtual detector
at $\mathbf{r}_{\textrm{out}}=(\bm{\rho}_{\textrm{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\mathrm{bar}{z}_{0}$, is only limited by diffraction, such that

$\delta\mathrm{bar}{z}_{0}\sim 2\lambda_{z}/\sin^{2}{\theta_{u}}$ (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}_{\bm{\rho}\bm{\rho}}$ can be derived at each depth by
coherently summing the monochromatic matrices over the frequency bandwidth:

${\mathbf{R}}_{\bm{\rho}\bm{\rho}}(z)=\int_{f^{-}}^{f^{+}}df~{}\overline{\mathbf{R}}_{\bm{\rho}\bm{\rho}}(z,f)$
(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.$ (S6)

In the single scattering regime, the coefficients of
$\mathbf{R}_{\bm{\rho}\bm{\rho}}(z)$ can be theoretically expressed as follows
(?):

${R}(\bm{\rho}_{out},\bm{\rho}_{in},z,f)=\int
d\bm{\rho}{H}(\bm{\rho},\bm{\rho}_{\textrm{out}},z){\gamma}(\bm{\rho},z){H}(\bm{\rho},\bm{\rho}_{\textrm{in}},z)$
(S7)

where $\gamma(\bm{\rho},z)$ is the medium reflectivity at depth $z$.
$H(\bm{\rho},\bm{\rho}_{\textrm{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}_{\textrm{in/out}}$. Its support
defines the characteristic size of each virtual source at $\mathbf{r_{in}}$
and detector at $\mathbf{r_{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 $\bm{\rho}_{c}=\bm{\rho}_{in}=\bm{\rho}_{\textrm{out}}$:

$\mathcal{I}(\bm{\rho}_{c},z)=R(\bm{\rho}_{c},\bm{\rho}_{c},z)$ (S8)

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

$\mathcal{I}(\bm{\rho}_{c},z)=\int
d\bm{\rho}{H}^{2}(\bm{\rho},\bm{\rho}_{c},z){\gamma}(\bm{\rho},z)$ (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

Figure S1: Scattering intensity vs. depth. The logarithm of the mean confocal
intensity, $\langle|I^{\prime}(\bm{\rho},z)^{2}\rangle_{\bm{\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.

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) (?, ?). 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 attenution of the
confocal image is thus mainly due to scattering and absorption. In a
statistically homogeneous disordered medium, the mean intensity,
$\langle|I(\bm{\rho},z)|^{2}\rangle$, shall scale as
$\exp(-2z/\ell_{\textrm{ext}})$. $\ell_{\textrm{ext}}$ is the extinction
length that combines the scattering and absorption losses as follows:
$\ell_{\textrm{ext}}^{-1}=\ell_{\textrm{s}}^{-1}+\ell_{\textrm{a}}^{-1}$, with
$\ell_{\textrm{s}}$, the scattering mean free path and $\ell_{\textrm{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(\bm{\rho},z)|^{2}\rangle_{\bm{\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 intenity is fitted by an exponential curve whose decay provides an
estimation of $\ell_{\textrm{ext}}$ reported in Tab. S1.

Depth range | Extinction length
---|---
0 - 0.5 km | 1665 m
0.5 - 1.9 km | 310 m
1.9 - 3.8 km | 765 m
3.8 -10.5 km | 3070 m
Table S1: Extinction length. Estimation of the extinction length
$\ell_{\textrm{ext}}$ from the depth decay of the confocal intensity displayed
in Fig. S1.

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_{\textrm{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_{\textrm{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_{\textrm{ext}}$
increases drastically: $\ell_{\textrm{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}(\bm{\rho},z)=\exp[\beta(z)/2]\mathcal{I}(\bm{\rho},z),$ (S10)

with $\mathcal{I}_{N}(\bm{\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(\bm{\rho},\bm{\rho}_{\textrm{in/out}},z)$ can be expressed using
the real Green’s function $G(\bm{\rho},\mathbf{g}_{i},z)$ between the
geophones and the focused basis (?):

$H(\bm{\rho},\bm{\rho}_{\textrm{in/out}},z)=\sum_{\mathbf{g}_{i}}G(\bm{\rho},\mathbf{g}_{i},z,f_{0})G_{0}^{*}(\bm{\rho}_{\textrm{in/out}},\mathbf{g}_{i},z,f_{0})$
(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
$\mathbf{g}_{i}$:

$H(\bm{\rho},\bm{\rho}_{\textrm{in/out}},z)=\int
d\mathbf{u}\overline{O}(\mathbf{u})G(\bm{\rho},\mathbf{u},z,f_{0})G_{0}^{*}(\bm{\rho}_{\textrm{in/out}},\mathbf{u},z,f_{0})$
(S12)

with
$\overline{O}(\mathbf{u})=N_{g}^{-1}\sum_{i}\delta(\mathbf{u}-\mathbf{g}_{||,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}(\bm{\rho},\bm{\rho}_{\textrm{in/out}},z)=\exp\left[i\frac{k_{0}}{2z}\left(|\bm{\rho}|^{2}-|\bm{\rho}_{\textrm{in/out}}|^{2}\right)\right]O\left(\frac{\bm{\rho}-\bm{\rho}_{\textrm{in/out}}}{\lambda
z}\right)$ (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 $\overline{O}$, such that
$O\left(\frac{\bm{\rho}-\bm{\rho}_{\textrm{in/out}}}{\lambda z}\right)=\int
d\mathbf{u}\overline{O}(\mathbf{u})\exp\left[i\frac{k_{0}}{z}\mathbf{u}.(\bm{\rho}-\bm{\rho}_{\textrm{in/out}})\right]$.
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)]$ (S14)

with $\theta_{u}$, the mean angle under which the geophone network is seen by
the focusing point $\mathbf{r}_{\textrm{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
$\overline{A}(\mathbf{u},z)$ at the Earth surface, such that

$G(\bm{\rho},\mathbf{u},z,f_{0})=\overline{A}(\mathbf{u},z)G_{0}(\bm{\rho},\mathbf{u},z,f_{0})$
(S15)

Equation S12 then simplifies into:

$H(\bm{\rho},\bm{\rho}_{\textrm{in/out}},z)=\int
d\mathbf{u}\overline{F}(\mathbf{u},z)G_{0}(\bm{\rho},\mathbf{u},z,f_{0})G_{0}^{*}(\bm{\rho}_{\textrm{in/out}},\mathbf{u},z,f_{0})$
(S16)

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

$H(\bm{\rho},\bm{\rho}_{\textrm{in/out}},z)=\exp\left[i\frac{k_{0}}{2z}\left(|\bm{\rho}|^{2}-|\bm{\rho}_{\textrm{in/out}}|^{2}\right)\right]\underbrace{\int
d\mathbf{u}\overline{F}(\mathbf{u},z)\exp\left[i\frac{k_{0}}{z}\mathbf{u}.(\bm{\rho}-\bm{\rho}_{\textrm{in/out}})\right]}_{F\left(\frac{\bm{\rho}-\bm{\rho}_{\textrm{in/out}}}{\lambda
z}\right)}$ (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 $\overline{A}$:

$F\left(\frac{\bm{\rho}-\bm{\rho}_{\textrm{in/out}}}{\lambda
z}\right)=O\otimes A\left(\frac{\bm{\rho}-\bm{\rho}_{\textrm{in/out}}}{\lambda
z},z\right)$ (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. (?) and Touma et al. (?) showed that the
amplitude distribution along each antidiagonal of
$\mathbf{R}_{\bm{\rho}\bm{\rho}}(z)$ provides a key quantity that we will
refer to as the reflection point-spread function (RPSF):

$RPSF(\Delta\bm{\rho},\bm{\rho}_{c},z)=R(\bm{\rho}_{c}-\Delta\bm{\rho},\bm{\rho}_{c}+\Delta\bm{\rho})$
(S19)

Along an antidiagonal of $\mathbf{R}_{\bm{\rho}\bm{\rho}}(z)$, all couple of
points on a given antidiagonal share the same midpoint
$\bm{\rho}_{c}=(\bm{\rho}_{out}+\bm{\rho}_{in})/2$ but with a varying relative
position $\Delta\bm{\rho}=(\bm{\rho}_{out}-\bm{\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 $\bm{\rho}_{s}$ scales as:

$RSPF(\Delta\bm{\rho},\bm{\rho}_{s},z)=\exp\left(i\frac{k_{0}}{z}|\Delta\bm{\rho}|^{2}\right)F\left(\frac{\Delta\bm{\rho}}{\lambda
z}\right)F\left(-\frac{\Delta\bm{\rho}}{\lambda z}\right).$ (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 ($\mathbf{u}$) and focused basis
($\mathbf{r}$).

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

$\mathbf{R}_{\mathbf{u}\bm{\rho}}(z)=\mathbf{G}_{0}^{\top}(z,f_{0})\times\mathbf{R}_{\bm{\rho}\bm{\rho}}(z)$
(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 $\mathbf{R}_{\mathbf{u}\bm{\rho}}$ and its ideal counterpart
$\mathbf{G}_{0}(z)$:

$\mathbf{D}_{\mathbf{u}\bm{\rho}}(z)=\mathbf{G}_{0}^{{\dagger}}(z,f_{0})\circ\mathbf{R}_{\mathbf{u},\bm{\rho}}(z)$
(S22)

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

$D(\mathbf{u}_{out},\bm{\rho}_{in},z)=\underbrace{\overline{F}(\mathbf{u}_{out},z)}_{\mbox{transmittance}}\int
d\bm{\rho}\underbrace{{\gamma}(\bm{\rho},z)\exp\left(i\frac{k_{0}}{z}\left|\bm{\rho}\right|^{2}\right)}_{\mbox{apparent
reflectivity}}\underbrace{F\left(\frac{\bm{\rho}-\bm{\rho}_{\textrm{in}}}{\lambda
z}\right)}_{\mbox{transmit
PSF}}\underbrace{\exp\left(i\frac{k_{0}}{z}\mathbf{u}_{out}.(\bm{\rho}-\bm{\rho}_{\textrm{in}})\right)}_{\mbox{angular
de-scan}}$ (S23)

In previous works, the aberration transmittance was extracted through a
singular value decomposition (SVD) of $\mathbf{D}_{\mathbf{u}\bm{\rho}}$ or,
equivalently, an eigenvalue decomposition (EVD) of the correlation matrix
$\mathbf{C}_{\mathbf{u}\mathbf{u}}=\mathbf{D}_{\mathbf{u}\bm{\rho}}\mathbf{D}_{\mathbf{u}\bm{\rho}}^{{\dagger}}$
(?, ?). 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}_{out},\bm{\rho}_{in},z)=\overline{F}(\mathbf{u}_{out}){\gamma}(\bm{\rho}_{in},z)\exp\left(i{k_{0}}|\bm{\rho}_{in}|^{2}/z\right).$

The $\mathbf{D}$-matrix is then of rank 1 and its first singular vector
$\mathbf{U}_{1}$ directly provides the aberration transmittance at depth $z$.

Physically, the first eigenvector $\mathbf{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 (?). Mathematically, the time
reversal invariant can be found by solving the following iterative relation
(?):

$\lambda^{(n+1)}_{1}\mathbf{U}_{1}^{(n+1)}=\mathbf{C}_{\mathbf{u}\mathbf{u}}\times\mathbf{U}_{1}^{(n)}$
(S24)

with $\mathbf{U}_{1}=\lim\limits_{n\rightarrow\infty}\mathbf{U}_{1}^{(n)}$ and
$\lambda_{1}=\lim\limits_{n\rightarrow\infty}\lambda_{1}^{(n)}$, the first
eigenvalue of $\mathbf{C}_{\mathbf{u}\mathbf{u}}$. ITR converges towards the
wave-front that maximizes the energy backscattered by the virtual guide star.
If this guide star is point like, $\mathbf{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 (?).

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

$\mathbf{W}^{(n+1)}=\exp\left(i\mbox{arg}\left\\{\mathbf{C}\times\mathbf{W}^{(n)}\right\\}\right)$
(S25)

with $\mathbf{W}^{(0)}=[1\cdots 1]^{T}$ chosen arbitrarily as a unit wave-
front. The resulting wave-front,
$\mathbf{W}=\lim\limits_{n\rightarrow\infty}\mathbf{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 $\mathbf{C}_{\mathbf{u}\mathbf{u}}$. In the speckle regime
(?) or in a sparse medium made of a few point-like reflectors, its
coefficients can be expressed as follows:

$C(\mathbf{u}_{out},\mathbf{u}^{\prime}_{out},z)=\overline{F}(\mathbf{u}_{out},z)\overline{F}^{*}(\mathbf{u}^{\prime}_{out},z)\left[\overline{F}\stackrel{{\scriptstyle\mathbf{u}_{out}}}{{\otimes}}\overline{F}\right](\mathbf{u}_{out}-\mathbf{u}^{\prime}_{out},z)$
(S26)

The correlation term of the right hand side results from the Fourier transform
of the input focal spot intensity distribution
$|F\left({\bm{\rho}-\bm{\rho}_{in}}/\lambda z\right)|^{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}_{out},z)=\exp\left(j\mbox{arg}\left\\{\overline{F}(\mathbf{u}_{out},z)\sum_{\mathbf{u}^{\prime}_{out}}W(\mathbf{u}^{\prime}_{out},z)\overline{F}^{*}(\mathbf{u}^{\prime}_{out},z)\left[\overline{F}\stackrel{{\scriptstyle\mathbf{u}_{out}}}{{\otimes}}\overline{F}\right](\mathbf{u}_{out}-\mathbf{u}^{\prime}_{out},z)\right\\}\right).$
(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}_{out},z)=\overline{F}(\mathbf{u}_{out},z).$ (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 $\mathbf{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:

$\mathbf{R}_{\bm{\rho}\bm{\rho}}(z)=\mathbf{G}_{0}^{{\dagger}}(z,f_{0})\times\left[\mathbf{G}_{0}(z,f_{0})\circ\mathbf{W}^{*}(z)\circ\mathbf{D}_{\mathbf{u}\mathbf{\rho}}(z)\right]$
(S29)

The whole process is then iterated to improve the estimation of the aberration
transmittance by alternating aberration correction at input and output (?). 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 $\mathbf{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 acompanying paper, this is nevertheless not sufficient to
have a contrasted image of volcano in depth.

##  S8 Iterative phase reversal driven from the $\mathbf{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
$\mathbf{R}-$matrix can be investigated from the $\mathbf{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:

$\mathbf{R}_{\mathbf{k}\bm{\rho}}(z)=\mathbf{T_{0}}\times\mathbf{R}_{\bm{\rho}\bm{\rho}}(z),$
(S30)

with $\mathbf{T}_{0}$, the Fourier transform operator

$T_{0}(\mathbf{k}_{||},\bm{\rho})=\exp\left(-i\mathbf{k}_{||}.\bm{\rho}\right).$
(S31)

Then, the $\mathbf{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 $\bm{r}_{in}$:

$\mathbf{D}_{\mathbf{k}\bm{\rho}}(z)=\mathbf{R}_{\mathbf{k}\bm{\rho}}(z)\circ\mathbf{T}_{0}^{*}(z)$
(S32)

To derive an expression for the $\mathbf{D}$-matrix coefficients in the
$\mathbf{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 $\overline{F}\equiv\overline{O}$), this
expression writes as follows:

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

with

$K(\mathbf{k}_{out},\bm{\rho},z)=\mathbb{1}_{|\mathbf{k}_{out}|<k_{0}}\left\\{\overline{O}(\bm{\rho}+\lambda
z\mathbf{k}_{out})\stackrel{{\scriptstyle\mathbf{k}_{out}}}{{\otimes}}\exp\left(i\pi\lambda
z|\mathbf{k}_{out}|^{2}\right)\right\\},$ (S34)

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

If we compare Eq. S33 with the $\mathbf{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 $\mathbf{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 $\overline{O}(\bm{\rho}+\lambda z\mathbf{k}_{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(\bm{\rho},z)=\sum_{s}\gamma_{s}\delta(\bm{\rho}-\bm{\rho}_{s}(z))$.
Under this assumption, Eq. S33 becomes:

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

The coefficients of the associated correlation matrix,
$\mathbf{C}_{\mathbf{k}\mathbf{k}}(z)=\mathbf{D}_{\mathbf{k},\bm{\rho}}(z)\mathbf{D}^{{\dagger}}_{\mathbf{k}\bm{\rho}}(z)$,
are given by:

$C(\mathbf{k}_{out},\mathbf{k}^{\prime}_{out},z)=\sum_{s}|\gamma_{s}|^{2}K(\mathbf{k}_{out},\bm{\rho}_{s},z)K^{*}(\mathbf{k^{\prime}}_{out},\bm{\rho}_{s},z)\left[\overline{O}\stackrel{{\scriptstyle\mathbf{k}_{out}}}{{\otimes}}\overline{O}\right](\mathbf{k}_{out}-\mathbf{k}^{\prime}_{out})$
(S36)

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

$C(\mathbf{k}_{out},\mathbf{k}^{\prime}_{out},z)\simeq\sum_{s}|\gamma_{s}|^{2}K(\mathbf{k}_{out},\bm{\rho}_{s})K^{*}(\mathbf{k^{\prime}}_{out},\bm{\rho}_{s})$
(S37)

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

$C(\mathbf{k}_{out},\mathbf{k}^{\prime}_{out})=\sum_{s}\sigma_{s}^{2}U_{s}(\mathbf{k}_{out})U_{s}^{*}(\mathbf{k}^{\prime}_{out})$
(S38)

with $\sigma_{s}$, the singular values of $\mathbf{D}_{\mathbf{k},\bm{\rho}}$,
and $\mathbf{U}_{s}=[U_{s}(\mathbf{k}_{out})]$, the eigenvectors of
$\mathbf{C}_{\mathbf{k}\mathbf{k}}$, or equivalently the output singular
vectors of $\mathbf{D}_{\mathbf{k},\bm{\rho}}(z)$. The identification of Eqs.
S37 and S38 shows a one-to-one association between each eigenstate of
$\mathbf{C}_{\mathbf{k}\mathbf{k}}$ and each scatterer. The singular values
$\sigma_{s}$ are proportional to each scatterer reflectivity $\gamma_{s}$,
while each eigenvector $\mathbf{U}_{s}$ provides the far-field transmittance
$\mathbf{K}(\mathbf{k}_{out},\bm{\rho}_{s})$ of each scatterer.

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 $\mathbf{D_{kr}}$. (B) Phase and (C) modulus of the first
three singular vectors $\mathbf{U}_{s}=[U_{s}(\mathbf{k}_{out})]$. (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 $\mathbf{U}_{s}$ as the focusing law. (G) Phase $\phi$
of the IPR invariant $\mathbf{W}$ and (H) associated angular spectrum obtained
by considering the modulus of $\mathbf{C_{kk}}\times\mathbf{W}$. (I) Confocal
image using the phase conjugate of $\mathbf{W}$ as the focusing law.

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

Note that this observation also enables to revisit the results of Ref. (?)
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
$\mathbf{k}-$space, the resulting focusing law $\mathbf{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 $\mathbf{C_{kk}}\times\mathbf{W}$ (Fig. S2H). The
comparison with the angular spectrum covered by the three first eigenvectors
of $\mathbf{C}_{\mathbf{k},\mathbf{k}}(z)$ (Fig. S2D) shows the benefit of IPR
for tailoring a focusing law operating over the whole $\mathbf{k}$-space. Also
applied at input of the $\mathbf{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 (?). Figure S3C illustrates
this fact by showing the fourth eigenstate $\mathbf{U}_{4}$ of
$\mathbf{D}_{\mathbf{k}\bm{\rho}}$. Its support emerges in the same angular
range as $\mathbf{U}_{2}$. The link between the second and fourth eigenstates
is confirmed by Fig. S3D that displays the phase difference between
$\mathbf{U}_{4}$ and $\mathbf{U}_{2}$. $\mathbf{U}_{4}$ is an higher-order
eigenstate associated with the same reflector as $\mathbf{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 $\mathbf{U}_{4}$ in Fig. S3 correspond to smaller singular
values and give rise to distorted PSFs (?). They cannot be used for imaging
but they pollute the singular value spectrum of
$\mathbf{D}_{\mathbf{k}\bm{\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).

Figure S3: Limit of iterative time reversal processing \- Illustration at
depth $z=6.9$ km. (A) Singular value histogram of the distortion matrix
$\mathbf{D_{kr}}$. (B) Modulus and (C) phase of the second and fourth singular
vectors, $\mathbf{U}_{2}$ and $\mathbf{U}_{4}$. (D) Phase of the Hadamard
product between $\mathbf{U}_{2}$ and $\mathbf{U}_{4}^{*}$.

##  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 $\mathbf{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(\mathbf{k}_{out},\rho_{s})$ in the
$\mathbf{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 frequency
components from $0$ to $k_{0}$. When properly realigned in phase, those
frequency components lead to a corrected PSF whose extension spans over
$\lambda/2$ instead of the usual aperture limited resolution
$\Delta\rho_{u}\sim\lambda z/d_{0}$.

##  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 (?).
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 the
interface between eruptive melt of the magma storage system and the deeper
host-rock.

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.