# Matrix imaging as a tool for high-resolution monitoring of deep volcanic plumbing systems with seismic noise

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

Komorowski, Alexandre Aubry

► To cite this version:

Elsa Giraudat, Arnaud Burtin, Arthur Le Ber, Mathias Fink, Jean-Christophe Komorowski, et al. Matrix imaging as a tool for high-resolution monitoring of deep volcanic plumbing systems with seismic noise. Communications Earth & Environment, 2024, 5 (1), pp.509. 10.1038/s43247-024-01659-2. hal-04703515

HAL Id: hal-04703515

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

Submitted on 23 Sep 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Matrix imaging as a tool for high-resolution monitoring of deep volcanic plumbing systems with seismic noise

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

Check for updates

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

Volcanic eruptions necessitate precise monitoring of magma pressure and inflation for improved forecasting. Understanding deep magma storage is crucial for hazard assessment, yet imaging these systems is challenging due to complex heterogeneities that disrupt standard seismic migration techniques. Here we map the magmatic and hydrothermal system of the La Soufrière volcano in Guadeloupe by analyzing seismic noise data from a sparse geophone array under a matrix formalism. Seismic noise interferometry provides a reflection matrix containing the signature of echoes from deep heterogeneities. Using wave correlations resistant to disorder, matrix imaging successfully unscrambles wave distortions, revealing La Soufrière’s internal structure down to 10 km with 100 m resolution. This method surpasses the diffraction limit imposed by geophone array aperture, providing crucial data for modeling and high-resolution monitoring. We see matrix imaging as a revolutionary tool for understanding volcanic systems and enhancing observatories’ abilities to monitor dynamics and forecast eruptions.

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<fnref n="1" /> 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 prohibitive 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 hundred 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 correlation functions (NCFs)¹³. 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 in reflection is directly proportional to the axial fluctuations of the acoustic impedance associated with length scales 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 true distribution that gives rise to a

<footnote n="1">¹Institut Langevin, ESPCI Paris, PSL University, CNRS, F-75005 Paris, France. ²Université Paris Cité, Institut de Physique du Globe de Paris, CNRS, F-75005 Paris, France. ✉e-mail: alexandre.aubry@espci.fr</footnote>
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's surface, thereby revealing a helical conduit in the upper part of the volcano. The sparsity issue is addressed by an iterative phase reversal process driven from the k-space¹⁹,²⁰ that resolves the deep reflectors with a transverse resolution of the order of a half-wavelength (~100 m), thereby breaking the free space diffraction limit usually limited by the array aperture. This is in contrast with a previous study on the Erebus volcano¹⁰ that provided a reflectivity image of its main structures (lava lake, magma chambers) but at a poor lateral resolution; no compensation for wave distortions was performed. Here, 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²¹⁻²³.

# Results

## Canonical reflection matrix

Figure 1 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_{\parallel} = 1300$ m and a vertical range $d_z = 500$ m. The impulse response $R(g_p, g_p, t)$ between each pair of stations $(i, j)$ (located at positions $g_i$ and $g_p$, respectively) is estimated by cross-correlation of ambient seismic noise⁷ (“Methods”). To avoid the detrimental effect of fumaroles which are extremely localized noise sources (Supplementary Note 2), only the anti-causal component of the NCFs has been considered to estimate the impulse response between each geophone (Supplementary Note 3). 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 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's surface (u) and the expected focal plane ($\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 modeling

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⁻¹ at shallow depth to $c_0 = 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 Fig. 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 $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¹⁰⁻¹²,²⁶ (“Methods”, Eq. (4)).

The result is a focused reflection matrix $R_{pp}(z) = [R(\rho_{out}, \rho_{in}, z)]$ at each depth $z$ (“Methods”, Eq. (5)) that contains the inter-element impulse responses between a set of virtual sources at $r_{in} = (\rho_{in}, z)$ and virtual receivers at $r_{out} = (\rho_{out}, z)$ mapping the inner structure of the volcano (Fig. 1b). Note that those responses are time-gated around the expected ballistic time ($t \sim 2z/c_0$) and thus frequency-averaged over the whole bandwidth. The diagonal elements of $R_{pp}(z)$ are associated with coincident input and output focusing points ($\rho_{in} = \rho_{out}$, see “Methods”). After compensation of wave attenuation with depth (“Methods”, Eq. (11)), 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 Fig. 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.

