JGR–Solid Earth

This paper aims to update our understanding of the carbon cycle in the Himalayas, the most intense collisional orogeny globally, by providing new insight into its impact on Cenozoic climate cooling through the use of isotopic variations in both organic and inorganic carbon and an isotopic mass balance model. Our results from 20 selected hot springs show that the relative contributions of dissolved carbon from the mantle, metamorphic decarbonization, aqueous dissolution, and soil organic matter are approximately 2%, 82%, 6%, and 10%, respectively. Approximately 87% ± 5% of CO2 generated in the deep crust precipitates as calcite, while approximately 5.5% ± 1% of this carbon is converted to biomass through microbial chemosynthesis at depths less than 2 km. Our random forests approach yielded a metamorphic carbon flux from the entire Himalayan orogenic belt of approximately 2.7 ∼ 4.5 × 1012 mol/yr. The minor CO2 released into the atmosphere (2.5 ∼ 4.2 × 1011 mol/yr) is comparable to the carbon consumption driven by Himalayan weathering. This paper provides new insights into deep carbon cycling, notably that approximately 93% of deeply sourced carbon is trapped in the shallow crust, rendering orogenic processes carbon neutral and possibly acting as one of the major triggers of long-term climate cooling in the Cenozoic.

The 2018 crisis of Kīlauea volcano stands as a major event in its evolution with a large down-rift effusive eruption that drained a shallow magma reservoir at the summit. The characterization of such active magmatic systems and associated hazardous events remains a necessity and a challenge. The summit area is hydrothermally active and strongly altered as indicated by geological mapping. A unique data set of geophysical measurements was collected around Halemaʻumaʻu crater before its collapse. Magnetic data are interpreted here in combination with geological information, temperature anomalies at the surface, self-potential measurements, and a model of electrical conductivity. 3D forward modeling shows that the main magnetic dipole-like anomaly observed around the crater is not only caused by the crater topography but suggests the presence of an important volume of weakly magnetic material beneath the crater, which may be caused by higher temperature and/or hydrothermal alteration. 3D inversion of the data allows us to explore the first order geometry of the magnetic structures. We complement this inversion with 2D forward modeling in order to refine the geometry of major structures. This study shows the presence of major geological structures in the 2018 collapsed area that may have been associated with mechanical weaknesses and could have played a role in the geometry of the collapse. Therefore, mapping magnetic anomalies and monitoring their temporal evolution are of great interest for constraining the nature and mechanical properties of the underlying formations and their temporal evolution in order to help predict future behavior.

Ultra-low velocity zones (ULVZs) are anomalous structures, generally associated with decreased seismic velocity and sometimes an increase in density, that have been detected in some locations atop the Earth's core-mantle boundary (CMB). A wide range of ULVZ characteristics have been reported by previous studies, leading to many questions regarding their origins. The lowermost mantle beneath Antarctica and surrounding areas is not located near currently active regions of mantle upwelling or downwelling, making it a unique environment in which to study the sources of ULVZs; however, seismic sampling of this portion of the CMB has been sparse. Here, we examine core-reflected PcP waveforms recorded by seismic stations across Antarctica using a double-array stacking technique to further elucidate ULVZ structure beneath the southern hemisphere. Our results show widespread, variable ULVZs, some of which can be robustly modeled with 1-D synthetics; however, others are more complex, which may reflect 2-D or 3-D ULVZ structure and/or ULVZs with internal velocity variability. Our findings are consistent with the concept that ULVZs can be largely explained by variable accumulations of subducted oceanic crust along the CMB. Partial melting of subducted crust and other, hydrous subducted materials may also contribute to ULVZ variability.

