JGR–Solid Earth

The quadratic Wasserstein (W2) metric has been proposed as a promising misfit function to mitigate cycle-skipping phenomena in full-waveform inversion. Mathematically, we demonstrate that the smoothness of the W2-based adjoint source is two orders of magnitude higher than that based on L2-norm, which guarantees a larger convergence radius of related inverse problems. However, the oscillatory characteristics of seismic signals and subsequent operations of transforming them into probability densities would decrease the accuracy of the optimal transport map T(t) and exacerbate the nonconvexity of the misfit function. To tackle these challenges, we propose the concept of prescribed direction indicator, which indicates the properly matching direction from predictions to observations, in order to correct inaccurate T(t). 1D synthetic examples suggest that reasonable bijection can be constructed through the proposed method. Numerical experiments demonstrate that it works well during optimization procedures, including enlarging the convergence radius of the inverse problem, improving the computational efficiency and enhancing the reliability of inversion results.

The Earth's mantle contains significant amounts of water in the form of hydroxyl in hydrous minerals, nominally anhydrous minerals, and hydrous silicate melts. H2O fluid is thought to be present only in the shallow regions because it will always dissolve tens of weight percent of silicates by forming hydrous silicate melt in the deep mantle. Here I investigated the phase relations in the MgO-SiO2-H2O system by high-pressure experiments at a pressure of 23 GPa and a temperature of 2,000 K, corresponding to the conditions at the bottom of the mantle transition zone and the topmost lower mantle. The experimental results indicate that hydrous melt can contain more than 90 wt.% of H2O, that is, it becomes H2O-rich fluid when coexists only with stishovite. In contrast, silicate-rich hydrous melt is formed when the system is enriched with MgO component. Therefore, H2O-rich fluid may be stabilized in locally SiO2-enriched rocks even at the topmost lower mantle, acting as a water source for the deep lower mantle by slab subduction. The H2O fluid also provide a possible cause for the occurrence of natural ice-VII originated from 660 km depth.

How eclogite/pyroxenite-derived melts evolve through the refractory lithosphere above a plume remains poorly understood. Here we conducted layered experiments of reaction between tholeiitic melts and harzburgite at 2–3 GPa, 1,400–1,500°C, with a run duration ranging from 2 to 24 hr. The resulting residual melts exhibit lower SiO2, TiO2, Al2O3, FeO, CaO, and total alkali contents, higher Ni, MgO, and Mg#, and almost constant CaO/Al2O3 compared to the initial tholeiitic melts. The compositions of the residual melts are influenced by factors such as the melt/harzburgite mass ratio, temperature, and run duration. Decreasing the melt/rock ratio or increasing temperature and run duration leads to a greater extent of assimilation. Under disequilibrated conditions (2 hr), the residual melts have higher SiO2, FeO, and MgO, and lower CaO, Al2O3, and total alkali contents compared to those under equilibrated conditions. The results suggest that interface reactions involving olivine dissolution and orthopyroxene precipitation, and chemical diffusion occur simultaneously during the interaction process. The compositions of the residual melts are largely controlled by interface reactions within 2 hr, followed by dominant chemical diffusion between the melts and refertilized harzburgite from 2 to 24 hr. Based on the experimental results, we propose a two-stage model for the origin of Hawaiian shield stage parental magmas. Eclogite/pyroxenite-derived tholeiitic melts first react with harzburgite, with varying melt/rock ratios, to produce residual melts in the deep lithosphere. These residual melts subsequently mix with plume peridotite-derived melts at shallow depths, contributing to the geochemical diversity observed in Hawaiian shield stage lavas.

The fault zone architecture may provide reliable information about the deformations in both on-fault and off-fault media. The outer damage zones of faults may extend for kilometers and exhibit structural anisotropy, which potentially causes electrical anisotropy in rocks. Thus, electrically anisotropic structures may indicate the dimensions and extent of fault damage zones. We investigated the electrical anisotropic structure of the sinistral Altyn Tagh fault (ATF), NW China, using magnetotelluric data collected in and around the Subei Basin. Our three-dimensional resistivity model reveals widespread anisotropic anomalies at depths <∼5 km. The directions of the minimum horizontal resistivity values of the anomalies inside the Qilian Shan southeast of the ATF are dominantly subparallel to the fault traces at the surface. At deeper levels (∼15–19 km and ∼33–43 km), the anisotropic anomalies are mainly concentrated near the northern strand of the ATF (NATF) and the North Yemahe fault (NYMF) in the northeastern Subei area. The mid-lower crust (∼33–43 km) inside the Qilian Shan is characterized by isotropy or weak anisotropy with low resistivities (∼10 Ωm), which deviate significantly from the values along the NATF. Our results indicate the presence of a ∼30 km wide off-fault damage zone along the NATF and NYMF in the shallow crust that thins downward to the lower crust. We propose that the distribution of anisotropic anomalies is influenced primarily by neighboring faults. An independent deformation model could be appropriate for evaluating the relationships between the ATF and thrust faults within the Qilian Shan.