## Assessing focusing quality

The focusing quality can actually be assessed by considering the off-diagonal elements of $R_{pp}(z)$ ($\rho_{in} \neq \rho_{out}$) that provide an estimator of the point spread

<figure><img src="image_4.png" /><figcaption>**Fig. 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 (orange) are used. Map data: Google, CNES, Airbus, 2023. b Covariance matrix of seismic noise acquired during 2 months is post-</figcaption></figure>

processed to obtain the impulse responses between a set of virtual geophones identified by their position $r_{in/out} = (\rho_{in/out}, z)$ and mapping the inside of the volcano. c 1D-velocity model²⁴ used for the seismic data redatuming process.
<figure><img src="image_5.png" /><figcaption>**Fig. 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 S5). b, c Horizontal slices (b) and associated reflection point spread functions (c) 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>

<figure><img src="image_6.png" /><figcaption>**Fig. 3 | Memory effect.** a Conventional memory effect: a plane wave coming from depth gives rise to a distorted wave-front at the Earth's surface due to wave velocity heterogeneities lying close to the Earth's surface (aberrating layer). If the incident plane wave direction is tilted, the same transmitted wave field is obtained but tilted by the same angle as the incident plane wave. b Tilt memory effect: If we now consider point-like sources at a given depth, an angular memory effect can be</figcaption></figure>

observed through a thin aberrating layer: The transmitted wavefronts display similar wave distortions but are tilted with respect to each other by an amount dictated by the position of point-like sources. c Shift memory effect: In a multi-layered medium, the transmitted wavefronts display similar wave distortions but, this time, laterally shifted with respect to each other.

function in reflection (RPSF) as a function of the relative position $\Delta\rho = \rho_{out} - \rho_{in}$ ("Methods"). Figure 2c displays the evolution of RPSF for different depths and highlights a drastic spreading of the back-scattered energy over off-diagonal elements of $R_{pp}(z)$. This is a direct manifestation of the gap between the wave velocity model and its true distribution in the volcano. Note that an imperfect convergence of the NCFs towards the Green's functions can also lead to an additional incoherent background on the RPSF (Supplementary Note 1). At this stage, it is difficult to discriminate between these different phenomena but we will see that the roughness of the wave velocity model is the main issue in the present case.

In absence of aberration, all the back-scattered energy would be con-
tained in a diffraction-limited confocal spot (white circle in Fig. 2c) whose
size is governed by the angle $\theta_u = \tan^{-1}(d_{\parallel}/2z)$ under which the geo-
phone 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 strong 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²⁸.
The memory effect ensures that a pattern of random phase shifts imparted

to a plane wave-front by an aberrating layer (wave velocity heterogeneity) keeps its general “shape” but is “tilted” if the orientation of the incoming plane wave is also tilted (Fig. 3a). This implies that the distorted wave-front coming from a point-like source inside the medium will be tilted if the source is laterally shifted (Fig. 3b). In the present case, we do not have any real source inside the volcano but virtual sources produced by the reda- tuming process. In a similar way, a shift of the virtual source will imply a tilt of the reflected wave-front at the Earth’s surface as sketched by Fig. 4a.

To exploit this tilt memory effect, our strategy is thus the following (“Methods”): (i) project the reflection matrix between the focused basis (r) and the Earth surface basis (u) (Fig. 4a); (ii) highlight the angular correlations of the reflected wave-field by building a dual-basis matrix (the distortion matrix D) that connects any input focal point in the medium with the distortion exhibited at the Earth surface by the corresponding reflected wavefront (Fig. 4b); (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. 4c); (iv) phase conjugate the resulting transmittance to tailor adaptive focusing laws that shall compensate for the volcano’s heterogeneities (Fig. 4d).

Vertical and horizontal cross-sections of the resulting confocal
image are displayed in Figs. 4e 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
<figure><img src="image_7.png" /><figcaption>**Fig. 4** | 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</figcaption></figure>

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 S5). f Corresponding horizontal slices at depths z = 1.6 km, 6.9 km, and 9.1 km. g, h Modulus (g) 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.

<figure><img src="image_8.png" /><figcaption>**Fig. 5** | 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 echo. Under the Fraunhofer 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₀. For sake of clarity, only the</figcaption></figure>

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 4d). This image is shown after depth compensation of seismic wave attenuation ("Methods", Eq. (11)). e, f Corresponding horizontal slices (e) and RPSFs (f) at the same depths as in (c). The spatial extension δρ<sub>u</sub> (Eq. (1)) of the aperture-limited focal spot is denoted as a white circle.

