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

The study addresses the challenge of high-resolution seismic imaging of complex volcanic and fault-zone structures using sparse arrays, where conventional migration and inversion methods struggle due to strong heterogeneities and limited array aperture. The research question is whether matrix imaging, adapted from ultrasonics and optics, can robustly reconstruct deep volcanic structures without requiring an accurate velocity model and despite sparse, noisy data. The context is the need to understand deep magma storage and plumbing systems for hazard assessment and eruption forecasting. The purpose is to develop and apply a passive seismic matrix imaging workflow to La Soufrière volcano using ambient noise, demonstrating improved resolution and robustness relative to traditional approaches. The significance lies in achieving detailed images of deep magmatic systems, potentially transforming volcano monitoring and forecasting.

Prior work in passive seismology using ambient noise has largely focused on extracting surface wave properties from noise correlation functions (NCFs). Matrix approaches in seismology have been demonstrated on Erebus volcano and the San Jacinto Fault Zone, but suffered from limited lateral resolution and lacked explicit aberration compensation. In wave physics, memory effects and distortion-matrix concepts have been used to correct aberrations and enhance resolution in ultrasound and optical imaging. Full waveform inversion offers detailed models but requires dense data and accurate starting models. This study builds on these advances, leveraging ambient noise interferometry to estimate the array response (reflection) matrix, and importing distortion-matrix and iterative phase reversal (IPR) methods to compensate for aberrations and diffraction, extending previous geophysical applications to deep volcanic imaging with higher resolution.

Data acquisition: A temporary nodal array of 65 three-component geophones (Zland 3C Gen2, 5 Hz natural frequency, 500 sps) and 6 permanent three-component broadband stations (flat response 1–50 Hz, 100 sps) operated by OVSG-IPGP were deployed on La Soufrière (Guadeloupe). Only vertical components were used. The temporary array operated from mid-November 2017 to mid-January 2018 in two sessions; 5 nodes were moved between sessions, yielding a virtual network of 76 sites. The array spans ~1300 m laterally and ~500 m vertically. Noise correlation processing: For each hourly vertical record: detrend and demean; deconvolve instrument response; band-pass 1–20 Hz; resample to 100 Hz; spectral whitening; 1-bit normalization. Compute NCFs by cross-correlating hourly records for lags −30 to +30 s; stack hourly to daily NCFs after quality control (correlation coefficient threshold 0.5); stack daily NCFs over two months. To suppress localized fumarolic sources, only the anti-causal component is retained to estimate impulse responses. The set of 2850 vertical impulse responses defines the canonical reflection matrix R(g_i, g_j, t). Velocity model: Only P waves (vertical component) are considered. A depth-dependent homogeneous P-wave velocity model c0(z) is defined from a 4-layer large-scale model, ranging from 2000 m/s near surface to 4300 m/s at 10 km depth. This rough model is sufficient within the matrix-imaging framework. Confocal redatuming: Transform R(t) to frequency domain R(f) over 10–20 Hz. Propagate to a focal plane at depth z at both emission and reception using free-space Green’s propagators G0(z, f); form focused reflection matrices R_pp(z, f) = G^† R(f) G. Coherently sum over frequency to obtain broadband R_pp(z), effectively time-gating singly scattered echoes near t ≈ 2z/c0(z), which improves axial resolution (Δz ≈ c0(z)/Δf). The confocal image I(ρ, z) is obtained from the diagonal elements R_pp(ρ, ρ, z). A depth gain compensation is applied based on an exponential fit of mean confocal intensity decay to estimate extinction length l_ext across depth ranges, and normalize images: I_x(ρ, z) = exp[β(z)/2] I(ρ, z). Focusing quality metric: The reflection point-spread function (RPSF) is computed along antidiagonals of R_pp(z) to assess lateral focusing and aberrations; ideally, energy concentrates within the aperture-limited focal spot of width δρ ≈ λ/(2 sin θ), with θ the viewing angle from array aperture. Strong off-diagonal spreading indicates aberrations and model mismatch. Aberration compensation via memory effect: Because near-surface heterogeneities distort wavefronts, the tilt memory effect is exploited. Project R_pp into a dual basis mixing the Earth surface (u) and focused basis (p) to build a distortion matrix D_up(z), representing deviations from ideal wavefronts for virtual sources at depth. Apply an iterative phase reversal (IPR) algorithm to the correlation matrix C_uu = D_up^T D_up to estimate an aberration transmittance W(u, z). Phase-conjugate W to adaptively refocus at input and output, iterating twice to obtain an updated R_pp(z) and improved confocal images. Diffraction compensation in k-space: To overcome the aperture-imposed diffraction limit at depth, project the updated focused reflection matrix to plane-wave (k) basis at output (and input), R_kp(z) = T_0 R_pp(z), where T_0 is the Fourier operator. Build D_kp(z) by comparing observed k-space wavefields to the ideal guide-star response. Due to Fresnel curvature, each reflector yields Fresnel ring patterns that populate the full diffraction disk (radius k0), not limited to array aperture. Apply IPR to C_kk = D_kp^* D_kp to estimate a diffraction transmittance W_k, and phase-conjugate to realign spatial frequency components, shrinking the focal spot to the ultimate diffraction limit δρ_0 ≈ λ/2 across depths. Iterate at input and output to update R_pp(z) and extract confocal images. Implementation notes and parameters: Depth slices from near surface to ~10.5 km were imaged. Bandwidth used: 10 Hz centered at 15 Hz. Lateral resolution achieved after k-space correction: ~λ/2 ≈ 100 m. Extinction lengths estimated from depth decay of mean confocal intensity: 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. The approach assumes sparse dominant reflectors per depth for effective super-resolved focusing.

