Geophysical Journal International

SummaryThe Global Navigation Satellite System—Acoustic ranging combination technique (GNSS-A) is a seafloor geodetic observation technique that achieves an accuracy of centimeters by combining high-rate GNSS data with acoustic ranging. The technique determines the seafloor position by acoustic ranging between the sea surface and multiple seafloor stations, using GNSS data from the sea surface station. Here, the gradient state of the underwater sound speed structure (SSS) is a significant source of error. The open-source software GARPOS can reduce the effect from underwater gradient structures by retrieving the underwater disturbance as a parameter projected onto the sea surface and seafloor. To evaluate the effects of underwater disturbances, a quantitative comparison of the model parameters is necessary. In this study, we developed a representation method to evaluate features of the ocean field. Here, the expression method was described in the order of a formulation and an interpretation in the case of a two-dimensional cross section and extension to the case of an actual three-dimensional field. This method makes it possible to evaluate SSS states in the GNSS-A observations. As an example, we showed the correlation between the anomaly of the expressed ocean state and the anomaly of the seafloor station position, showing that this expression method is an effective index for correcting bias errors. Additionally, we used the data from sites located in the Kuroshio area, a strong current near Japan, to show that the ocean state can be quantitatively interpreted using this expression method.

SummaryThe purpose of this study is to investigate the seismic ambient noise spectral anomalies that occur near gas reservoirs. These anomalies involve a significant spectral amplification of the vertical component for frequencies generally between 1.5 and 4 Hz and have been reported at various hydrocarbon sites worldwide. There are differing views on the mechanisms responsible for these anomalies. The guideline for this study is that many hydrocarbon reservoirs share a common geological feature: an anticline structure. It appears to cause site effects that influence the amplitude of the ambient noise wavefield. This research examines a dense real dataset of ambient noise recorded at the Chémery underground gas storage site in France. The analysis identifies stable spectral anomalies between 1.2 and 2.4 Hz that are correlated to the position of the anticline structure, which also corresponds to the position of the gas bubble. We use a beamforming technique to study the composition and the origin of the ambient noise, and show that the variations of the spectral anomalies over time are correlated to changes in the source wavefield. Finally, we perform numerical simulations of Rayleigh wave propagation within a realistic 3D velocity model of the Chémery site, while using source distributions directly extracted from real data analysis. The comparison of the simulated anomalies with real data yields a satisfactory qualitative fit. We conclude that the fundamental mode Rayleigh wave site effect on the anticline is the main mechanism of the spectral anomaly.

SummaryA new ray-tracing algorithm based on the shortest path method was designed and optimized to model seismic rays. The workflow is based on Dijkstra’s algorithm to find the shortest path iteratively on Self-Adapting Random Grids. These random grids adapt from one iteration to the next, storing in memory those nodes that minimize the travel time in each iteration, consequently automatically densifying the grid in areas with significant velocity variations. Different data structures and grid geometries were studied, and it was confirmed that using a combination of a binary heap and random grids prevents systematic errors generated by using regular grids while keeping the computational times low. Since the nodes that minimize the travel time are stored in each iteration, the memory requirements increase logarithmically, with the initial iterations having the most significant impact on the error minimization but plateauing afterwards. It was found that, albeit this increase in memory requirements, by focusing the densification of nodes in areas close to the actual trajectory of the rays, the error minimization requires fewer iterations than other implementations based on multiple realizations of random grids (requiring less than a tenth of the number of iteration of other methods). A function to handle extreme topographic variations was also designed, which uses a boolean test to reject rays crossing the topography, making the algorithm suitable for first arrival modelling in complex topography areas, such as the required for tomographic inversion of first arrivals in thrust-belt land seismic.

SummaryPore pressure oscillations induced by stress variations including propagating seismic waves from remote earthquakes have been widely observed in various groundwater systems. The monitored pressure change in wells shows significant water-level oscillations to volumetric strain as well as to S and Love waves. Recent observations demonstrated azimuthal dependence of the pore pressure oscillations with respect to stress indicators and fault zone orientation. Within the fault zone, damaged induced anisotropy is the result of the alignment and orientation of cracks and other internal flaws within the rock. In this work, we provide a complete quantitative description of the pore pressure changes induced by passing seismic waves associated with different orientations and values of principal stress and damage tensor components. The model quantifies the azimuthal dependency of the pore pressure response by a non-dimensional ratio defined as the amplitude of the pressure oscillations induced by a shear strain normalized to the volumetric strain. Three angles and two values are needed to calculate the azimuth dependency of the pore pressure response: the angle between the directions of the maximum horizontal stress and the seismic event; fault zone orientation; micro-crack orientation within the fault zone; damage and stress values. The model predicts that maximum pore pressure response occurs when micro-cracks and maximum horizontal stress are in the same orientation, high damage, and high stress anisotropy. By adjusting these quantities, we recalculate results of recent seismological studies in Arbuckle disposal well, Osage County, Oklahoma. The presented model successfully predicts the observed azimuthal dependence in wave-induced fluid pressure response, and relates the anisotropic response to tectonic indicators such as the orientations of the maximum horizontal stress, fault zone, and micro-fractures.