RPSFs (Figs. 2c and 4g) 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 (Fig. 4e, f). The RPSFs still exhibit secondary lobes, a manifestation of residual aberrations (Fig. 4g). 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¹. 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. 5a and Supplementary Note 9). However, the focal spots also exhibit a parabolic phase law scaling as $exp(jk_0||\Delta\rho||^2/z)$, resulting from the curvature of focused wavefronts (see Fig. 4h and “Methods”). Projected in the 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. 5b). Iterative phase reversal applied to the **D**-matrix expressed in the **k**-space leads to a focusing law (Fig. 5c) 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$ (Supplementary Note 10). It leads to a new confocal image whose several cross-sections are displayed in Fig. 5d and e. In particular, a complex multi-lens melt reservoir is revealed by Fig. 5d beyond a depth of 5 km. The comparison with the previous image (Fig. 4e 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. 5f. Compared to their previous version (Fig. 4g), the diffuse background has been suppressed by a compensation of residual axial wave distortions²⁹ exhibiting a shift memory effect²⁸ (Fig. 3c). More importantly, the RPSF extension is now of the order of $\lambda/2 \sim 100$ m 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. 5e). 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. This explanation shed a new light onto the high-resolution image of the San Jacinto Fault Zone provided by matrix imaging in a previous study¹¹. While super-resolution was initially accounted for by physical phenomena such as multiple scattering or wave channeling, medium sparsity and wave-front curvatures most likely explain the striking performance of matrix imaging in that seismic configuration as well. The phenomenon described above is therefore not just a curiosity of La Soufrière; it can actually be exploited in many other seismic situations.

## Unveiling the plumbing system of La Soufrière

Figure 6 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 upper 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 host rock 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. We interpret them as sub-horizontal lenses containing an unknown volume

of potentially eruptive magma. The presence of the superimposed magma storage zones is also highlighted by the depth-dependence of unnormalized scattering signal displayed in Supplementary Fig. 12. The magma storage 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 host rock 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/host rock interface. The 3D view of the internal structure of La Soufrière volcano displayed in Fig. 6 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³⁰⁻³² and that can be observed in the field in old dissected magmatic systems.

Transcrustal magmatic systems consist of vertically-arranged piles of lenses of magmatic mushes more or less ductile (intricate network of crystals and interstitial 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. 6 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 extend 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²⁷ that 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

Passive seismic matrix imaging 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, inaccuracy of the wave velocity model, presence of spurious arrivals in the NCFs (Supplementary Note 4),

<figure><img src="image_9.png" /><figcaption>Fig. 6 | 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 (“Methods”, Eq. (11)).</figcaption></figure>
and their imperfect convergence towards the Green's functions (Supplementary Note 1). Despite the seemingly random feature of the computed NCFs (Supplementary Note 3), matrix imaging provides a coherent image of the internal structure of the volcano at an unprecedented resolution to the best of our knowledge. This robustness with respect to data quality and migration model inaccuracy is a great asset with respect to more computationally intensive methods such as full waveform inversion.

In the future, matrix imaging will be combined with time-lapse ability resulting from reiteration surveys at active unresting volcanoes and can be coupled with multi-parameter 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³⁷.

## Methods

### Acquisition of seismic data

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 (Fig. 1a), for which we can apply the computation of seismic NCFs.

### Noise correlation processing

The procedure to compute the seismic NCF mainly follows the stages detailed by Bensen et al.³⁸. 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 deconvolve the instrument response to homogenize the seismic signals and we apply 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 station pair for time delays ranging from –30 to +30 s. To increase the signal-to-noise ratio of an 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). Finally, the average over each daily NCF estimated during the 2 months of nodal array deployment, Γ(gᵢ, gⱼ, t), provides an estimation of the impulse response R(gᵢ, gⱼ, t) between each couple of geophones i and j by only considering the anti-causal component to avoid the detrimental effect of fumaroles (Supplementary Note 3): R(gᵢ, gⱼ, t) = Γ(gᵢ, gⱼ, –t) with t spanning from 0 to +30 s. The set of the estimated 2850 vertical impulse responses forms the canonical reflection matrix R<sub>gg</sub>(t) = [R(gᵢ, gⱼ, t)].