Many natural faults are believed to consist of velocity weakening (VW) patches surrounded by velocity strengthening (VS) sections. Numerical studies routinely employ this framework to study earthquake sequences including repeating earthquakes. In this laboratory study, we made a VW asperity, of length L, from a bare Poly(methyl methacrylate) PMMA frictional interface and coated the surrounding interface with Teflon to make VS fault sections. Behavior of this isolated asperity was studied as a function of L (ranging from 100 to 400 mm) and the critical nucleation length, h∗ ${h}^{\ast }$, which is inversely proportional to the applied normal stress (2–16 MPa). Consistent with recent numerical simulations, we observed aseismic slip for L/h∗ $L/{h}^{\ast }$ < 2, periodic slip for 2 < L/h∗ $L/{h}^{\ast }$ < 6, and non-periodic slip for 10 < L/h∗ $L/{h}^{\ast }$. Furthermore, we compared the experiments where L was contained by VS material to standard stick-slip events where L was bounded by free surfaces (i.e., L = the total sample length). The free surface case produced ∼10 times larger slip during stick-slip events compared to the contained fault ruptures, even with identical L/h∗ $L/{h}^{\ast }$. This disparity highlights how standard, complete-rupture stick-slip events differ from contained events expected in nature, due to both the free surface conditions and the heterogeneous normal stress along the fault near the free ends, as confirmed by Digital Image Correlation analysis. This study not only introduces the Teflon coating experimental technique for containing laboratory earthquake ruptures, but also highlights the utility of L/h∗ $L/{h}^{\ast }$ as a predictive parameter for earthquake behavior.

Changes in Earth's magnetic field during the Cretaceous Normal Superchron (CNS) spanning ∼121 Ma to ∼84 Ma hold important clues about the geodynamo evolution. Canonical models predict a persistently strong geomagnetic field with low variability during CNS, which, however, has not been observed in the available absolute paleointensity data and seafloor marine magnetic anomaly (MMA) records. The lack of relative paleointensity (RPI) data across CNS further impedes tests of model predictions. Here, we present a ∼9-Myr (∼94–∼85 Ma) RPI record from a Turonian to Santonian hemipelagic succession from IODP Site U1512 offshore southern Australia. Detailed paleomagnetic and rock magnetic analyses demonstrate that the ratio of natural remanent magnetization (NRM) demagnetized at 20 mT over magnetic susceptibility (MS), that is, NRM20mT/MS, as a reliable proxy for the RPI of the Upper Cretaceous succession. The new RPI record shows marked changes in both intensity and variability at ∼90.8 Ma. Also, the 6 Myr-long (∼94–∼88 Ma), near-continuous, ∼1.2 kyr-resolution RPI record exhibits a strong positive correlation between field intensity and variability. Assuming this correlation holds for the entire CNS, an extrapolated RPI curve for the entire CNS is obtained by integrating the positive correlation with field variability estimates from the MMA data. The extrapolated RPI curve shows a strong and highly variable field in the middle CNS but a weak and stable field at its beginning and ending. These features imply a much more dynamic geodynamo than previously thought, and provide crucial benchmarks for unraveling the geodynamo evolution during CNS.

I employ an elasticity-based method to invert a geodetically derived surface velocity field in the western US using for present-day surface strain rate fields with uncertainties. The method uses distributed body forces in a thin elastic sheet and allows for discontinuities in velocity across creeping faults using the solution for dislocations in a thin elastic plate. I compare the strain rate fields with previously published stress orientations and moment rates from geological slip rate data and previous geodetic studies. Geologic and geodetic moment rates are calculated using slip rate and off-fault strain rates from the 2023 US National Seismic Hazard Model (NSHM) deformation models. I find that computed total geodetic moment rates are higher than NSHM summed moment rates on faults for all regions of the western US except the highest deforming rate regions including the Western Transverse Ranges and the northern and southern San Andreas Fault (SAF) system in California. Computed geodetic moment rates are comparable to the moment rates derived from the geodetically based NSHM deformation models in all regions. I find systematic differences in orientations of maximum horizontal shortening rate and maximum horizontal compressive stress in the Pacific Northwest region and along much of the SAF system. In the Pacific Northwest, the maximum horizontal stress orientations are rotated counterclockwise 40–90° relative to the maximum horizontal strain rate directions. Along the SAF system, the maximum horizontal stresses are rotated systematically 25–40° clockwise (closer to fault normal) relative to the strain rates.

