Precise monitoring of magma pressure and inflation within volcanic plumbing systems is crucial for accurate eruption forecasting and hazard assessment. However, imaging these deep systems is exceptionally challenging due to the complex geological heterogeneities that often disrupt conventional seismic migration techniques. Standard methods struggle with the irregular topography and complex velocity structures typically found in volcanic regions, often requiring large, dense networks of seismometers, which can be impractical and costly. This research addresses these limitations by applying a matrix imaging approach, borrowed from other fields like medical ultrasonics and optical microscopy. Matrix imaging offers several advantages over existing seismic methods, particularly full waveform inversion. It doesn't rely on a precise, often uncertain, wave velocity model and it is robust to the often-poor quality of seismic data frequently encountered in seismology. The study focuses on the La Soufrière volcano in Guadeloupe, a site known for its complex geological structure, to demonstrate the effectiveness of matrix imaging in providing high-resolution images of deep volcanic systems. The success of this method has important implications for improving our understanding of volcanic systems and enhancing the capabilities of volcano observatories to monitor dynamic changes and potentially predict eruptions more accurately.

Previous studies have utilized seismic noise interferometry to extract information about subsurface structures. Many studies focused primarily on surface wave properties from ambient noise correlation functions (NCFs). However, NCFs also contain information from body waves reflected by deeper structures and fluid reservoirs. Existing work on volcanic imaging, such as that performed on the Erebus volcano, showed reflectivity images of main structures, but lacked sufficient lateral resolution due to the absence of wave distortion compensation. Other studies, like the one on the San Jacinto Fault Zone, demonstrated high-resolution imaging with matrix imaging, but the cause of this super-resolution was not fully understood. This research builds upon these previous efforts, offering a refined methodology to improve both resolution and accuracy.

The study employed seismic noise data recorded by a sparse geophone network deployed at La Soufrière volcano. The network comprised 76 virtual geophones (a combination of temporary and permanent stations) spanning 1300 m laterally and 500 m vertically. Cross-correlation of ambient seismic noise provided the impulse response between each geophone pair, forming a time-dependent response matrix. To mitigate the influence of localized noise sources like fumaroles, only the anti-causal component of the noise correlation functions (NCFs) was used. This canonical reflection matrix was the foundation for subsequent analysis. The methodology involved several key steps: First, a homogeneous P-wave velocity model was adopted to perform a confocal redatuming process. This involved back-propagating the recorded echoes to retrieve reflectivity information at various depths. The quality of the initial confocal image was assessed using the off-diagonal elements of the focused reflection matrix, which revealed significant wave distortions caused by the mismatch between the simplified velocity model and the complex velocity structure of the volcano. To address these distortions, the researchers utilized the concept of the memory effect in wave physics. The memory effect states that a pattern of phase shifts imparted to a wavefront by an aberrating layer retains its shape, even when the wavefront's orientation is tilted. By leveraging this effect, the researchers developed an iterative phase reversal algorithm to estimate and compensate for the wave distortions. Two key improvements were implemented: first, an analysis of the wave distortions in the geophone basis allowed the compensation of wave distortions from the Earth's surface which improved the image quality at shallow depths. Second, a k-space analysis was used to overcome the diffraction limit imposed by the geophone array aperture. This involved projecting the reflection matrix into the k-space (plane-wave basis), where the diffraction pattern from each reflector appears as Fresnel rings. An iterative phase reversal process, applied in k-space, enabled the realignment of phase components, ultimately reducing the focal spot size to the diffraction limit of approximately half a wavelength (≈ 100 m).

The matrix imaging approach successfully revealed the internal structure of La Soufrière volcano with unprecedented detail down to 10 km depth and a resolution of about 100 meters. The high-resolution images revealed a complex three-dimensional structure. The upper part (depth < 5km) of the volcano shows a clear signature of a tortuous conduit. The deeper structure (5-8.5 km) exhibits a diffuse scattering pattern consistent with a vertical succession of several sub-horizontal, irregular, coalescing structures. These sub-horizontal structures, extending laterally over approximately 8 km, are interpreted as sub-horizontal lenses containing potentially eruptible magma. These are connected by narrower, sub-vertical structures. The reflectivity data shows variations in signal strength indicative of compositional differences within the magmatic system. Higher reflectivity zones are interpreted as areas with higher impedance contrasts, potentially due to the presence of gases and/or hydrothermal fluids. The overall three-dimensional image strongly supports the current model of transcrustal magmatic systems, showing a vertical arrangement of magma mush lenses, eruptible melt, and magmatic fluids. The depth range of the main magmatic plumbing system (5-8.5 km) aligns with independent petrological studies. This imaging technique achieves a remarkable improvement over previous studies, providing higher resolution and more detailed information about the complex geometry, layering, and connectivity of the volcano's internal structure. The observed high resolution, exceeding the diffraction limit imposed by the geophone array aperture, is largely attributed to the sparsity of reflectors at each depth and the compensation of the phase distortions and diffraction effects.

The results significantly advance our understanding of deep volcanic plumbing systems. The high-resolution images provide crucial insights into the geometry, dimensions, layered structure, and connectivity of the magma storage zones and the upper eruptive conduit. The impedance contrasts observed within the image offer potential for differentiating between zones of mush, eruptible melt, and their relative volume and position within the system. Further analysis could potentially yield estimates of pressure, temperature, volatile saturation, density contrast, and connectivity to the surface—parameters that are fundamental to understanding volcanic eruptions. The robustness of the matrix imaging method to the sparsity of the geophone array, inaccuracies in the velocity model, and imperfections in the noise correlation functions represents a major advance over computationally expensive methods like full waveform inversion. This suggests that the method could be effectively deployed in various seismic environments.

This study demonstrates that matrix imaging of seismic noise is a powerful tool for high-resolution monitoring of deep volcanic plumbing systems. The approach overcomes limitations of traditional seismic imaging techniques and provides unprecedented detail on the internal structure of volcanoes. The method’s robustness to data quality issues and its ability to exceed the diffraction limit open new avenues for volcano monitoring and eruption forecasting. Future research should focus on integrating matrix imaging with time-lapse monitoring to capture dynamic changes in volcanic systems and combine it with multi-parameter data from other monitoring networks.

While the study demonstrates remarkable results, it is important to acknowledge potential limitations. The use of a simplified homogeneous P-wave velocity model for redatuming might introduce some inaccuracies, although the method demonstrated robustness to model uncertainties. The interpretation of impedance contrasts relies on assumptions about the composition and properties of the subsurface materials, and further investigation may be needed to refine these interpretations. The study focused on a specific volcano, and further research is necessary to assess the generalizability of the findings to other volcanic systems with potentially different geological settings and seismic noise characteristics. Finally, the algorithm assumes single scattering; the effect of multiple scattering might influence the results, although it appears to be minimal in this particular geological setting.