### Broadband focused reflection matrix

To build the focused reflection matrix a temporal Fourier transform is first applied to R<sub>gg</sub>(t) to get the set of monochromatic canonical reflection matrices R<sub>gg</sub>(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 G₀(z, f) whose coefficients write:

$$G_0(\rho, \mathbf{g}, z, f) = \frac{e^{-i2\pi f \sqrt{\|\rho - \mathbf{g}\|}^2 + |z - \mathbf{g}_z|^2/c_0(z)}}{4\pi \sqrt{\|\rho - \mathbf{g}\|}^2 + |z - \mathbf{g}_z|^2} \quad (2)$$

G₀(z, f) describes the causal 3-D propagation of waves between any geophone g = (gₓ, gᵧ, g_z) = (g||, g_z) and any focusing point ρ = (x, y) in the focused basis at depth z in a supposed homogeneous medium with a wave velocity c₀(z). The evolution of the wave velocity c₀(z) with respect to depth is provided in Fig. 1c. Within the framework of matrix imaging, the projection of R<sub>gg</sub>(f) at each depth z is described by the following matrix product:

$$\bar{R}_{\rho\rho}(z, f) = G_0^*(z, f) \times \bar{R}_{gg}(f) \times G_0^\dagger(z, f) \quad (3)$$

or in terms of matrix coefficients:

$$\bar{R}(\rho_{\text{out}}, \rho_{\text{in}}, z, f) = \sum_{\mathbf{g}_{\text{out}}} G_0^*(\rho_{\text{out}}, \mathbf{g}_{\text{out}}, z, f) \sum_{\mathbf{g}_{\text{in}}} \bar{R}(\mathbf{g}_{\text{out}}, \mathbf{g}_{\text{in}}, f) G_0^*(\rho_{\text{in}}, \mathbf{g}_{\text{in}}, z, f) \quad (4)$$

where the symbols *, †, and × stand for phase conjugate, transpose conjugate, and matrix product respectively. It leads to the set of monochromatic focused reflection matrices R<sub>ρρ</sub>(z, f) = [R(ρ<sub>out</sub>, ρ<sub>in</sub>, z, f)]. Physically, each coefficient of R<sub>ρρ</sub>(z, f) contains the inter-element impulse response between a virtual source located at r<sub>in</sub> = (ρ<sub>in</sub>, z) and a virtual detector at r<sub>out</sub> = (ρ<sub>out</sub>, z) (see Fig. 1b).

In order to enhance this axial resolution, a broadband focused reflection matrix R<sub>ρρ</sub> can be derived at each depth by coherently summing the monochromatic matrices over the frequency bandwidth:

$$R_{\rho\rho}(z) = \int_{f^{-}}^{f^{+}} df \bar{R}_{\rho\rho}(z, f) \quad (5)$$

with f<sup>±</sup> = f₀ ± Δf/2, f₀ = 15 Hz, and Δf = 10 Hz. The operation amounts to a ballistic time gating of singly-scattered echoes at times t ~ 2z/c₀(z). Thanks to this operation, the axial dimension of virtual geophones is greatly reduced and only limited by the frequency bandwidth Δf: δz₀ ~ c₀(z)/Δf. In the single scattering regime, the coefficients of R<sub>ρρ</sub>(z) can be theoretically expressed as follows²⁶:

$$R(\rho_{\text{out}}, \rho_{\text{in}}, z) = \int d\rho H(\rho, \rho_{\text{out}}, z) \gamma(\rho, z) H(\rho, \rho_{\text{in}}, z) \quad (6)$$

where γ(ρ, z) is the medium reflectivity at depth z. H(ρ, ρ<sub>in/out</sub>, z) corresponds to the point-spread-function (PSF), that is to say, the spatial amplitude distribution of the focal spot around the focusing point r<sub>in/out</sub>. Its support defines the characteristic size of each virtual source at r<sub>in</sub> = (ρ<sub>in</sub>, z) and detector at r<sub>out</sub> = (ρ<sub>out</sub>, z).

Under the Fresnel approximation, the transmit PSF H is the product of a parabolic phase law that results from the curvature of focused wave fronts and a focusing function F (Supplementary Note 6):

$$H(\rho, \rho_{\text{in/out}}, z) = \frac{1}{z^2} \exp \left[ -i \frac{k_0}{2z} \left( \| \rho \|^2 - \| \rho_{\text{in/out}} \|^2 \right) \right] F \left( \frac{\rho - \rho_{\text{in/out}}}{\lambda z} \right) \quad (7)$$

The focusing function F results from the convolution between the network PSF O that accounts for diffraction and an aberration PSF A that results from the mismatch between the wave velocity model and the true wave velocity distribution in the volcano (Supplementary Note 6):

$$F\left(\frac{\rho - \rho_{\text{in/out}}}{\lambda z}\right) = O \otimes A\left(\frac{\rho - \rho_{\text{in/out}}}{\lambda z}, z\right) \quad (8)$$

where the symbol ⊗ stands for the convolution product.
<table><caption>**Table 1 | Extinction length**</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>

Estimation of the extinction length $\ell_{ext}$ from the depth decay of the confocal intensity displayed in Supplementary Fig. 9.

## Confocal Imaging

The 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_{in} = \rho_{out}$:

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

At each depth $z$, each line of the confocal image results from the convolution between sample reflectivity $\gamma$ and the confocal PSF $H^2$ (Eqs. (6) and (9)):

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

The confocal image is displayed in Fig. 2a but note that a time gain compensation has been priorly applied to get a homogeneous contrast over the whole depth range, as described below.

## Depth gain compensation of the confocal image

The raw confocal image actually exhibits a strong amplitude drop with depth (Supplementary Fig. 9). 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 inelastic attenuation)$^{39,40}$. 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 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 Supplementary Fig. 9. 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 Table 1.