The collisional history between Greater India and the Eurasian plate has been well constrained by the study of exhumed Barrovian metamorphic sequence (BMS) rocks in the Himalayan Range. However, in the southeastern Tibetan Plateau, the collisional records have been obscured by intense, regional-scale strike-slip overprinting and recrystallization. Here, in BMS rocks from the Ailao Shan–Red River shear zone (ARSZ), we report the first discovery of a >250 km long, high-pressure (high-P) granulite belt (>1.0 GPa), identified by the presence of relict kyanite and associated decompression reaction textures. Petrological phase equilibrium modeling showed that exposed micaschists in the region represent exhumed middle crust (20–25 km, 600–670°C), while the high-P granulite rocks are remnants of thickened lower crust (45–55 km, 800–850°C). This indicates that the northeast edge of the ARSZ experienced an additional ∼25 km of uplift and exhumation compared to the southwest side, facilitated by brittle thrusting/imbrication along the Ailao Shan fault (micaschists) and ductile extrusion along the Red River fault (granulite). Geochronological study shows that the upper portion of the BMS preserves older metamorphic ages (52–34 Ma) than the lower portion (32–29 Ma), which was attributed to spatial variation in cooling rates. Using calculated P–T–t–d paths, we also examined variation in density and seismic wave speeds for BMS in the ARSZ. Our data correlate with fieldwork conducted elsewhere within the Himalayan Range indicating that the kyanite to sillimanite transition zone may serve as a “cap” for the horizontal migration of melt within the lower crust.

Elevated pore fluid pressure is proposed to contribute to slow earthquakes along shallow subduction plate boundaries. However, the processes that create high fluid pressure, disequilibrium compaction and dehydration reactions, lead to different effective stress paths in fault rocks. These paths are predicted by granular mechanics frameworks to lead to different strengths and deformation modes, yet granular mechanics do not predict their effects on fault stability. To evaluate the role of fluid overpressure on shallow megathrust strength and slip behavior, we conducted triaxial shear experiments on chlorite and celadonite rich gouge layers. Experiments were conducted at constant temperature (130 and 100°C), confining pressure (130 and 140 MPa), and pore fluid pressures (between 10 and 120 MPa). Fluid overpressure due to disequilibrium compaction was simulated by increasing confining and pore fluid pressure synchronously without exceeding the target effective pressure, whereas overpressure due to dehydration reactions was simulated by first loading the sample to a target isotropic effective pressure and then increasing pore fluid pressure to reduce the effective pressure. We find that the effects of fluid pressure and stress path on the mechanical behavior of the chlorite and celadonite gouges can generally be described using the critical state soil mechanics (CSSM) framework. However, path effects are more pronounced and persist to greater displacements in chlorite because its microstructure is more influenced by stress path. Due to its effects on microstructure, the stress path also imparts greater control on the rate-dependence of chlorite strength, which is not predicted by CSSM.

The Seiland Igneous Province (SIP) is a large province of mafic and ultramafic (UM) complexes interpreted to be relics of a giant plumbing system feeding the Ediacaran Central Iapetus Magmatic Province. The Reinfjord Ultramafic Complex (RUC) is one of the four major ultramafic complexes of the SIP. The RUC has a younger dunite core surrounded by wehrlite and lherzolite embedded in country rocks consisting of layered gabbros with sub-horizontal layering and metamorphosed sedimentary rocks. Here, we develop a 3D subsurface model for the RUC using high-resolution magnetic and gravity data and extensive petrophysical measurements from oriented surface samples and drill core samples. Our model indicates that the RUC narrows in depth, extending a minimum of 1.4 km below sea level, and plunges eastwards below the country rock. This model allows us to decipher the lithologic heterogeneities, and the depth and lateral extent of ultramafic rocks, which we interpret in the context of the geologic history of the area. The RUC is spatially separated from other UM complexes of the SIP and the result of this study indicates a smaller depth extent. Combining these findings with the previously reported distribution of the SIP rocks based on the regional gravity data, we propose that the uplift of the crustal block hosting the RUC is larger than for ultramafic complexes in the northwestern part of the SIP.

We report laboratory experiments and numerical simulations demonstrating that the anisotropic characteristics of rocks play a major role in the elongation of hydraulic fractures (HFs) propagating in a plane perpendicular to the rocks' inherent layering (the bedding planes in sedimentary rocks and foliation planes in metamorphic rocks). Transverse anisotropy leads to larger HF extension in the parallel-to-layers/divider direction compared to the perpendicular-to-layers/arrester direction. This directly promotes vertical containment of HFs in most sedimentary basins worldwide even in the absence of any favorable in-situ stress contrasts or other material heterogeneities. More importantly, the ratio of the energy dissipated in fluid viscous flow in the fracture to the energy dissipated in the creation of new surfaces is found to play a critical role on fracture elongation, with fracture-energy dominated HFs being the most elongated while the viscous dominated ones remain more circular. These results open the door to a better engineering and control of HFs containment at depth in view of the competition between material anisotropy (both elastic stiffnesses and fracture toughness anisotropy) and injection parameters (fluid viscosity and rate of injection).

