JGR–Solid Earth

The Hengill geothermal field, located in southwest Iceland, is host to the Hellisheiði power plant, with its 40+ production wells and 17 reinjection wells. Located in a tectonically active area, the field experiences both natural and induced seismicity linked to the power plant operations. To better manage the risk posed by this seismicity, the development of robust and informative forecasting models is paramount. In this study, we compare the forecasting performance of a model developed for fluid-induced seismicity (the Seismogenic Index model) and a class of well-established statistical models (Epidemic-Type Aftershock Sequence). The pseudo-prospective experiment is set up with 14 months of initial calibration and daily forecasts for a year. In the timeframe of this experiment, a dense broadband network was in place in Hengill, allowing us to rely on a high quality relocated seismic catalog. The seismicity in the geothermal field is characterized by four main clusters, associated with the two reinjection areas, one production area, and an area with surface geothermal manifestations but where no operations are taking place. We show that the models are generally well suited to forecast induced seismicity, despite some limitations, and that a hybrid ETAS model accounting for fluid forcing has some potential in complex regions with natural and fluid-induced seismicity.

Compressive stress generated at collision fronts can propagate over long distances, inducing deformation within the continent's interior. Nevertheless, the factors governing the partitioning of intracontinental deformation remain enigmatic. The Altai Mountains serve as a type-example of ongoing intracontinental deformation. Here, we investigate the crustal architecture of the Chinese Altai Mountains, using receiver functions obtained from newly deployed dense seismic nodal arrays. The new seismic results reveal distinct crustal features, including (a) a negative polarity discontinuity beneath Chinese Altai Mountains, suggesting a low-velocity layer; (b) a north-dipping mid-crustal structure beneath the suture zone between East Junggar and Chinese Altai, indicating underthrusting of East Junggar's lower crust beneath the Chinese Altai Mountains; (c) a double Moho structure beneath East Junggar, revealing a high-velocity lower crustal layer. In conjunction with constraints from previous multi-disciplinary regional studies, the double Moho structures are interpreted as mafic restite from Late Paleozoic magma underplating. The addition of mafic materials can significantly enhance the rheological strength of East Junggar's crust, causing it to function as an indenter that thrust beneath the Chinese Altai Mountains during the subsequent convergence process. As a consequence, significant deformation occurs in the Chinese Altai region, resulting in the emergence of decollements, as evident by the negative polarity discontinuity. The presence of pre-existing decollements makes the Altai Mountains region more susceptible to deformation, thereby facilitating the concentration of intracontinental deformation. These findings illuminate the evolution history of the Chinese Altai Mountains and highlight the great impacts of ancient tectonics on intracontinental deformation partitioning.

Microcracking due to thermal stresses affects the mechanical and flow properties of rocks, which is significant for thermally dynamic environments such as volcanoes and geothermal reservoirs. Compared with other crustal rocks like granite, volcanic rocks have a complex and variable response to temperature; it remains unclear how thermal microcracks form and how they are affected by temperature. We heated and cooled samples of low-porosity basalts containing different amounts of microcracks and a porous andesite over three cycles, whilst monitoring microstructural changes by acoustic emission (AE) monitoring and measurement of P-wave velocity (v P ; up to 450°C) and thermal expansion coefficient (TEC; up to 700°C). During the second and third cycles, the TEC was positive throughout and the rate of detected AE was low. In contrast to studies on granite, we measured a strong and reversible increase in v P with increasing temperature (by 15%–40% at 450°C), which we interpret as due to microcrack closure. During the first cycle, AE and v P measurements indicated thermal microcracking within the andesite and the basalt with a low initial microcrack density. For these samples, strong inflexions in the TEC indicated stress relaxation during heating, preceding significant thermal microcracking during cooling. The basalt with a high initial microcrack density underwent little microcracking throughout all cycles. Our results and a review of the literature relate the initial microstructure to the occurrence of thermal microcracking and explore the potentially significant influence of temperature on volcanic rock properties.