- Passive matrix imaging of ambient seismic noise from a sparse array (76 sites; ~1.3 km lateral aperture; ~0.5 km vertical range) reveals La Soufrière’s internal structure down to ~10–10.5 km depth. - After aberration correction (tilt memory effect, distortion matrix, IPR) and k-space diffraction compensation (Fresnel ring realignment), the achieved lateral resolution reaches the ultimate diffraction limit δρ ≈ λ/2 ≈ 100 m over the full depth range, surpassing the aperture-limited resolution of conventional methods. - Shallow structure (< ~4–5 km): a twisted/helical sub-vertical conduit is resolved, consistent with an upper eruptive conduit connecting to the surface. - Deep structure (~5–8.5 km): a complex, multi-lens magma storage system is imaged, consisting of several sub-horizontal, irregular globular lenses extending laterally up to ~8 km, interconnected by narrow sub-vertical diffuse structures, consistent with a transcrustal mush-based system. - Reflectivity variations with depth: enhanced reflectivity above ~3.5–5 km (outer carapace) likely due to gas/liquid/supercritical fluids in porous host rock; lower reflectivity within 5–8.5 km consistent with extended magma volumes; increased reflectivity at ~8.5–10 km attributed to strong impedance contrast at eruptible melt/host rock interface. - Depth range of main magma storage constrained to ~5–8.5 km, matching independent petrological constraints for the 1530 CE eruption (top 5.6–7 km; base ≤ 8.5 km). - Robust imaging despite rough velocity model (c0 increasing 2000 to 4300 m/s) and imperfect NCF convergence; RPSF analysis confirms drastic reduction of aberrations and background after corrections. - Quantitative attenuation characterization via extinction length estimates: 1665 m (0–0.5 km), 310 m (0.5–1.9 km), 765 m (1.9–3.8 km), 3070 m (3.8–10.5 km). - Demonstrates that medium sparsity and wavefront curvature enable super-resolved focusing in seismology, explaining prior high-resolution results (e.g., San Jacinto) without invoking multiple scattering or channeling.

The findings demonstrate that passive seismic matrix imaging can overcome the twin challenges of aberrations from near-surface heterogeneities and aperture-limited diffraction that traditionally hamper deep imaging with sparse arrays. By leveraging memory effects, distortion matrices, and iterative phase reversal, the method corrects wavefront distortions without requiring an accurate velocity model. The k-space treatment harnesses Fresnel phase curvature to realign spatial frequencies beyond the array’s nominal aperture limit, achieving λ/2 resolution. These advances directly address the research goal: to image, at high resolution and depth, the magmatic plumbing of a volcano from ambient noise. The recovered architecture of La Soufrière’s transcrustal magmatic system matches independent petrological and conceptual models, validating the method and providing detailed constraints on geometry, connectivity, and depth extent of storage regions. This level of detail has important implications for hazard assessment: the impedance contrasts and structural connectivity observed can, with further analysis, inform estimates of pressure, temperature, volatile saturation, and melt/mush distribution, all key to eruption forecasting. The method’s robustness to sparse arrays and uncertain models offers practical advantages over more computationally intensive approaches like full waveform inversion, particularly in rugged volcanic terrains where dense deployments are impractical.

This work introduces and validates a passive seismic matrix imaging framework that uses ambient noise to build a reflection matrix, compensates aberrations via memory-effect-based distortion matrices and IPR, and beats aperture-limited diffraction through k-space focusing. Applied to La Soufrière, it resolves a twisted shallow conduit and a deep, multi-lens magma storage system from ~5 to ~8.5 km, with λ/2 (~100 m) lateral resolution down to ~10 km depth, in agreement with independent constraints. The approach offers a robust, high-resolution alternative for imaging complex volcanic systems with sparse arrays and imperfect velocity models. Future directions include time-lapse (4D) matrix imaging through repeated surveys to monitor temporal changes, integration with multi-parameter observatory data for joint interpretation, and application to other volcanic and tectonic settings where sparsity and wavefront curvature can be exploited for super-resolved monitoring.

- Performance relies on relative sparsity of dominant reflectors at each depth; if too many reflectors coexist, contrast and focusing may degrade. - Residual aberrations remain at larger depths even after correction, as indicated by secondary lobes in RPSFs prior to k-space compensation. - Use of a simplified, depth-dependent homogeneous P-wave velocity model may introduce focusing errors; S-wave and anisotropy effects are not modeled. - Only vertical components and P-wave responses are used; multi-component information could improve discrimination and robustness. - NCFs may not perfectly converge to Green’s functions; spurious arrivals and localized noise sources (e.g., fumaroles) can affect data quality despite mitigation (anti-causal component selection). - The method requires sufficiently broad bandwidth and adequate signal-to-noise; attenuation with depth necessitates compensation and may limit imaging in highly absorptive media. - Absolute physical property inference (e.g., melt fraction, temperature) from impedance contrasts requires additional modeling and multi-parameter constraints.