The overall fitting curve, $\exp[-\beta(z)]$, displayed in Supplementary Fig. 9, 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 (11) $$

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

## Reflection point spread function

Interestingly, the focused reflection matrix can provide a local assessment of the focusing quality. Lambert et al.$^{26}$ showed that the amplitude distribution along each antidiagonal of $R_{pp}(z)$ provides a key quantity that we will refer to

as the reflection point-spread function (RPSF):

$$ RPSF(\Delta\rho, \rho_c, z) = R(\rho_c - \Delta\rho, \rho_c + \Delta\rho, z) \qquad (12) $$

Along an antidiagonal of $R_{pp}(z)$, all couple of points on a given antidiagonal share the same midpoint $\rho_c = (\rho_{out} + \rho_{in})/2$ but with a varying relative position $\Delta\rho = (\rho_{out} - \rho_{in})/2$. In the vicinity of an isolated scatterer at $(\rho_s, z)$, the RSPF is a direct indicator of the local focusing quality (Supplementary Note 7):

$$ RSPF(\Delta\rho, \rho_s, z) = \exp\left(i\frac{k_0}{z} \|\Delta\rho\|^2\right) F\left(\frac{\Delta\rho}{\lambda z}\right) F\left(-\frac{\Delta\rho}{\lambda z}\right). \qquad (13) $$

Therefore, the energy spreading in the vicinity of each scatterer position shall enable one to probe the spatial extension of the transmit PSF. As the scatterer positions are a priori unknown, the RSPF is, in practice, probed by considering the antidiagonal whose common mid-point exhibits the maximum confocal signal.

## Compensation of aberrations

As highlighted above, the focused basis is the adequate framework for imaging and quantification of focusing quality. However, the reflection matrix shall be investigated into a dual basis to analyse and compensate for aberrations.

The relevant basis for aberration correction depends on the nature of the medium heterogeneities. In a statistically random medium, the most adequate correction plane lies at $z/3$, a plane where the tilt and shift memory effect can be combined to maximize isoplanicity$^{28}$. For a medium displaying wave velocity heterogeneities localized in depth, the best correction plane is the one conjugated with this aberrating layer. In the present case, we are in the latter situation with wave velocity heterogeneities arising close to the Earth's surface.

Therefore, the most adequate basis is here a plane lying on the Earth surface, i.e., at the depth origin $z=0$ defined by the average elevation of the seismic stations $g_i$. This plane is described by the coordinate vector $\mathbf{u}$. The broadband matrix can be projected in the Earth's surface basis, first at its output, to yield the dual reflection matrix $R_{up}(z) = [R(\mathbf{u}_{out}, \rho_{in}, z)]$. This projection can be performed by performing the following matrix product:

$$ R_{up}(z) = G_0^\top(z, f_0) \times R_{pp}(z) \qquad (14) $$