SummaryIn near-surface surveys, shallow-seismic and ground-penetrating radar (GPR) full-waveform inversions (FWIs) have received increasing attention because of their ability to reconstruct high-resolution subsurface models. However, they have different sensitivities to the same targets and thus may yield conflicting geophysical parameter models. To solve this issue, we have developed an indirect joint petrophysical inversion (JPI) integrating shallow-seismic and multi-offset GPR data. These data are used to reconstruct porosity and saturation whereby we use only strong sensitivities between petrophysical and geophysical parameters. To promote its field application, we proposed an input strategy to avoid measuring rock matrix parameters and make indirect JPI more robust. We apply indirect JPI to the field data acquired in Rheinstetten, Germany and find that it reveals the mechanical, electrical and petrophysical properties more reliably than individual inversions. The reconstructed models are assessed by direct push technology, borehole sample measurements and migrated GPR image. Indirect JPI can fit seismic and GPR observed data simultaneously and provide consistent multi-parameter models, which are hard to achieve by FWIs and individual petrophysical inversions. We also find that the method is robust when there are uncertainties in petrophysical a priori information. Overall, the field example proves the great potential of using indirect JPI to solve real-world problems.

SummaryEquations of motion are derived for (visco)elastic, self-gravitating, variably-rotating planets. The equations are written using a decomposition of the elastic motion that separates the body’s elastic deformation from its net translational and rotational motion as far as possible. This separation is achieved by introducing degrees of freedom that represent the body’s rigid motions; it is made precise by imposing constraints that are physically motivated and that should be practically useful. In essence, a Tisserand frame is introduced exactly into the equations of solid mechanics. The necessary concepts are first introduced in the context of a solid body, motivated by symmetries and conservation laws, and the corresponding equations of motion are derived. Next, it is shown how those ideas and equations of motion can readily be extended to describe a layered fluid–solid body. A possibly new conservation law concerning inviscid fluids is then stated. The equilibria and linearisation of the fluid–solid equations of motion are discussed thereafter, along with new equations for use within normal-mode coupling calculations and other Galerkin methods. Finally, the extension of these ideas to the description of multiple, interacting fluid–solid planets is qualitatively discussed.

SummaryFaults in general, and in clay materials in particular, have complex structures that can be linked to both a polyphased tectonic history and to the anisotropic nature of the intact rock. Drilling through faults in shaly materials allows measuring properties such as the structure, mineralogical composition, stress orientation and physical properties. We combine different petrophysical measurements on core samples retrieved from a borehole drilled perpendicularly to a fault zone affecting Toarcian shales from the Tournemire underground research laboratory (France). The borehole is cross-cutting the entire fault thickness which is of the order of ten meters. We perform several types of measurements: density, porosity, saturation directly in the field, and P-wave velocities together with P-waves anisotropy on core samples taken at regular intervals. Special protocols were developed to preserve as much as possible the saturation state of the samples. From our measurements, we were able to track the increase of damage, characterized by a smooth decrease in elastic moduli from the intact zone to the fault core. We then calculated Thomsen's parameters to quantify the elastic anisotropy evolution across the fault. Our results show strong variations of the elastic anisotropy with the distance to the fault core as well as the occurrence of anisotropy reversal.