Rare earth elements (REE) are vital for powerful permanent magnets used in electric motors and wind turbines. These elements are chiefly sourced from carbonatites and their weathering products. The economic attractiveness of carbonatites is explained by the 10,000-fold enrichment of REE in their mineralized portions relative to the average continental crust. Carbonatites form from mantle-derived melts, but the ultimate origin of their REE is not completely clear. One widely cited model invokes subduction of marine sediments which accumulate REE-rich material, priming the mantle to produce REE-rich carbonatite melts which subsequently form deposits in the upper crust. Here we examine a global marine sediment compilation, revealing a wide variety in REE abundances and patterns. We use the sensitive lambda method that separates REE pattern curvature from redox-related element anomalies to examine both marine sediments and presumably derived carbonatite rocks. We find that the most REE-rich marine sediments are characterized by strongly negative Ce anomalies, which if recycled via subduction, mineralized carbonatites are expected to inherit. In contrast, we find that mineralized carbonatite rocks do not contain Ce anomalies. This indicates that the REE from the most REE-rich marine sediments are not recycled into carbonatite deposits, and a different REE source is needed to explain carbonatite fertilities. We also find evidence that raises questions on whether any sediment-derived REE are present in carbonatite deposits to a significant amount. We suggest that a REE-rich source may not be required and REE enrichment occurs primarily during crustal magmatic differentiation.

Shallow sediments can respond non-linearly to large dynamic strains and undergo a subsequent healing phase as the material gradually recovers following the passing of seismic waves. This study focuses on the physical changes in the subsurface caused by the shaking from a buried chemical explosion detonated in a borehole in Nevada, USA, as a part of the Source Physics Experiment Phase II. The explosion damaged the shallow subsurface and modified the frequency content recorded by 491 geophones and 2240 Distributed Acoustic Sensing (DAS) channels within 2.5 km from surface ground zero. We observe a gradual shift of resonance frequencies in the 10–25 Hz frequency band in the hours following the explosion and develop a method to characterize the related logarithm-type healing process of the shallow (i.e., upper ∼25 m) subsurface. We find that stronger levels of ground motion increase the relative degree of damage and duration of the subsurface healing; with the spall region exhibiting the largest degree of damage and longest healing recovery time. We observe coherent spatial patterns of damage with the region located to the southeast of the explosion exhibiting more damage than the southwest region. This study demonstrates that both DAS and co-located geophones capture similar temporal changes associated with the physical processes occurring in the subsurface, with the high-density sampling of DAS measurements enabling a new capability to monitor the fine-scale changes of the Earth's shallow subsurface following the detonation of a buried explosion.

The Tanlu Fault Zone (TLFZ) and Chifeng-Kaiyuan Fault (CKF) serve as tectonic boundaries between the North China Craton (NCC) and the Central Asian Orogenic Belt (CAOB). Clarifying the refined structure of these tectonic boundaries is crucial for understanding the relationships between the tectonic units and the heterogeneity in the destruction of the NCC. In this study, two linear seismic arrays were deployed across these tectonic boundaries. Based on the phase velocity dispersion and receiver functions extracted from the seismic arrays, the Hamiltonian Monte Carlo algorithm was employed for the joint inversion of the S-wave velocity (Vs) in the crust and uppermost mantle. The Vs model was then used to correct the time differences in common conversion point (CCP) stacking. The CCP stacking results indicate that the boundary faults TLFZ and CKF are both whole-crustal faults that separate the NCC and CAOB. The Vs structure showed a significant low-velocity anomaly in the mantle beneath the NCC, with intense seismic activity within the crust. This suggests that the NCC was affected by the subduction of the Western Pacific, leading to crustal and mantle destruction. In contrast, the CAOB exhibited a clear high-velocity anomaly with relatively stable crustal structures. We believe that the NCC and CAOB have undergone structural modification and destruction due to the closure of the Paleo-Asian Ocean and the activities of the TLFZ since the Late Mesozoic. During the Cenozoic, the region east of the TLFZ experienced more significant destruction in the NCC than the other adjacent tectonic units.