Genuine non-double-couple (non-DC) components of a seismic source, defined here as the non-DC components that are not due to summation of pure double-couple (DC) components, provide important insight into special physical processes in non-earthquake sources such as explosion, volcano eruption and collapse etc. Yet they remain challenging to be resolved. To address the issue and explore the physical mechanism of those special events, we develop a waveform-polarity-based moment tensor (WPMT) inversion method and employ it to study physical process in the 2014–2015 Bárðarbunga volcano event sequence. The WPMT method incorporates P-wave polarity data and seismic waveforms in the source inversion, designs a source simplicity test to check possible complex rupture in the seismic source, and employs a simulated annealing algorithm to search the best source solution. The simplicity test checks consistency of the source processes in the initiation stage of the event and the major energy release process of the event, thus ensuring that the inferred non-DC components are genuine to the seismic source. Real event and synthetic tests indicate that the WPMT method can identify and resolve genuine non-DC components in a seismic source. The WPMT inversions of the Bárðarbunga sequence yield many genuine non-DC source components and reveal that the eruptions are accompanied by seismic activities in depths of 1–5 km with magma migrations out of chambers, collapses of conduits, failures of normal faults, and a magma recharge at a depth of 9 km accompanied by a failure on a nearby normal fault.

We construct a thin-sheet dynamical model of the New Zealand plate boundary that includes faults. Our model fits fault slip rates, style of distributed deformation, and is constrained by relative plate boundary motion. We assume a pseudo-plastic rheology and achieve a best fit to slip rate observations with a deviatoric stress magnitude of 20 MPa. Modeled local forces are significant at Puysegur and Hikurangi subduction zones, and smaller forces are related to mantle downwelling beneath South Island and Havre Trough mantle upwelling. Modeled tractions on faults are mostly 5–20 MPa, similar to or slightly smaller than stress magnitudes adjacent to faults. Modeled shear tractions are generally 2–10 MPa, comparable to stress drops during earthquakes. Modeled stress orientations and fault dips suggest that many faults are not optimally oriented for their style of faulting. Notably small traction magnitudes of <5 MPa and shear tractions of <0.5 MPa are modeled for faults in the central North Island Dextral Fault Belt (NIDFB), which we infer to be very poorly oriented. Friction coefficients on faults (ratio of shear stress to effective normal stress) are in the range 0.1–0.3 for major crustal faults such as the Alpine Fault and Marlborough faults, but subduction zones and the NIDFB have values <0.1. We propose that low values of long-term fault strength, shear stress resolved onto the fault, and overall magnitudes of deviatoric stress in the crust are a consequence of dynamic weakening of faults during fault slip.

Paleomagnetic secular variation (PSV) records provide important information for the dynamic processes of the Earth's geomagnetic field, and also can be used for regional stratigraphic correlation. We conducted a paleomagnetic study on a high sedimentation rate Holocene loess section (the Minle section) with precise 14C age constraints in Northwest China. Rock magnetic results indicate that single domain and pseudo-single-domain magnetites are the main magnetic carriers of the natural remanent magnetization. Combining the alternating field demagnetization and thermal demagnetization results, we obtained a high-resolution PSV record with both absolute declination and inclination since ∼13 ka. Combining with previously published sedimentary records in East Asia, we developed an East Asian PSV reference curve over the past 14 ka (EASed14k) by relocating all the data to a common location (30ºN, 108ºE) and averaging them through Fisher statistics. We further generated a regional geomagnetic directional model spanning the past 14 ka for East Asia (SCHA.EAS14k), using the in situ PSV records and applying the revised spherical cap harmonic analysis in space and cubic B-spline in time. The predicted directional curves from this model at the common location are consistent with the EASed14k curves, which represent the first-order variation of the geomagnetic field in this area. The established reference curve and regional model for East Asia in this study will be used for regional and global comparison of the geomagnetic field.