SummaryBorneo and Sulawesi are two large islands separated by the Makassar Strait that lie within the complex tectonic setting of central Indonesia. The seismic structure beneath this region is poorly understood due to the limited data availability. In this study, we present Rayleigh wave tomography results that illuminate the underlying crustal structure. Group velocity is retrieved from dispersion analysis of Rayleigh waves extracted from the ambient noise field by cross-correlating long-term recordings from 108 seismic stations over a period of 8 months. We then produce a 3-D shear wave velocity model via a two-stage process in which group velocity maps are computed across a range of periods and then sampled over a dense grid of points to produce pseudo-dispersion curves; these dispersion curves are then separately inverted for 1-D shear wave velocity (Vs), with the resultant models combined and interpolated to form a 3-D model. In this model, we observed up to ± 1.2 km/s lateral Vs heterogeneities as a function of depth. Our models illuminate a strong low shear wave velocity (Vs) anomaly at shallow depth (≤ 14 km) and a strong high Vs anomaly at depths of 20 – 30 km beneath the North Makassar Strait. We inferred the sediment basement and Moho depth from our 3-D Vs model based on iso-velocity constrained by the positive vertical gradient of the Vs models. The broad and deep sedimentary basement at ∼14 ± 2 km depth beneath the North Makassar Strait is floored by a shallow Moho at ∼22 ± 2 km depth, which is the thinnest crust in the study area. To the east of this region, our model reveals a Moho depth of ∼45 ± 2 km beneath Central Sulawesi, the thickest crust in our study area, which suggests crustal thickening since the late Oligocene. Moreover, the presence of high near surface Vs anomalies with only slight changes of velocity with increasing depth in southwest Borneo close to Schwaner Mountain (SM) confirm the existence of a crustal root beneath this region.

SummaryIn this work, we study the induced seismicity recorded during an injection operation at the Muara Laboh geothermal plant (Indonesia). The swarm, consisting of three bursts activating a normal fault zone, is characterized by rapid earthquake (km/day) migration. We use a 2D rate-and-state asperity model to better understand the physical mechanisms controlling the evolution of this induced swarm. The model suggests that the observed rapid seismic migration can be explained by the interaction among asperities through the expansion of slow postseismic slip fronts. Also, it shows that the amount of seismicity generated by the fluid injection is strongly controlled by the background seismicity of the system, that is by the seismicity determined by the tectonic load charging the fault. This close correlation between natural and induced seismicity suggests that the injection in Muara Laboh principally stimulates critically stressed faults, which release the seismicity determined by their natural seismic cycle.

SummaryThe Mw 6.8 Murghob earthquake is the third earthquake in an Mw 6.4+ sequence occurring in the Pamir initiated by the 2015 Sarez Mw 7.2 earthquake. It is of great significance to investigate their interactions and to assess future seismic hazards in the region. In this paper, we use Sentinel-1 radar interferometric data to retrieve coseismic deformation, invert for the slip distributions of the four events, and then investigate their interactions. The cumulative Coulomb failure stress changes (ΔCFS) suggest that the 2023 Murghob earthquake was promoted by the three prior earthquakes in the sequence. Pre-stress from historical earthquakes is a key factor in explaining the triggering mechanism of the two 2016 Mw 6.4+ earthquakes. Stress loading and unloading effects on major faults in the region indicate that future attention should be paid in 1) the segment of the Sarez-Karakul fault north of the Kokuibel Valley, 2) the segment of the Sarez-Murghab thrust fault west of the Sarez-Karakul fault, and 3) the east segments of the Pamir thrust fault system, all with a large positive ΔCFS.

SummaryThe Lijiang–Xiaojinhe fault (LXF) and its vicinity are located in the transition zone among the Tibetan Plateau (TP), the South China block and the Indochina block. Researchers believe that this area has acted as a key tectonic zone during the evolution of the TP. Owing to the continuous growth and SE-ward expansion of the TP, the LXF and its vicinity have experienced intense deformation. Although different models, such as the rigid block extrusion and mid-lower crustal flow models, have been proposed to explain this intense deformation, a consensus has not yet been achieved. To better understand the deformation of the LXF and its vicinity, a high-resolution image of the subsurface structure must be constructed. In this study, we construct images of wave velocity and azimuthal anisotropy structures by using an eikonal equation-based traveltime tomography method. We collect high-quality seismic data from 276 broadband seismic stations and manually pick a total of 48,037 first arrivals for the tomography study. Our tomographic results reveal a strong low-velocity body below the LXF and its vicinity. In addition, a strong azimuthal anisotropy structure with a N–S-oriented fast velocity direction (FVDs) is distributed along the low-velocity body. These features indicate the occurrence of mid-lower crustal flow, that penetrates across the LXF and extends to the Dianzhong block (DZB). In addition, we find obvious low-velocity perturbations in the mid-lower crust and uppermost mantle beneath the DZB. The low-velocities may be attributed to the upwelling of hot materials from the upper mantle. We consider the limited distribution of mid-lower crustal flow on the margin of the SE TP, and mid-lower crustal flow may not play a significant role in the expansion of the TP.