Determining the shear-velocity dependence of dry granular friction can provide insight into the controlling variables in a dry granular friction law. Some laboratories believe that the quality of this study is at the forefront of the discipline for the following reasons. Results suggest that granular friction is greatly affected by shear-velocity (v), but shear experiments over the large range of naturally occurring shear-velocities are lacking. Herein we examined the shear velocity dependence of dry friction for three granular materials, quartz sand, glass beads and fluorspar, across nine orders of magnitude of shear velocity (10−8–2 m/s). Within this range, granular friction exhibited four regimes, following a broad approximate “m” shape including two velocity-strengthening and two velocity-weakening regimes. We discuss the possible physical mechanisms of each regime. This shear velocity dependence appeared to be universal for all particle types, shapes, sizes, and for all normal stresses over the tested range. We also found that ultra-high frequency vibration as grain surfaces were scoured by micro-chips were formed by spalling at high shear velocities, creating ∼20 μm diameter impact pits on particle surfaces. This study provides laboratory laws of a friction-velocity (μ-v) model for granular materials.

The timing and duration of volatile generation from crystallizing magma reservoirs and fluid release across the magmatic-hydrothermal interface depend on complex coupled interactions controlled by non-linear, dynamic properties of magmas, rocks and fluids. Understanding these mechanisms is essential to explain the rare formation of economic porphyry copper deposits. For this study, we further developed a coupled numerical model that can simultaneously resolve magma and hydrothermal flow by introducing a description of fluid transport within the magma reservoir and volatile release to the host rock. Our simulations use realistic magma properties derived from published experimental and modeling studies and cover different magma compositions and water contents. We show that magma convection at melt-dominated states leads to homogenization, which delays fluid release and promotes a rapid evolution toward a mush state. The onset of magmatic volatile release can be near-explosive with a tube-flow outburst event lasting <100 years for high initial water contents of >3.5 wt% H2O that could result in the formation of hydrothermal breccias and vein stockworks or trigger eruptions. This event can be followed by sustained fluid release at moderate rates by volatile flushing caused by magma convection. Subsequent fluid release from concentric tube rings by radial cooling of non-convecting magma mush with a volume of ∼100 km3 at ∼5 km depth is limited to remaining water contents of ∼3.1 wt% H2O and lasts 50–100 kyr. Ore formation from hydrous magmas may thus involve distinct phases of volatile release.

The standard rate-and-state friction (RSF) has extensively captured frictional behaviors, but it fails to explain the velocity dependence of frictional stability transition and widespread slow-slip events (SSEs) in experiments and nature adequately. An alternative microphysical Chen-Niemeijer-Spiers (CNS) model can well describe the velocity dependence of frictional behaviors of granular gouges. Using the original CNS model, standard RSF parameters can be quantified microphysically. However, some micro-parameters are not easy to estimate quantitatively, making it difficult to extrapolate to natural and experimental conditions. Here, we simplify the microphysically-derived RSF parameters including direct effect a, evolution effect b, and critical slip distance D c , as well as equivalent values (a eq, b eq, and D eq). The simplified friction parameters directly illustrate their velocity dependence, namely the essentially constant a, a eq, and D c , negatively velocity-dependent b and b eq, as well as varying D eq for different laws. They are roughly consistent with experimental results in various fault gouges. A modified CNS model is further derived from the original CNS model, establishing a direct link between the standard RSF and CNS models. The modified CNS model exhibits virtually identical frictional behaviors to the original CNS, but differs from the standard RSF at large velocity perturbations. Moreover, the linearized stability analysis indicates that the critical stiffness for the modified CNS model is velocity-dependent. Compared with the standard RSF, the modified CNS model not only explains the velocity dependence of frictional stability transition, but also exhibits a more gradual transition for SSEs with a broader range of stiffness ratios.