Open conduit volcanoes erupt with the highest frequency on Earth. Their activity is characterized by an outgassing flux that largely exceeds the gas that could be released by the erupted magma; and by frequent small explosions intercalated by larger events that pose a significant risk to locals, tourists, and scientists. Thus, identifying the signs of an impending larger explosion is of utmost importance for the mitigation of volcanic hazard. Larger explosive events have been associated with the sudden ascent of volatile rich magmas, however, where and why magma accumulates within the plumbing system remains unclear. Here we show that the interaction between CO2-rich fluids and magma spontaneously leads to the accumulation of volatile-rich, low density and gravitationally unstable magma at depth, without the requirement of permeability barriers. CO2-flushing forces the exsolution of water and the increase of magma viscosity, which proceeds from the bottom of the magma column upward. This rheological configuration unavoidably leads to the progressive thickening of a gas-rich and low density (i.e., gravitationally unstable) layer at the bottom of the feeding system. Our calculations account for observations, gas monitoring and petrological data; moreover, they provide a basis to trace the approach to deeply triggered large or paroxysmal eruptions and estimate their size from monitoring data. Our model is finally applied to Stromboli volcano, an emblematic example of open conduit volcano, but can be applied to any other open conduit volcano globally and offers a framework to anticipate the occurrence of unexpectedly large eruptions.

In South Iceland, populated and agricultural areas are at risk of earthquakes due to their location within the South Iceland Seismic Zone (SISZ). In 2008, two moderate-sized earthquakes (M5.8 and M5.9) occurred in Ölfus, the western end of this highly active transform zone. We analyze temporal seismic velocity variations (dv/v) related to the Ölfus earthquake doublet, using cross-correlations of ambient noise in the frequency range of 0.1–3.0 Hz. The two mainshocks decrease the average velocity by 0.8% at the nearest stations. The co-seismic changes are most noticeable from 0.7 to 1.7 Hz and affect a 40 km wide region. We present a first-time comparison of dv/v to crustal deformation, seismicity, co-seismic volumetric stress changes and reported PGA distribution for the Ölfus doublet. Ground accelerations caused by mainshocks at intermediate distances suggest that strong shaking-related damage may contribute to the co-seismic dv/v decrease. A rapid velocity increase (0.3%) in a month after the co-seismic drop indicates crustal rock healing. We find 3-months of post-seismic decorrelation, followed by a nearly permanent velocity decrease (0.2%) confined to a shallow layer (1 km) until the end of the observation period. Afterslip and pore fluid effects in the near-source region are likely to influence post-seismic dv/v. We demonstrate that seismic interferometry can contribute to future fault-zone monitoring operations in the SISZ by detecting small changes in velocity.

We develop a method to estimate relative seismic moments M 0 and corner frequencies f c of acoustic emission events recorded in laboratory experiments from amplitude spectra of signal's coda composed of reverberated and scattered waves. This approach has several advantages with respect to estimations from direct waves that are often clipped and also are difficult to separate in experiments performed on small samples. Also, inversion of the coda spectra does not require information about the source locations and mechanisms. We use the developed method to analyze the data of two experiments: (a) on granite from the Voronezh crystal massif and (b) on Berea sandstone. The range of absolute corner frequencies estimated in both experiments is around 70 − −700 kHz. The range of relative seismic moments covers 103.5. The relation between f c and M 0 observed on the first stages of both experiments, consisted of increasing isotropic confining pressure, approximately follow M0∼fc−3 ${M}_{0}\sim {f}_{c}^{-3}$ scaling and the b-value of the Gutenberg-Richter distribution was found close to 1. This can be interpreted as rupturing of preexisting material defects with a nearly constant stress-drop and has a similarity with observations of “natural” earthquakes. Deviations from this “earthquake-like” behavior observed after applying axial loading and initiation of sample damaging can be interpreted as changes in stress-drop. Lower stress-drops prevail for sandstone and higher for granite sample respectively that can be related to the strength of corresponding material.

Pockmarks are morphological depressions commonly observed in ocean and lake floors. Pockmarks form by fluid (typically gas) seepage thorough a sealing sedimentary layer, deforming and breaching the layer. The seepage-induced sediment deformation mechanisms, and their links to the resulting pockmarks morphology, are not well understood. To bridge this gap, we conduct laboratory experiments in which gas seeps through a granular (sand) reservoir, overlaid by a (clay) seal, both submerged under water. We find that gas rises through the reservoir and accumulates at the seal base. Once sufficient gas over-pressure is achieved, gas deforms the seal, and finally escapes via either: (a) doming of the seal followed by dome breaching via fracturing; (b) brittle faulting, delineating a plug, which is lifted by the gas seeping through the bounding faults; or (c) plastic deformation by bubbles ascending through the seal. The preferred mechanism is found to depend on the seal thickness and stiffness: in stiff seals, a transition from doming and fracturing to brittle faulting occurs as the thickness increases, whereas bubble rise is preferred in the most compliant, thickest seals. Seepage can also occur by mixed modes, such as bubbles rising in faults. Repeated seepage events suspend the sediment at the surface and create pockmarks. We present a quantitative analysis that explains the tendency for the various modes of deformation observed experimentally. Finally, we connect simple theoretical arguments with field observations, highlighting similarities and differences that bound the applicability of laboratory experiments to natural pockmarks.