where the symbol $\top$ stands for matrix transpose. An angular de-scan of the input focusing points as sketched in Fig. 4b can be performed by computing the Hadamard product between $R_{up}$ and its ideal counterpart $G_0(z)$:

$$ D_{up}(z) = G_0^\dagger(z, f_0) \circ R_{up}(z) \qquad (15) $$

Each column of the resulting distortion matrix, $D_{up}(z) = [D(\mathbf{u}_{out}, \rho_{in}, z)]$, maps the phase-distortions with respect to the ideal wave-front that would be obtained for a point-like source at $(\rho_{in}, z)$.

An estimator $W_u$ of the aberration transmittance is then extracted through an iterative phase reversal (IPR) process applied to the correlation matrix $C_{uu} = D_{up} \times D_{up}^\dagger$ (see below). The phase conjugate of the estimator $W_u$ is then used as a focusing law to compensate (partially) for wave distortions. An updated focused reflection matrix is obtained through the following relation:

$$ R_{pp}(z) = G_0^\dagger(z, f_0) \times [G_0(z, f_0) \circ W^*(z) \circ D_{up}(z)] \qquad (16) $$

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 [Fig. 4e and f]. The fine compensation of wave distortions is highlighted by the RSPF (Eq. (12)) deduced from the updated focused R-matrix (Fig. 4g). As expected theoretically, compensation of aberrations in the geophone basis enables the recovery of a resolution only limited by the geophone aperture. As shown by Fig. 4e, this is nevertheless not sufficient to have a contrasted image of the volcano in depth.

## Compensation of diffraction

To go beyond and beat diffraction, the parabolic phase law exhibited by the RPSF (Fig. 4h) can be exploited (Supplementary Note 7). To that aim, the R-matrix shall be investigated between the focused basis and the k-space. The focused reflection matrix is thus projected in the plane wave basis, first at output, such that:

$$R_{kp}(z) = T_0 \times R_{pp}(z), \quad (17)$$

with $T_0$, the Fourier transform operator

$$T_0(k_{||}, \rho) = \exp(-ik_{||} \cdot \rho). \quad (18)$$

Then, a new distortion matrix can be built in the plane wave basis by comparing each reflected wave field, $R(k_{out}, \rho_{in}, z)$, in the k space with its reference counterpart, $T_0(k_{out}, \rho_{in})$, that would be obtained for a point-like guide star at $(\rho_{in}, z)$:

$$D_{kp}(z) = R_{kp}(z) \circ T_0^*(z) \quad (19)$$

From this matrix $D_{kp}$, a diffraction transmittance $W_k$ (Fig. 5c) can be extracted by applying the IPR process to the correlation matrix $C_{kk} = D_{kp} \times D_{kp}^\dagger$ (see below).

The phase conjugate of $W_k$ is then used as a focusing law to compensate for wave diffraction (Supplementary Note 9). An updated focused reflection matrix is obtained through the following relation:

$$R_{pp}(z) = T_0^\dagger \times [T_0 \circ W_k^*(z) \circ D_{kp}(z)] \quad (20)$$

The whole process is then repeated at input to compensate for the diffraction phenomena undergone by the down-going wave fields. At the end of the process, a novel confocal image is obtained by considering the diagonal elements of the updated focused reflection matrix [Fig. 5d and e]. The beating of diffraction is highlighted by the RSPF (Eq. (12)) deduced from the updated focused R-matrix (Fig. 5f). The resolution reaches the ultimate diffraction limit (~$\lambda/2$) and is no longer limited by the geophone network aperture (white circle in Fig. 5f).

## Iterative phase reversal

The IPR algorithm is a computational process that provides an estimator of the aberration and/or diffraction transmittance at each depth $z$ extracted from the correlation matrices computed at the Earth surface ($C_{uu}$) and in the plane wave basis ($C_{kk}$), respectively. Mathematically, the algorithm is based on the following recursive relation:

$$W_x^{(n+1)} = \exp(i\arg\{C_{xx} \times W_x^{(n)}\}) \quad (21)$$

with $x = u$ or $k$, the coordinate vector in the correction basis. $W_x^{(0)} = [1 \dots 1]^T$ is chosen arbitrarily as a unit wave-front. The resulting wave-front at the end of the IPR process, $W_x = \lim_{n \to \infty} W_x^{(n)}$, provides an estimator of the aberration transmittance in the Earth basis (Supplementary Note 8) and of the diffraction transmittance in the plane wave basis (Supplementary Note 9).