Inflation of viscous magma intrusions in Earth's shallow crust often induces strain and fracturing within heterogeneous host rocks and dome-shaped ground deformation. Most geodetic models nevertheless consider homogeneous, isotropic, and linear-elastic media wherein stress patterns indicate the potential for failure, but without simulating actual fracturing. We present a two-dimensional Discrete Element Method (DEM) application to simulate magma recharge in a pre-existing laccolith intrusion. In DEM models, fractures can propagate during simulations. We systematically investigate the effect of the host rock toughness (resistance to fracturing), stiffness (resistance to deformation) and the intrusion depth, on intrusion-induced stress, strain, displacement and spatial fracture distribution. Our results show that the spatial fracture distribution varies between two end-members: (a) for high stiffness or low toughness host rock or a shallow intrusion: extensive cracking, multiple vertical surface fractures propagating downward and two inward-dipping highly cracked shear zones that connect the intrusion tip with the surface; and (b) for low stiffness or high toughness host rock or a deeper intrusion: limited cracking, one central vertical fracture initiated at the surface, and two inward-dipping fractures at the intrusion tips. Abrupt increases in surface displacement magnitude occur in response to fracturing, even at constant magma injection rates. Our modeling application provides a novel approach to considering host rock mechanical strength and fracturing during viscous magma intrusion and associated dome-shaped ground deformation, with important implications for interpreting geodetic signals at active volcanoes and the exploitation of geothermal reservoirs and mineral deposits.

Altered oceanic crust (AOC) plays a critical role in geochemical recycling in subduction zones. However, identifying contributions of subducted AOC to arc magmas remains a conundrum due to the lack of effective tracers. Here, we investigate the Ba-Sr-Nd isotopic compositions of lavas from the Mariana arc and back-arc. Based on a statistical analysis of the Sr-Nd isotopes for global arc volcanoes, we confirm that AOC-derived fluid (or hydrous melt), rather than sediment-derived melt or fluid, is responsible for the Sr-Nd isotope decoupling (i.e., 87Sr/86Sr is “excessively” enriched relative to 143Nd/144Nd when compared to the “normal” mantle derivates) observed in island arc lavas. We show that the arc lavas with increasingly decoupled Sr-Nd isotopes generally have heavier Ba isotope ratios, which is also a characteristic feature of AOC-derived fluids. Thus, these results establish an intimate link between subducted AOC, heavy Ba isotope compositions, and Sr-Nd isotope decoupling signature in island arcs, which provides a powerful tool to trace the AOC recycling in subduction zones. Furthermore, a similar correlation is observed between Sr-Nd isotope decoupling and heavy B isotope ratios for global arc lavas, implying that the recycling of AOC component is generally linked to serpentinite dehydration in subduction zones.

Phase relation and equation of state of iron-titanium oxyhydroxide with α-PbO2-type crystal structure (hydrous α-phase) was investigated at pressure from middle upper mantle to mantle transition zone. The experiments were performed using multi-anvil apparatus with in situ synchrotron X-ray diffraction. The metal-diamond sample container was used to maintain a closed system with respect to water while allowing X-ray transmittance. Starting materials were mixture of reagent grade goethite and anatase with Fe:Ti = 1:1 and 1:3. The hydrous α-phase was synthesized at 11–14 GPa and 800–900°C. The X-ray diffraction data were obtained over a wide range of pressures from 5 to 22 GPa and temperatures from room temperature to 930°C. Using pressure-volume-temperature data of the hydrous α-phase with Fe:Ti = 1:1 collected at room temperature to 800°C and at 8.9–20.8 GPa, we determined the isothermal bulk modulus (K T0 = 183 (1) GPa) and the thermal expansivity (α 0 = 3.29 (1) × 10−5K−1) at ambient condition. The stability field of the hydrous α-phase and phase relation of FeOOH-TiO2 system at 900°C was well constrained. It was found that the hydrous α-phase decomposes into baddeleyite-type TiO2 + ε-FeOOH at pressure of approximately 20–21 GPa, and into ilmenite + rutile at 5–6 GPa. This stability field equivalent to depth of 180–600 km in the subduction zone. Our results suggest that the hydrous α-phase is an important water reservoir in the middle upper mantle to mantle transition zone.

Compaction bands are a type of localized deformation that can occur as diffuse or discrete bands in porous rocks. While modeling of shear bands can replicate discrete and diffusive bands, numerical models of compaction have so far only been able to describe the formation of discrete compaction bands. In this study, we present a new thermodynamic approach to model compaction bands that is able to capture both discrete and diffuse compaction band growth. The approach is based on a reaction-diffusion formalism that includes an additional entropy flux. This entropic velocity regularizes the solution, by introducing a characteristic diffusion length scale and controlling the mode change from discrete to diffusive post-localisation growth. The approach is used to model compaction band growth in highly porous carbonates. The model can replicate the areas of material damage exhibiting reduced porosity which are often observed as nuclei for the growth of compaction bands in experiments. The model also has the versatility to predict the formation of diffuse compaction bands, which is a significant advance in the field of compaction band modeling. The method can potentially be used for investigating the effect of material heterogeneities on compaction band growth and is heuristic for developing new methodologies for forecasting compaction band formation.