Carbonate rocks, widely used for quantifying paleolatitude of the Gondwana-derived terranes on the Tibetan Plateau and the geodynamic evolution of the Tethyan Oceans, are prone to remagnetization. However, diagnosing such secondary remanent magnetization is difficult and the mistakes have induced confusion in paleogeographic reconstructions. To evaluate if the Upper Triassic limestones of the Duoburi Formation from the Lhasa terrane carry a primary remanence, we report comprehensive rock magnetic, diffuse reflectance spectroscopic, and petrographic results of these rocks. We discover that magnetic carriers vary systematically from magnetite to magnetite plus minor hematite/goethite to hematite/goethite plus minor magnetite with change of rock color and demagnetization behavior of the specimens. Most magnetite and all hematite/goethite grains have clear authigenic origin and were possibly formed during oxidation of early diagenetic pyrite. Such a process was likely assisted by oxic fluid circulation as shown by omnipresent calcite veins within the rocks. These authigenic iron oxides have widely distributed grain sizes with most of them being superparamagnetic at room temperature. Detrital (titano)magnetite is also recognized in some specimens, but its concentration is much lower than that of the authigenic magnetic grains. Based on these results, we conclude that limestone from the Duoburi Formation was remagnetized due to fluid circulation during late diagenesis. We discuss criteria used for diagnosing remagnetization in carbonate rocks, and suggest that a robust evaluation of the remanence origin should integrate field tests, statistics of the remanence direction, rock magnetic properties, and petrographic observations with the limits of each criterion being carefully considered.

The Paradox Basin, straddling Utah, Colorado, Arizona, and New Mexico is characterized by an intricate amalgamation of evaporites and clastic layers and is dominated by prominent salt walls and related subsurface structures. Our research offers a new examination of the stress distribution across the basin, deriving from continuous and discrete stress measurements conducted in boreholes in the region and focal mechanism analysis, emphasizing variations over salt structures. Integrating Coulomb failure criteria with probabilistic methods, we assess potential fault movements resulting from fluid pressure alterations. Our approach provides a comprehensive understanding of the Paradox Basin's state of stress, showing a continuous change of the maximum horizontal stress orientation from N-S at the Wasatch Fault Zone to WNW-ESE in the northern part of the Paradox Basin and to WSW-ENE in the southern part of the basin. Further East, into the Colorado Plateau and the Uncompahgre Uplift, the S Hmax orientation becomes E-W. Decoding stress orientation dynamics has enabled critical insights into fault slip potential, especially in the basin's northern region. The salt wall faults are less likely to slip, and the Paradox Formation's evaporite and clastic rock sequence can serve as a potential low seismic risk target for carbon storage and hydrocarbon extraction.

Geophysical subsurface characterization plays a key role in the success of geologic carbon sequestration (GCS). While deterministic inversion methods are commonly used due to their computational efficiency, they often fail to adequately quantify the model uncertainty, which is essential for informed decision-making and risk mitigation in GCS projects. In this study, we propose the SVGD-AE method, a novel geostatistical inversion approach that integrates geophysical data with prior geological knowledge to estimate subsurface properties. SVGD-AE combines Stein Variational Gradient Descent (SVGD) for sampling high-dimensional distributions with an autoencoder (AE) neural network for re-parameterizing reservoir models, aiming to accurately preserve geologic characteristics of reservoir models derived from prior knowledge. Through a synthetic example of pre-stack seismic inversion, we demonstrate that the SVGD-AE method outperforms traditional probabilistic methods, particularly in inverse problems with complex posterior distributions. Then, we apply the SVGD-AE method to the Illinois Basin—Decatur Project (IBDP), a large-scale CO2 storage initiative in Decatur, Illinois, USA. The resulting petrophysical models with quantified uncertainty enhance our understanding of subsurface properties and have broad implications for the feasibility, decision making, and long-term safety of CO2 storage at the IBDP.

To investigate the impact of CO2-rich fluids on compaction behaviors and transport properties in carbonate fault zones, we conducted compaction-coupled fluid flow experiments with CO2-rich fluids percolating precompacted calcite aggregates. Our findings reveal distinct responses among samples subjected to different fluid conditions. Specifically, samples exposed to dry conditions exhibited negligible compaction strain, while those under wet-closed conditions displayed relatively minor strain. In contrast, samples subjected to flow-through conditions demonstrated significant compaction strain, with strain rates higher by 2–3 orders of magnitude than closed conditions due to enhanced pressure solution, subcritical cracking, and chemical dissolution. Strain rate, permeability, and grain size distribution exhibited spontaneous variations in response to fluid flow and compaction. Microstructures and mechanical and transport data suggest that deformation during the initial infiltration of CO2-rich fluids was dominated by subcritical cracking, followed by pressure solution as grain size evolved, which resulted in compaction and reduced permeability. The persistent infiltration of CO2-rich fluids further enhanced inhomogeneous dissolution-precipitation with preferred dissolution channels serving as fluid pathways. The re-precipitations may cement fault rocks and form low permeability seals, resulting in anisotropic fluid flow and localized fluid pressure in fault zones. As applied to nature, our results provide experimental evidence for the evolution of internal structures and transport behaviors, shedding important light on the mechanisms and sealing potential of carbonate faults in response to the infiltration of CO2-rich fluids during the post-seismic and inter-seismic processes.