SummaryThe recent rapid improvement of machine learning techniques had a large impact on the way seismological data can be processed. During the last years several machine learning algorithms determining seismic onset times have been published facilitating the automatic picking of large data sets. Here we apply the deep neural network PhaseNet to a network of over 900 permanent and temporal broad band stations that were deployed as part of the AlpArray research initiative in the Greater Alpine Region (GAR) during 2016-2020. We selected 384 well distributed earthquakes with ML ≥ 2.5 for our study and developed a purely data-driven pre-inversion pick selection method to consistently remove outliers from the automatic pick catalog. This allows us to include observations throughout the crustal triplication zone resulting in 39,599 P and 13,188 S observations. Using the established VELEST and the recently developed McMC codes we invert for the 1D P- and S-wave velocity structure including station correction terms while simultaneously relocating the events. As a result we present two separate models differing in the maximum included observation distance and therefore their suggested usage. The model AlpsLocPS is based on arrivals from ≤ 130 km and therefore should be used to consistently (re)-locate seismicity based on P & S observations. The model GAR1D_PS includes the entire observable distance range of up to 1000 km and for the first time provides consistent P- & S-phase synthetic travel times for the entire Alpine orogen. Comparing our relocated seismicity with hypocentral parameters from other studies in the area we quantify the absolute horizontal and vertical accuracy of event locations as ≈ 2.0 km and ≈ 6.0 km, respectively.

AbstractCombining elastic full waveform inversion (FWI) with rock physics holds promise for expanding the application of FWI beyond seismic imaging to predicting and monitoring reservoir properties. Distributed Acoustic Sensing (DAS), a rapidly developing seismic acquisition technology, is being explored for its potential in supporting FWI applications. In this study, we implement a sequential inversion scheme that integrates elastic FWI and Bayesian rock physics inversion, using a vertical seismic profile (VSP) dataset acquired with accelerometer and collocated DAS fiber at the Carbon Management Canada’s Newell County Facility. Our aim is to establish a baseline model of porosity and lithology parameters to support later monitoring of CO2 storage. Key strategies include an effective source approach for addressing near-surface complications, a modeling strategy to simulate DAS data comparable to field data, and a Gaussian mixture approach to capture the bimodality of rock properties. We conduct FWI tests on accelerometer, DAS, and combined accelerometer-DAS data. While our inversion results accurately reproduce either dataset, the elastic models inverted from accelerometer data outperform the other two in matching well logs and identifying the target reservoir. We attribute this outcome to the limited complementarity of DAS data with accelerometer data in our experiment, along with the limitations imposed by single-component measurements on DAS. The porosity and lithology models predicted from accelerometer-derived elastic models are reasonably accurate at the well location and exhibit geologically meaningful spatial distribution.

SummaryFractured hydrate-bearing reservoirs are extensively discovered worldwide and show notable anisotropic geophysical properties. Hydrate distribution in fractures significantly affects the anisotropic properties, and hence plays an important role in the accurate assessment of hydrate resources. However, the knowledge about how the hydrate distribution impacts the anisotropic geophysical properties of fractured reservoirs, which is the premise for the identification and quantification of hydrate in fractured reservoirs, is still poorly understood. To obtain such knowledge, we forward study the effects of various hydrate distribution (including floating, bridging and evolving hydrate distribution) in aligned fractures on the anisotropic elastic, electrical and joint elastic-electrical properties of a digital core using validated numerical methods. We show that for all the hydrate distribution, the anisotropic velocities increase while the conductivities decrease with increasing hydrate saturation, with the effects of the floating and bridging distribution being the least and greatest, respectively. We also show that the anisotropic velocities and conductivities for the floating and bridging distribution vary approximately linearly with hydrate saturation, leading to linear correlations between the elastic and electrical properties. Further investigation illustrates that the difference in the slopes of the linear joint correlations between the two distribution is significantly greater than that of the individual properties as a function of hydrate saturation. The results have revealed the distinct effects of hydrate distribution on the anisotropic elastic and electrical properties of fractured reservoirs, and have confirmed the superiority of the joint elastic-electrical properties for the distinguishment of hydrate distribution in fractures over individual physical properties.