Stress-driven melt segregation may have important geochemical and geophysical effects but remains a poorly understood process. Few constraints exist on the permeability and distribution of melt in deformed partially molten rocks. Here, we characterize the 3D melt network and resulting permeability of an experimentally deformed partially molten rock containing several melt-rich bands based on an X-ray microtomography data set. Melt fractions range from 0.08 to 0.28 in the ∼20-μm-thick melt-rich bands, and from 0.02 to 0.07 in the intervening ∼30-μm-thick regions. We simulated melt flow through subvolumes extracted from the reconstructed rock at five length scales ranging from the grain scale (3 μm) to the minimum length required to fully encompass two melt-rich bands (64 μm). At grain scale, few subvolumes contain interconnected melt, and permeability is isotropic. As the length scale increases, more subvolumes contain melt that is interconnected parallel to the melt bands, but connectivity diminishes in the direction perpendicular to them. Even if melt is connected in all directions, permeability is lower perpendicular to the bands, in agreement with the elongation of melt pockets. Permeability parallel to the bands is proportional to melt fraction to the power of an exponent that increases from ∼2 to 5 with increasing length scale. The permeability in directions parallel to the bands is comparable to that for an isotropic partially molten rock. However, no flow is possible perpendicular to the bands over distances similar to the band spacing. Melt connectivity limits sample scale melt flow to the plane of the melt-rich bands.

We calculated the pressure-volume relations of liquids Fe and Fe7C3 up to ∼360 GPa and 3000–8000 K based on the first-principles molecular dynamics simulations. The results demonstrate that liquid Fe7C3 is similar in compressibility and thermal expansivity to liquid Fe at the Earth's core pressure range. We then obtained a thermodynamic model of the liquidus phase relations in the Fe-C system at high pressures by using the thermal equation of state (EoS) of liquid Fe7C3 as well as that of liquid Fe3C estimated by interpolation between the volumes of Fe and Fe7C3 calculated in this study. The previously reported eutectic points in the Fe-Fe3C system are reproduced with interaction parameters that are dependent on temperature but independent on pressure at >∼50 GPa. The melting curve of Fe7C3 should be close to that of Fe, which leads to a change in the eutectic system from Fe-Fe3C to Fe-Fe7C3 above 250 GPa as observed by previous experiments. The thermodynamic model also suggests that the solidus and liquidus of Fe3C are relatively low. The Fe-Fe7C3 liquidus phase relations at the Earth's inner core boundary (ICB) pressure indicate the Fe7C3 inner core is unlikely. While the carbon content in the Earth's outer core may be small, the addition of 1 wt% C to liquid Fe drops the liquidus temperature to crystallize hcp-Fe by ∼450 K, suggesting that the presence of carbon in the outer core is preferable for a relatively low core-mantle boundary (CMB) temperature.

The thermochemical structure of the lithosphere controls melting mechanisms in the mantle, as well as the location of volcanism and ore deposits. Obtaining reliable images of the lithosphere structure, and its complex interactions with the asthenosphere, requires the joint inversion of multiple data sets and their associated uncertainties. In particular, the combination of seismic velocity and electrical conductivity, along with proxies for bulk composition and elusive minor phases, represents a crucial step toward fully understanding large-scale lithospheric structure and melting processes. We apply a novel probabilistic approach for joint inversions of 3D magnetotelluric and seismic data to image the lithosphere beneath southeast Australia. The results show a highly heterogeneous lithosphere with deep conductivity anomalies that correlate with the location of Cenozoic volcanism. In regions where the conductivities have been at odds with sub-lithospheric temperatures and seismic velocities, we observe that the joint inversion provides conductivity values consistent with other observations. The results reveal a strong relationship between metasomatized regions in the mantle and (a) boundaries of geological provinces, elucidating the subduction-accretion process in the region; (b) distribution of leucitite and basaltic magmatism; (c) independent geochemical data, and (d) a series of lithospheric steps which constitute areas prone to generating small-scale instabilities in the asthenosphere. This scenario suggests that shear-driven upwelling and edge-driven convection are the primary mechanisms for melting in eastern Australia, contrary to the conventional notion of mantle plume activity. Our study presents an integrated lithospheric model for southeastern Australia and provides valuable insight into the mechanisms driving surface geological processes.