Large Low Velocity Provinces (LLVPs) are situated oppositely in the lowermost mantle beneath the Pacific Ocean and Africa. Deciphering the detailed seismic structures at the edge of LLVPs can provide key information on the composition and dynamics in the deep Earth. Here, we provide a detailed seismic image at the western edge of the Pacific LLVP by dense recordings. Differential travel time residuals and amplitude ratios between ScS and S outline the S-wave western boundary of the Pacific LLVP, suggesting the complex structures including low/high-velocity patches in the lowermost mantle in our study region. We determine the 3D low-velocity structure by modeling the delayed ScS and high-velocity D″ layer structure by modeling the anomalous Scd, with tight constraints from multiple events data. The drastically varied waveforms in azimuth suggests a sharp transitional boundary among the complex structures. After comparing the velocity structures in adjacent regions, we propose that the 3D structures of the western edge of the Pacific LLVP are strongly influenced by the vigorous mantle flow associated with the actively subducted slab.

Basaltic fissure eruptions, which are the most common type of eruption on Earth, are fed by dykes which mediate magma transport through the crust. Dyke propagation processes are important because they determine the geometry of the transport pathway and the nature of any geophysical signals associated with magma ascent. Here, we investigate small-scale (mm–cm wide) banding features at the margins of dykes in the Teno Massif (Tenerife, Spain) and the Columbia River Basalt Province (CRBP) (USA). Similar marginal bands have been reported for dykes in numerous localities around the world. Dyke margins record valuable information about propagation because they are the first material to solidify against the host rock at the propagating dyke tip. We find that the marginal bands are defined by cyclic variations in phenocryst concentration and vesicularity, and we infer that these cyclic variations in texture are a product of cyclic variations in magma flow rates and pressures within the dyke tip. This indicates that dyke emplacement occurs in pulses, with propagation repeatedly hindered by the rapid cooling and solidification of magma in the narrow dyke tip. Using a 1D conduction model, we estimate the time taken for each band to cool and solidify, which provides a timescale of several minutes to tens of minutes for the pulses. The occurrence of similar bands in various volcanic settings suggests that pulsatory propagation is a common, if not ubiquitous, process associated with dyke emplacement.

All instrumented basaltic caldera collapses have generated Mw > 5 very long period earthquakes. However, previous studies of source dynamics have been limited to lumped models treating the caldera block as rigid, leaving open questions related to how ruptures initiate and propagate around the ring fault, and the seismic expressions of those dynamics. We present the first 3D numerical model capturing the nucleation and propagation of ring fault rupture, the mechanical coupling to the underlying viscoelastic magma, and the associated seismic wavefield. We demonstrate that seismic radiation, neglected in previous models, acts as a damping mechanism reducing coseismic slip by up to half, with effects most pronounced for large magma chamber volume/ring fault radius or highly compliant crust/compressible magma. Viscosity of basaltic magma has negligible effect on collapse dynamics. In contrast, viscosity of silicic magma significantly reduces ring fault slip. We use the model to simulate the 2018 Kı̄lauea caldera collapse. Three stages of collapse, characterized by ring fault rupture initiation and propagation, deceleration of the downward-moving caldera block and magma column, and post-collapse resonant oscillations, in addition to chamber pressurization, are identified in simulated and observed (unfiltered) near-field seismograms. A detailed comparison of simulated and observed displacement waveforms corresponding to collapse earthquakes with hypocenters at various azimuths of the ring fault reveals a complex nucleation phase for earthquakes initiated on the northwest. Our numerical simulation framework will enhance future efforts to reconcile seismic and geodetic observations of caldera collapse with conceptual models of ring fault and magma chamber dynamics.