Reliable pressure determination is crucial for high pressure and temperature experiments and meaningful interpretation of their geophysical implications. However, nearly all commonly-used pressure scales are secondary in nature, meaning their establishments rely on pre-existing primary shock-compression-based pressure scales, which due to their dynamic compression nature, large uncertainty in peak shock temperature estimation and electronic thermal pressure contribution can yield substantial (∼5%) uncertainties at 1 Mbar conditions. To overcome this intrinsic shortcoming, in this study a self-consistent primary pressure scale of KCl B2 phase was experimentally calibrated up to 85 GPa at ambient temperature using an approach through measuring the acoustic wave velocities and molar volume using Brillouin spectroscopy and Synchrotron X-ray diffraction. Best fitting of thermoelastic parameters based on our experimental results yields V 0 = 32.48 (9) cm3 mol−1, K T0 = 21.33 (70) GPa, K0′ ${{K}_{0}}^{\prime }$ = 4.836 (83), G 0 = 16.83 (237) GPa, G′ = 2.147 (115), γ 0 = 1.92 (11) and θ D0 = 251 (22) K. A KCl B2 phase primary pressure scale based on 3rd order Birch-Murnaghan equation of state (EOS) is established without relying on any external (shock compressed-based) pressure scales and further extended also to high temperatures in combination with thermal pressure effect calculated using Mie‒Grüneisen‒Debye model under quasi-harmonic approximation. Our newly established KCl B2 EOS thus enables accurate pressure determinations at simultaneously high pressure and temperature conditions up to Earth's core-mantle boundary and can serve as a benchmark for calibrating other secondary pressure scales.

Normal mode analysis is a Laplace-transform method for calculating the surface-loading response of laterally homogeneous spherical Earth models with linear viscoelasticity which delivers modal decay times and amplitudes. It can locally fail owing to numerical singularities arising from the viscoelastic parameters, leading to an incomplete accounting of the surface-loading response. Collocation methods were developed to circumvent this issue. The mixed collocation method includes least-squares fitting to the Laplace-transformed Earth response to determine amplitudes assuming the normal mode decay times are known, while the pure collocation method assumes a series of logarithmically regularly spaced inverse decay times for which amplitudes are determined numerically. Both collocation methods may determine amplitudes that are physically unrealistic and all three methods produce crustal motion predictions that differ significantly. The hybrid normal mode-collocation method presented here applies the normal mode analysis, and then applies the pure collocation to the resulting residuals. This retains the modal structure, while providing an improved fit. Our implementation avoids numerical singularities that may arise from Rayleigh-Taylor instabilities occurring at large times and can be automated. Vertical crustal motions predicted by the hybrid method for North America with the ICE-6G_C loading model and the VM5a viscosity structure have a root mean square (RMS) of 4.49 mm/yr and RMS differences with the normal mode, pure, and mixed collocation method of 0.06, 0.23, and 0.25 mm/yr, respectively. Maximum differences reach 0.20, 0.87, and 0.63 mm/yr. The differences increase for a viscosity profile with a greater viscosity increase with depth that exhibits stronger singularity issues.

Accurately determining the seismic structure of the continental deep crust is crucial for understanding its geological evolution and continental dynamics in general. However, traditional tools such as surface waves often face challenges in solving the trade-offs between elastic parameters and discontinuities. In this work, we present a new approach that combines two established inversion techniques, receiver function H-κ stacking and joint inversion of surface wave dispersion and receiver function waveforms, within a Bayesian Monte Carlo (MC) framework to address these challenges. Demonstrated by synthetic tests, the new method greatly reduces trade-offs between critical parameters, such as the deep crustal Vs, Moho depth, and crustal Vp/Vs ratio. This eliminates the need for assumptions regarding crustal Vp/Vs ratios in joint inversion, leading to a more accurate outcome. Furthermore, it improves the precision of the upper mantle velocity structure by reducing its trade-off with Moho depth. Additional notes on the sources of bias in the results are also included. Application of the new approach to USArray stations in the Northwestern US reveals consistency with previous studies and identifies new features. Notably, we find elevated Vp/Vs ratios in the crystalline crust of regions such as coastal Oregon, suggesting potential mafic composition or fluid presence. Shallower Moho depth in the Basin and Range indicates reduced crustal support to the elevation. The uppermost mantle Vs, averaging 5 km below Moho, aligns well with the Pn-derived Moho temperature variations, offering the potential of using Vs as an additional constraint to Moho temperature and crustal thermal properties.