Contact metamorphism of carbonate rocks in response to fluid infiltration releases carbon dioxide (CO2), potentially affecting Earth's carbon budget over geological timescales. This process was notably active in the Cretaceous when the closure of the Neo-Tethys Ocean led to extensive arc magmatism in southern Tibet, coincident with a greenhouse climate interval. These arc plutons intersected carbonate sequences and formed widespread calcsilicate rocks. The mineral assemblages of these rocks record both the progression and duration of the CO2 release. In this study, we focus on a representative aureole in southern Tibet that was metamorphosed by a ∼64 Ma granitoid pluton of the Gangdese Batholith. The aureole presents lithologic zonation, with the inner garnet-wollastonite zone indicating equilibrium with a water-rich fluid at ∼525°C. Large garnet porphyroblasts feature cross-cutting veins of secondary garnets that illustrate several stages of interface-coupled dissolution-reprecipitation processes facilitated by fluid infiltration. Diffusion modeling for Mg concentration profiles across vein-host garnet yields a brief high-temperature stage (25–150 kyr)—a duration still significantly overestimated considering factors that may accelerate diffusion in garnet, such as lattice parameters and hydration. Such short timescale aligns with conduction modeling of a cooling pluton, which, combined with the bulk-rock mass balance analysis, indicates a local carbon flux comparable to fluxes at other modern tectonic settings. The metamorphic decarbonation, potentially episodic, in the continental arc may have been an important endogenic carbon source in the late Cretaceous.

Large earthquakes can trigger smaller seismic events, even at significant distances. The process of earthquake triggering offers valuable insights into the evolution of local stress states, deepening our understanding of the mechanisms of earthquake nucleation. However, our ability to detect these triggered events is limited by the quality and spatial density of local seismometers, posing significant challenges if the triggered event is hidden in the signal of a nearby larger earthquake. Distributed acoustic sensing (DAS) has the potential to enhance the monitoring capability of triggered earthquakes through its high spatial sampling and large spatial coverage. Here, we report on an uncatalogued magnitude (M) 5.1 event in northeast Turkey, which was likely dynamically and instantaneously triggered by the 2023 M7.8 earthquake in southeast Turkey, located 400 km away. This event was initially discovered on ∼1,100 km of active DAS recordings that are part of an 1,850-km linear array. Subsequent validation using local seismometers confirmed the event's precise time, location, and magnitude. Interestingly, this dynamically triggered event exhibited precursory signals preceding its P arrivals on the nearby seismometers. It can be interpreted as the signal from other nearby, uncatalogued, smaller triggered events. Our results highlight the potential of high-spatial-density DAS in enhancing the local-scale detection and the detailed analysis of earthquake triggering.

Large earthquakes that occur away from plate boundaries (i.e., large intraplate earthquakes), though rare, could cause heavy damage. Understanding their causes could help with seismic hazard assessment. In this study, we image the seismogenic structure under the 1969 Ms6.4 Yangjiang intraplate earthquake within the stable South China Block with a high-resolution three-dimensional shear-wave velocity model constructed from ambient noise tomography based on a dense nodal seismic array. The model images relatively low velocities at a fault intersection region from the surface to at least 13 km depth in the Yangjiang seismogenic zone. Given spatial links among the intersection zone of faults, surface seawater, and low velocities, we interpret the low velocities to reflect the presence of seawater-filled, highly fractured rocks created by the intersection of faults. We infer that the seawater infiltration through densely-fractures could have elevated pore pressure in the long term, as evidenced by low earthquake b values at 8–13 km depths. We hypothesize that this long-term seawater infiltration could have lowered the stress threshold of earthquake occurrence and contributed to the generation of the 1969 Ms6.4 Yangjiang earthquake. We propose that the large intraplate earthquakes within stable plates tend to occur in regional weakening regions (e.g., fault intersections) favorable to stress buildup. Given sufficient tectonic stress accumulation, the long-term hydrologically driven crustal stress variations could play an essential role in triggering large intraplate earthquakes. This study could potentially contribute to improving seismic hazard assessment in the Guangdong-Hong Kong-Macao-Greater-Bay-Area, China.