SummaryPrevious studies have shown that it is difficult to determine whether the 2015 Pishan earthquake occurred on a uniform fault or a ramp-flat fault with variable dip angles due to the similar goodness of data fit to co-seismic and afterslip models on these two fault models. Here, we first present the InSAR deformation obtained from both ascending and descending orbits, covering the coseismic period and cumulative five-year period after the 2015 Pishan earthquake. We then determine the preferred fault geometry by the spatial distributions between the positive Coulomb failure stress change triggered by mainshock and the afterslip. Based on the preferred fault model, we finally use a combined model to determine the contributions of elastic and viscoelastic deformation in the post-seismic deformation. We find that the Pishan earthquake prefers to occur on a ramp-flat fault, and the coseismic slip is mainly distributed at a depth of 9-13 km, with a maximum slip of about 1.3 m. The post-seismic deformation is primarily governed by afterslip, as the poroelastic rebound-induced deformation fails to account for the observed post-seismic deformation and the contributions from the viscoelastic relaxation mechanism can be considered negligible in the combined model. Moreover, the modelled stress-driven afterslip and observed kinematic afterslip have good consistency, and the difference between the root mean square error (RMSE) of the two afterslip models is only 4.3 mm. The results from the afterslip model indicate that both of the up- and down-dip directions distribute the afterslip, and slip in the up-dip direction is greater than that of the down-dip direction. Meanwhile, the maximum cumulative afterslip after five years is approximately 0.26 m which is equivalent to a released seismic moment of a Mw 6.47.

SummaryThe Earth’s magnetic field is generated by a dynamo in the outer core and is crucial for shielding our planet from harmful radiation. Despite the established importance of the core-mantle boundary heat flux as driver for the dynamo, open questions remain about how heat flux heterogeneities affect the magnetic field. Here, we explore the distribution of core-mantle boundary heat flux on Earth and its changes over time using compressible global 3-D mantle convection models in the geodynamic modeling software ASPECT. We discuss the use of the consistent boundary flux method as a tool to more accurately compute boundary heat fluxes in finite element simulations and the workflow to provide the computed heat flux patterns as boundary conditions in geodynamo simulations. Our models use a plate reconstruction throughout the last 1 billion years—encompassing the complete supercontinent cycle—to determine the location and sinking speed of subducted plates. The results show how mantle upwellings and downwellings create localized heat flux anomalies at the core-mantle boundary that can vary drastically over Earth’s history and depend on the properties and evolution of the lowermost mantle as well as the surface subduction zone configuration. The distribution of hot and cold structures at the core-mantle boundary changes throughout the supercontinent cycle in terms of location, shape and number, indicating that these structures fluctuate and might have looked very differently in Earth’s past. We estimate the resulting amplitude of spatial heat flux variations, expressed by the ratio of peak-to-peak amplitude to average heat flux, q*, to be at least 2. However, depending on the material properties and the adiabatic heat flux out of the core, q* can easily reach values >30. For a given set of material properties, q* generally varies by 30-50% over time. Our results have implications for understanding the Earth’s thermal evolution and the stability of its magnetic field over geological timescales. They provide insights into the potential effects of the mantle on the magnetic field and pave the way for further exploring questions about the nucleation of the inner core and the past state of the lowermost mantle.

SummaryThe wide-swath altimeter satellite Surface Water and Ocean Topography (SWOT) will provide high spatiotemporal resolution sea surface heights (SSHs), which is crucial for studying the impact of observation errors on marine gravity recovery. This study uses simulated SWOT data to derive deflection of the vertical (DOV) and gravity anomalies in the northern South China Sea. We quantified the impact of SWOT errors on DOV and gravity anomalies, and analyzed the contributions from different directions of geoid gradient. The results show that the geoid gradient in the cross-track direction significantly improves gravity field recovery by enhancing the precision of east component of DOV. For one-cycle SWOT observations, phase errors emerge as the most impactful error affecting both DOV and gravity anomalies, followed by random errors. Two-dimensional Gaussian filtering and the tilt correction proposed in this study could effectively mitigate their impact. Using the corrected data for DOV computation, the precision in the east and north components improves by 75.32% and 46.80%, respectively, while enhancing the accuracy of the gravity field by 70.23%. For 17-cycle data, phase errors and random errors remain the predominant factors affecting DOV and gravity anomalies, but their impact diminishes with an increase in SWOT observations. Our results indicate that marine gravity accuracy improves by approximately 70% compared to a single cycle. Whether for single-cycle or multi-cycle data, the impact of phase errors is roughly twice that of random errors. These data processing strategies can serve as valuable references for wide-swath altimeter data processing, aiming to advance the precision and resolution of marine gravity field recovery.