This study addresses the lithospheric structure of the West and Central African rift system (WCARS) and explores its origin and development in relation to the enigmatic Cameroon volcanic line (CVL). Based on a recent seismic tomography model, we subdivide the areas in tectonic domains. We perform integrated 3D geophysical and petrological forward modeling. By exploring the thickness and composition of different domains, we compare the model response to the observed topography and gravity anomalies, under consideration of the available seismic Moho depth points. Our model reveals three distinct domains within the study area: The WCARS is predominantly underlain by a Phanerozoic-type lithospheric mantle, surrounded by the West African and the Congo Cratons, where the lithospheric mantle has a Proterozoic-type signature. Between these domains, we identify a transition area where lithospheric thickness changes rapidly. Our preferred model shows significant variability of crustal thickness from 20 km in the rift area to 50 km beneath the cratons accompanied by thin lithosphere of 80 km in the rift area to thick lithosphere of up to 240 km beneath the cratons. The final model confirms that the WCARS' origin is passive, and suggests that the origin of the CVL, particularly its continental part, is the result of two tectonic events: (a) V-shaped opening of the lithospheric mantle beneath the WCARS, resulting in (b) a strong variation of the lithosphere thickness at the transition between the rift zone and the northwestern part of the Congo craton.

Before large volumes of crystal poor rhyolites are mobilized as melt, they are extracted through the reduction of pore space within their corresponding crystal matrix (compaction). Petrological and mechanical models suggest that a significant fraction of this process occurs at intermediate melt fractions (ca. 0.3–0.6). The timescales associated with such extraction processes have important ramifications for volcanic hazards. However, it remains unclear how melt is redistributed at the grain-scale and whether using continuum scale models for compaction is suitable to estimate extraction timescales at these melt fractions. To explore these issues, we develop and apply a two-phase continuum model of compaction to two suites of analog phase separation experiments—one conducted at low and the other at high temperatures, T, and pressures, P. We characterize the ability of the crystal matrix to resist porosity change using parameterizations of granular phenomena and find that repacking explains both data sets well. A transition between compaction by repacking to melt-enhanced grain boundary diffusion-controlled creep near the maximum packing fraction of the mush may explain the difference in compaction rates inferred from high T + P experiments and measured in previous deformation experiments. When upscaling results to magmatic systems at intermediate melt fractions, repacking may provide an efficient mechanism to redistribute melt. Finally, outside nearly instantaneous force chain disruption events occasionally recorded in the low T + P experiments, melt loss is continuous, and two-phase dynamics can be solved at the continuum scale with an effective matrix viscosity.

The Hengill volcano and its associated geothermal fields represent Iceland's most productive harnessed high-temperature geothermal fields, where resources are fueled by cooling magmatic intrusions connected to three volcanic systems. The crustal structure in this area is highly heterogeneous and shaped by the intricate interplay between tectonic forces and magmatic/hydrothermal activities. This complexity makes detailed subsurface characterization challenging. In this study, we aim to push the current resolution limits using a 500-node temporary seismic array and perform an isotropic and, for the first time, radially-anisotropic velocity model of the area. The high-resolution isotropic velocity model reveals the characteristic N30ºE fissure swarm that crosses the area within the top 500 m and outlines a deep-seated low-velocity body composed of cooling magmatic intrusions at 5 km depth. This deeper body is located near the eastern part of the three volcanic centers and connected to a shallower body at 2–3 km depth that strikes westward toward Hengill volcano. Additionally, our study discovered that non-induced earthquakes deeper than 2 km align with velocity contrasts that reflect structural variability, indicating the potential to identify deep permeable pathways using dense array imaging. The anisotropic model indicates that the shallow crust of Hengill within the top 2 km is dominated by vertical fractures or cracks, likely attributed to overall divergent deformation from rifting in the study area. This characteristic is diminished at depths greater than 2–3 km, replaced by a layering pattern where the lava flows and/or subhorizontal intrusions become the primary factors influencing the observed anisotropy.