Active seismic surveys are routinely employed by academia to study geological structure of the crust and upper mantle. Wavefields generated during these surveys are sampled at the receiver locations, but the wave-paths traveled from a source to a sensor remains unknown. Although seismic acquisition layouts designed to investigate complex crustal-scale environments are often two-dimensional, the seismogram recorded at the receiver location represents information gathered along the three-dimensional wavepaths that might offset from the 2D source/receiver profile along its transverse direction. This so-called 3D-effect distorts the results of 2D seismic imaging, which is unable to handle the out-of-plane propagation. Despite the numerous 2D seismic imaging case studies, the assessment of this issue is often overlooked. However, the problem exists - especially for crustal-scale profiles, where seismic energy propagates over distances of hundreds of kilometers and probes different crustal units. In this work we investigate the impact of 3D-effect on the results of 2D velocity model building from the academic ocean-bottom seismometer data. We show with polarization analysis how the 3D-effect can manifest itself in the data domain. Using various scenarios of acquisition we evaluate the imprint of the out-of-plane propagation on the data and the results of full-waveform inversion. We show that 2D velocity model building from the seismic profiles acquired in the complex geological setting can lead to wrong solution. Looking for the remedy to this issue we couple different configurations of acquisition geometries with 3D full-waveform inversion that allow to handle the 3D effect and provide correct model reconstruction.

Numerous studies have reported that surface hydrological loading can seasonally modulate seismicity rates at crustal depths. For example, substantial winter snow accumulation occurs across the Japanese Islands, and these snowy regions appear to have seasonally modulated the occurrence of previous large inland earthquakes. Therefore, it is important to investigate the impact of seasonal stress changes on crustal seismicity to deepen our understanding of earthquake generation. Here we constrain seasonal changes in the surface load across northeastern Japan using Global Navigation Satellite System surface displacements and evaluate the potential relationship between temporal trends in inland seismicity and estimated seasonal stress changes. The spatial distribution of the seasonal surface load is consistent with snow depth along the Sea of Japan. The inland seismicity beneath northeastern Japan is modestly modulated by the seasonal stress changes that are induced by the annual snow load. However, this seasonal response is weaker than that in other regions. This weak modulation may be due to the small surface-load-induced stress perturbation relative to the long-term-averaged stressing rate and/or the limited presence of crustal fluids to trigger seismicity in Japan.

To integrate structural subsurface models and smooth seismic velocity models, they need to share common features and resolutions. Here, we propose a new approach, Depth Assessment from Rayleigh Wave Ellipticities (DARE), for estimating the depth of sudden velocity changes from ambient-noise multi-mode Rayleigh waves applicable to a wide range of frequencies. At frequencies where multi-mode Rayleigh waves have an extremum in ellipticity, the phase velocity can be used to estimate the depth of sudden velocity changes. We test our approach theoretically, numerically, and on real data from two geothermal sites by extracting Rayleigh wave ellipticities and phase velocities from three-component beamforming of ambient noise using the python code package B3AMpy. For a small-scale array, our approach validates the depth of quaternary sediments predicted by geological models. For deeper velocity changes, high uncertainties remain but the general trend of inclining boundaries can be recovered well. We demonstrate that, if impedance contrasts are larger than three, our approach is valid for multiple layers, laterally heterogeneous models, and a wide range of Poisson ratios.

The purpose of this study was to prove the direct correlation of a successions of gravel-clay beds recovered in borehole with the melt-water pulses associated with the sea-level oscillations indicated in the δ 18 O record in the time-span preceding the 100-kyr glacial cycles. Aimed at this scope, we provide combined 40 Ar/39 Ar and paleomagnetic constraints to a set of seven aggradational successions recovered from a 120 m deep borehole drilled in the buried Paleo-Tiber delta in Rome (central Italy). The geochronologic constraints enable the correlation of each aggradational succession, characterized by a sudden transition from coarse gravel at the base to sandy clay sediments, with periods of sea-level rise indicated by the δ 18 O curve encompassing MIS 37 through MIS 19, from 1250 to 780 ka. This stratigraphy, provides a unique and unprecedented well-dated evidence of glacial/deglacial events, matching the global benthic δ 18 O stack during this time frame. Furthermore, this study validates the hypothesis that gravel deposition in the catchment basin and the delta of the main rivers in central Italy is triggered by the melting of glaciers in the Apennines Mountain range. It demonstrates the significant potential of these deglaciation proxies to be used worldwide to unravel the chronology of glacio-eustatic events.