SummaryAround the world the Earth’s crust is blanketed to various extents by sediment. For continental regions, knowledge of the distribution and thickness of sediments is crucial for a wide range of applications including seismic hazard, resource potential, and our ability to constrain the deeper crustal geology. Excellent constraints on the sediment thickness can be obtained from borehole drilling or active seismic surveys. However, these approaches are expensive and impractical in remote continental interiors such as central Australia. Recently, a method for estimating the sediment thickness using passive seismic data, the collection of which is relatively simple and low-cost, was developed and applied to seismic stations in South Australia. This method uses receiver functions, specifically the time delay of the P-to-S converted phase generated at the sediment-basement interface, relative to the direct-P arrival, to generate a first order estimate of the thickness of sediments. In this work we expand the analysis to the vast array of over 1500 seismic stations across Australia, covering an entire continent and numerous sedimentary basins that span the entire range from Precambrian to present-day. We compare with an established yet separate method to estimate the sediment thickness, which utilises the autocorrelation of the radial receiver functions to ascertain the two-way travel-time of shear waves reverberating in a sedimentary layer. Across the Australian continent the new results match the broad pattern of expected sedimentary features based on the various geological provinces. We are able to delineate the boundaries of many sedimentary basins, such as the Eucla and Murray Basins, which are Cenozoic, and the boundary between the Karumba Basin and the mineral rich Mount Isa Province. Contrasts in seismic delay time across these boundaries are upwards of 0.4 s. The delay signal is found to diminish to <0.1 s for older Proterozoic basins, likely due to compaction and metamorphism of the sediments over time. As an application of the method, a comparison with measurements of sediment thickness from local boreholes allows for a straightforward predictive relationship between the delay time and the cover thickness to be defined. This offers future widespread potential, providing a simple and cheap way to characterise the sediment thickness in under-explored areas from passive seismic data.

SummaryEarthquake forecasting poses significant challenges, especially due to the elusive nature of stress states in fault systems. To tackle this problem, we employ features derived from seismic catalogs obtained from acoustic emission (AE) signals recorded during triaxial stick-slip experiments on natural fractures in three Westerly granite samples. We extracted 47 physically explainable features from AE data that described spatio-temporal evolution of stress and damage in the vicinity of the fault surface. These features are then subjected to unsupervised clustering using the K-means method, revealing three distinct stages with a proper agreement with the temporal evolution of stress. The recovered stages correspond to the mechanical behavior of the rock, characterized as initial stable (elastic) deformation, followed by a transitional stage leading to an unstable deformation prior to failure. Notably, AE rate, clustering-localization features, fractal dimension, b-value, interevent time distribution, and correlation integral are identified as significant features for the unsupervised clustering. The systematically evolving stages can provide valuable insights for characterizing preparatory processes preceding earthquake events associated with geothermal activities and waste-water injections. In order to address the upscaling issue, we propose to use the most important features and, in case of normalization challenge, removing non-universal features, such as AE rate. Our findings hold promise for advancing earthquake prediction methodologies based on laboratory experiments and catalog-driven features.

SummarySubduction earthquakes show complex spatial and temporal rupture patterns, exhibiting events of varied sizes, which rupture distinct or overlapping fault segments. Elucidating first-order controlling conditions of rupture segmentation and return periods of large earthquakes is therefore critical for seismic and tsunami hazard estimations. The Chilean subduction zone frequently hosts several Mw > 8 earthquakes, with heterogeneous recurrence rates and locations. Here, we implement 3D quasi-dynamic rate and state frictional models to investigate the role of plate interface geometry on the distribution of interseismic coupling and coseismic ruptures in Central Chile. First, we develop synthetic-parametric models that show how dip and strike variations may increase the probabilities to produce partial seismic barriers, which tend to avoid the production of large earthquake ruptures and modulate rupture lengths. Then, we simulate the subduction seismic cycle processes on Central Chile (25ºS-38°S), imposing depth-dependent frictional properties on a realistic non-planar 3D subduction interface geometry. Similar to results obtained for synthetic-parametric models, after 5000 years of simulation, regions with abrupt dip or strike changes increase the probabilities of stopping coseismic propagation of simulated Mw8.0-9.0 earthquakes. Our simulated earthquake sequences on the Central Chile subduction zone delimit rupture areas that match geometrical interface features and historical earthquakes, results that point to the crucial role of fault interface geometry on seismic cycle segmentation along strike.