## Data 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 are available at Zenodo⁴¹ (https://zenodo.org/record/10066910).

## Code availability

The codes used to post-process the noise correlation data are available at Zenodo⁴¹ (https://zenodo.org/record/10066910).

Received: 21 December 2023; Accepted: 28 August 2024; Published online: 16 September 2024

## References

1. Mosk, A. P., Lagendijk, A., Lerosey, G. & Fink, M. Controlling waves in space and time for imaging and focusing in complex media. *Nat. Photonics* **6**, 283–292 (2012).

2. Moustakas, A. L., Baranger, H. U., Balents, L., Sengupta, A. M. & Simon, S. H. Communication through a diffusive medium: coherence and capacity. *Science* **287**, 287–290 (2000).

3. Chau, G., Jakovljevic, M., Lavarello, R. & Dahl, J. A locally adaptive phase aberration correction (LAPAC) method for synthetic aperture sequences. *Ultrason. Imaging* **41**, 3–16 (2019).

4. Lambert, W., Cobus, L. A., Frappart, T., Fink, M. & Aubry, A. Distortion matrix approach for ultrasound imaging of random scattering media. *Proc. Nat. Acad. Sci. USA* **117**, 14645–14656 (2020).

5. Kang, S. et al. High-resolution adaptive optical imaging within thick scattering media using closed-loop accumulation of single scattering. *Nat. Commun.* **8**, 2157 (2017).

6. Badon, A. et al. Distortion matrix concept for deep optical imaging in scattering media. *Sci. Adv.* **6**, eaay7170 (2020).

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

8. Weaver, R. L. & Lobkis, O. I. Ultrasonics without a source: thermal fluctuation correlations at MHz frequencies. *Phys. Rev. Lett.* **87**, 134301 (2001).

9. Campillo, M. & Paul, A. Long-range correlations in the diffuse seismic coda. *Science* **299**, 547–549 (2003).

10. Blondel, T., Chaput, J., Derode, A., Campillo, M. & Aubry, A. Matrix approach of seismic imaging: application to the Erebus volcano, Antarctica. *J. Geophys. Res. Solid Earth* **123**, 10936–10950 (2018).

11. Touma, R., Blondel, T., Derode, A., Campillo, M. & Aubry, A. A distortion matrix framework for high-resolution passive seismic 3-D imaging: Application to the San Jacinto fault zone, California. *Geophys. J. Int.* **226**, 780–794 (2021).

12. Touma, R., Le Ber, A., Campillo, M. & Aubry, A. Imaging the crustal and upper mantle structure of the North Anatolian Fault: A transmission matrix framework for local adaptive focusing. *J. Geophys. Res. Solid Earth* **128**, e2023JB026704 (2023).

13. Shapiro, N. M., Campillo, M., Stehly, L. & Ritzwoller, M. H. High-resolution surface-wave tomography from ambient seismic noise. *Science* **307**, 1615–1618 (2005).

14. Roux, P., Sabra, K. G., Gerstoft, P., Kuperman, W. A. & Fehler, M. C. P-waves from cross-correlation of seismic noise. *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. Burtin, A. et al. Dense seismic monitoring of La Soufrière de Guadeloupe hydrothermal system, https://www.fdsn.org/networks/detail/ZK_2017 (2018).

17. Schuster, G. T. et al. Chap. 5. Redatuming. In *Seismic Interferometry*, 331–448 (Society of Exploration Geophysicists, 2008).
18. Sentenac, A. & Mertz, J. Unified description of three-dimensional optical diffraction microscopy: From transmission microscopy to optical coherence tomography: tutorial. *J. Opt. Soc. Am. A* **35**, 748 (2018).

19. Najar, U. et al. Harnessing forward multiple scattering for optical imaging deep inside an opaque medium. *Nat. Commun.* **15**, 7349 (2024).

20. Bureau, F. et al. Three-dimensional ultrasound matrix imaging. *Nat. Commun.* **14**, 6793 (2023).

21. Moretti, R. et al. The 2018 unrest phase at La Soufrière de Guadeloupe (French West Indies) andesitic volcano: Scrutiny of a failed but prodromal phreatic eruption. *J. Volcanol. Geotherm. Res.* **393**, 106769 (2020).

22. Metcalfe, A. et al. Magmatic processes at La Soufrière de Guadeloupe: Insights from crystal studies and diffusion timescales for eruption onset. *Front. Earth Sci.* **9**, 617294 (2021).

23. Metcalfe, A., Moune, S., Komorowski, J.-C. & Moretti, R. Bottom-up vs top-down drivers of eruption style: petro-geochemical constraints from the Holocene explosive activity at La Soufrière de Guadeloupe. *J. Volcanol. Geotherm. Res.* **424**, 107488 (2022).

24. Dorel, J., Eschenbrenner, S. & Feuillard, M. Coupes sismiques des structures superficielles dans les petites Antilles: I: Guadeloupe. *Pure Appl. Geophys.* **117**, 1050–1069 (1979).

25. Berkhout, A. J. & Wapenaar, C. P. A unified approach to acoustical reflection imaging. II: the inverse problem. *J. Acoust. Soc. Am.* **93**, 2017–2023 (1993).

26. Lambert, W., Cobus, L. A., Couade, M., Fink, M. & Aubry, A. Reflection matrix approach for quantitative imaging of scattering media. *Phys. Rev. X* **10**, 021048 (2020).

27. Pichavant, M., Poussineau, S., Lesne, P., Solaro, C. & Bourdier, J. L. Experimental parametrization of magma mixing: application to the ad 1530 eruption of La Soufrière, Guadeloupe (Lesser Antilles). *J. Petrol.* **59**, 257–282 (2018).

28. Osnabrugge, G., Horstmeyer, R., Papadopoulos, I. N., Judkewitz, B. & Vellekoop, I. M. Generalized optical memory effect. *Optica* **4**, 886–892 (2017).

29. Lambert, W., Cobus, L. A., Robin, J., Fink, M. & Aubry, A. Ultrasound matrix imaging—Part II: The distortion matrix for aberration correction over multiple isoplanatic patches. *IEEE Trans. Med. Imag.* **41**, 3921–3938 (2022).

30. Cashman, K. V., Sparks, R. S. J. & Blundy, J. D. Vertically extensive and unstable magmatic systems: A unified view of igneous processes. *Science* **355**, eaag3055 (2017).

31. Cassidy, M., Manga, M., Cashman, K. & Bachmann, O. Controls on explosive-effusive volcanic eruption styles. *Nat. Commun.* **9**, 2839 (2018).

32. Edmonds, M., Cashman, K. V., Holness, M. & Jackson, M. Architecture and dynamics of magma reservoirs. *Philos. Trans. Royal Soc. A* **377**, 20180298 (2019).

33. Combier, V. et al. Three-dimensional geometry of axial magma chamber roof and faults at Lucky Strike volcano on the mid-Atlantic ridge. *J. Geophys. Res. Solid Earth* **120**, 5379–5400 (2015).

34. Edmonds, M., Kohn, S., Hauri, E., Humphreys, M. & Cassidy, M. Extensive, water-rich magma reservoir beneath southern Montserrat. *Lithos* **252-253**, 216–233 (2016).

35. Koulakov, I. et al. Causes of volcanic unrest at Mt. Spurr in 2004–2005 inferred from repeated tomography. *Sci. Rep.* **8**, 17482 (2018).

36. Winslow, H. et al. Insights for crystal mush storage utilizing mafic enclaves from the 2011–12 Cordón Caulle eruption. *Sci. Rep.* **12**, 9734 (2022).

37. Papale, P. & Marzocchi, W. Volcanic threats to global society. *Science* **363**, 1275–1276 (2019).

38. Bensen, G. et al. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements. *Geophys. J. Int.* **169**, 1239–1260 (2007).

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

40. Aki, K. Analysis of the seismic coda of local earthquakes as scattered waves. *J. Geophys. Res.* **74**, 615–631 (1969).

41. Giraudat, E. et al. Unveiling the deep plumbing system of a volcano by a reflection matrix analysis of seismic noise [data, code]. Zenodo. https://zenodo.org/record/10066910 (2023).

## Acknowledgements

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. 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 de la transition écologique et de la cohésion des territoires (MTECT) 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.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-024-01659-2.

**Correspondence** and requests for materials should be addressed to Alexandre Aubry.

**Peer review information** *Communications Earth & Environment* thanks the anonymous reviewers for their contribution to the peer review of this work. Primary Handling Editor: Carolina Ortiz Guerrero. A peer review file is available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
**Open Access** This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

© The Author(s) 2024