The top of the subducting plate is a thick, complex zone where the heterogeneous structure likely controls earthquake rupture processes. Imaging this heterogeneous channel typically involves active-source methods with limited depth penetration, or low-resolution teleseismic methods. To access short wavelengths at greater depth, we use high-frequency P-to-S (PS, 1–15 Hz) mode-converted arrivals from nearby earthquakes >50 km deep to image the plate interface at vertical scales <1 km. We use 37 broadband stations in southcentral Alaska between 2007 and 2008 at 10–15 km spacing, spanning the great 1964 earthquake rupture zone and adjacent deeper slow slip and tremor regions. The central 21 stations record high-amplitude PS arrivals converting from the megathrust region, at depths corresponding to the top of a prominent low-velocity zone (LVZ) in receiver function images. The PS/P amplitude ratio (APS_P) varies along strike and with depth of the conversion point but is independent of earthquake location and varies slowly between adjacent stations. APS_P changes with slab depth, indicating changes in lithology or fluid content of the plate interface, consistent with transitions in slip behavior from locked to slow slip. High APS_P cannot be explained by a velocity step or a single low-velocity, high Vp/Vs layer, but requires several alternating high and low-velocity layers. These observations indicate that the LVZ is a highly heterogeneous channel at multiple scales, resembling a subduction channel or sheared zone of metasediment and altered crust as observed in many exhumed subduction zones.

Extensive investigation of continental rift systems has been fundamental for advancing the understanding of extensional tectonics and modes of formation of new ocean basins. However, current rift classification schemes do not account for conjugate end members formed by Large Igneous Province crust, referring to thick mafic crust, sometimes including continental fragments. Here, we investigate the rifting of William's Ridge (Kerguelen Plateau) and Broken Ridge, components of the Kerguelen Large Igneous Province now situated in the Southeast Indian Ocean, and incorporate these end members into the deformation migration concept for rifted margins. We use multichannel seismic reflection profiles and data from scientific drill cores acquired on both conjugate margins to propose, for the first time, a combined tectono-stratigraphic framework. We interpret seismic patterns, tectonic features, and magnetic anomaly picks to determine an across-strike structural domain classification. This interpretation considers the rift system overall to be “magma-poor” despite being located proximal to the Kerguelen plume but suggests that syn-rift interaction between the Kerguelen mantle plume and the lithospheric structure of William's Ridge and Broken Ridge has controlled the along-strike segmentation of both conjugates. We integrate seismic reflection and bathymetric data to test the hypothesis of predominantly transform motion, between the Australian and Antarctic plates, in Late Cretaceous and Paleogene time.

Fault reactivations induced by deep excavation can pose significant challenges to underground construction or resource extraction. Laboratory experiments on rock faults demonstrate that unloading-induced fault reactivations obey the Coulomb failure criterion derived from loading-induced events. However, the effect of stress path during unloading on the failure criterion and rupture dynamics of fault reactivations remains poorly understood. Here, we present findings from a series of laboratory experiments aimed at elucidating the effect of the unloading path on the failure criterion and rupture dynamics of fault reactivations. We conducted experiments under various stress conditions, examining two cases of unloading paths. In Case I, we unloaded the minimum principal stress, while in Case II, the maximum principal stress was unloaded. Strain gauges and high-speed photography were employed to capture the transient dynamic rupture process. Our investigations have yielded new insights into the effect of unloading path on the rupture dynamics when the fault is reactivated. In Case I, we observed fault reactivations resembling those loading-induced events characterized by forward sliding. Conversely, in Case II, fault reactivations associated with stress reversal produce mild reversed sliding with lower stress drop and rupture velocity. Furthermore, we find that there is a remarkable reduction in static friction for reversed sliding, indicating that the failure criterion for fault reactivation is influenced by the stress path. We demonstrate that enhanced stress heterogeneity, caused by stress reversal, serves as a mechanism for reduced static friction. These findings contribute to our understanding of the mechanisms underlying fault reactivations, particularly those involving reversed sliding.