Geophysical Journal International

SummarySeismic source models that use an elastic relation between pressure decrease, compaction, and stress change have been shown to successfully reproduce induced seismicity in producing natural gas reservoirs undergoing differential compaction. However, this elastic relation is inconsistent with observations of non-linear reservoir compaction in the Groningen field. We utilize critical state mechanics theory to derive a 3D stress-strain framework that is able to house 1D non-linear stress-strain relations typically used for subsidence models, without the need for recalibration of the subsidence model parameters. This is used to adapt the elastic thin sheet stress model that is currently in use as the state-of-the-art for seismicity predictions as part of the hazard and risk assessment of the Groningen gas field. The new thin sheet model has one additional model parameter that modulates the impact of inelastic deformation on fault loading, whilst keeping the intended function of the model calibration from the original elastic thin sheet model intact. The resulting elastic-viscoplastic thin sheet stress model is consistent with previously reported non-linear rate-dependent reservoir compaction in Groningen found from inverting subsidence data and from rock deformation experiments. Our elastic-viscoplastic thin sheet stress model is able to predict ongoing stress increase, and therefore ongoing seismicity, in areas where pressure does not decrease anymore due to shut-in. A pseudo-prospective forecasting exercise indeed shows that the elastic-viscoplastic stress model performs better than the linear elastic stress model. This model addition ensures that the Groningen seismic source model is well suited for predicting seismicity in the post shut-in phase.

SummaryThis paper presents a novel framework for full-waveform seismic source inversion using simulation-based inference (SBI). Traditional probabilistic approaches often rely on simplifying assumptions about data errors, which we show can lead to inaccurate uncertainty quantification. SBI addresses this limitation by learning an empirical probabilistic relationship between the parameters and data, without making assumptions about the data errors. This is achieved through the use of specialised machine learning models, known as neural density estimators, which can then be integrated into the Bayesian inference framework. We apply the SBI framework to point-source moment tensor inversions as well as joint moment tensor and time-location inversions. We construct a range of synthetic examples to explore the quality of the SBI solutions, as well as to compare the SBI results with standard Gaussian likelihood-based Bayesian inversions. We then demonstrate that under real seismic noise, common Gaussian likelihood assumptions for treating full-waveform data yield overconfident posterior distributions that underestimate the moment tensor component uncertainties by up to a factor of 3. We contrast this with SBI, which produces well-calibrated posteriors that generally agree with the true seismic source parameters, and offers an order-of-magnitude reduction in the number of simulations required to perform inference compared to standard Monte Carlo sampling techniques. Finally, we apply our methodology to a pair of moderate magnitude earthquakes in the North Atlantic. We utilise seismic waveforms recorded by the recent UPFLOW ocean bottom seismometer array as well as by regional land stations in the Azores, comparing full moment tensor and source-time location posteriors between SBI and a Gaussian likelihood approach. We find that our adaptation of SBI can be directly applied to real earthquake sources to efficiently produce high quality posterior distributions that significantly improve upon Gaussian likelihood approaches.

SummaryThe thermal structure of the continental crust plays a critical role in understanding its elastic and rheologic properties as well as its dynamic processes. Thermal parameter datasets on continental scales have been used to constrain the crustal thermal structure, including both the direct (e.g., temperature, heat flux, and heat conductivity measured at the surface) and indirect (e.g., seismically derived Mohorovičić discontinuity (Moho) temperature, geomagnetically derived Curie depth) observations. In this study, we present a new continental scale crustal heat generation model with additional information from seismologically-inferred crustal composition. Together with previous direct and indirect thermal parameter datasets in the conterminous United States, we use the new crustal heat generation model to construct a 3-dimensional (3-D) crustal temperature model under a newly developed Bayesian framework. Specifically, we first derive profiles of crustal heat generation based on an empirical geochemical relationship at 1683 locations where seismologically derived crustal composition information is available. Then for each of these locations, the average heat generation values in the upper, middle, and lower crust are combined with other thermal parameters through a Markov Chain Monte-Carlo inversion for a conductive, vertically smooth temperature profile. The results, posterior distributions of temperature profiles, are used to generate a 3-D crustal thermal model with the uncertainties systematically assessed. The new temperature model overall exhibits similar patterns to that from the U.S. Geological Survey National Crustal Model, but also reduces possible biases and the model's dependence on a single thermal parameter.

SummaryThe ∼300km-long rupture of the February 6 2023 Kahramanmaraş earthquake began in the Narlı section of the Karasu trough, a pull-apart basin sandwiched and sheared between the two major strike-slip faults of the region, the East Anatolian Fault (EAF) on the west and the Dead Sea Fault (DSF) on the east. Rupture started where the northern segment of the DSF enters the Narlı Basin with Mw7.0 sub-event and propagated across the basin before making its junction with the EAF. In the seven months preceding the earthquake this basin was the seat of anomalous seismic activity. This activity occurred in bursts interweaved with periods of quiescence. It started near-concomitantly in two clusters located on the opposite edges of the pull-apart basin ∼20 km apart. The organization of this seismicity into families of numerous repeating earthquakes suggests an aseismic process linked to fault healing and rapid reloading in a critically stressed zone. By December 2022, two months before the earthquake, activity had migrated to a cluster located along the path that rupture was to follow during the initial stage of the earthquake. These observations show that the pull-apart basin where rupture started was progressively deforming in a succession of bursts before the earthquake. This regional-scale deformation is closely linked with the transitional nature of geodynamics and kinematics influenced by large-scale fault interactions in the surrounding area. The location of the epicenter near the northern termination of the rupture of the 1822 M7.4 earthquake suggests that the ∼45 km long Narlı sub-rupture which constituted the first stage of this giant earthquake was closing a long-present seismic gap between the DSF and the EAF.

SummaryOcean-bottom seismic acquisitions are gaining widespread popularity across a variety of subsurface applications. However, the high cost of these systems often necessitates receiver geometries with large intervals between ocean-bottom cables or nodes. The upside-down Rayleigh-Marchenko (UD-RM) method has been recently proposed as an effective solution for accurate redatuming and imaging of sparse seabed data. In this paper, we present the first successful application of the UD-RM method to field data. We demonstrate that in the presence of a shallow seabed, an improved data pre-processing workflow is necessary to generate more accurate input wavefields compared to the one produced by the workflow presented in the original paper. To validate the proposed processing workflow, the UD-RM method is initially tested on a synthetic dataset that mimics the Volve field data (referred to as the Volve synthetic dataset); this is followed by its application to a 2D line of the Volve ocean-bottom cable dataset. Subsequently, the field dataset is subsampled by retaining only 25 percent of the total receivers to numerically validate the UD-RM method’s capability to handle sparse receiver arrays. The resulting images reveal that the UD-RM method, when paired with our enhanced data processing workflow, can effectively handle surface-related multiples, internal multiples, and sparse receiver arrays, producing accurate imaging results without the need for costly and labor-intensive multiple removal processes. Finally, we provide theoretical insights and numerical evidence supporting the necessity of source-side deghosting prior to redatuming. While a pre-processing workflow that omits source-side deghosting can offer some practical advantages, we show that this ultimately produces blurrier images compared to those obtained using source-side deghosted input data.

SummaryFull waveform inversion (FWI) has the potential to provide high-resolution insights into subsurface structures. However, its adoption, particularly in 3D multiparameter applications, has been limited by high computational costs. This study addresses this challenge by introducing an optimized experimental design (OED) method that simultaneously optimizes source placement and model parameterization. The result is an optimized survey design and a compressed model representation that maximizes the information content. By reducing the source layout by approximately 50% and compressing the model by approximately 90%, this approach significantly reduces computational demands, allowing the use of fast convergence inversion algorithms such as the Gauss-Newton method. The OED calculation is reduced from a typical $\mathcal {O}(n^3)$ complexity, as in eigenvalue-based criteria, to $\mathcal {O}(n \log _2n)$ with the newly introduced wavelet transform-based criterion. Additionally, a post-acquisition source-receiver pair optimization method is developed, demonstrating that while random selection captures high information content, the proposed OED criterion effectively minimizes the number of required simulations. This approach further reduces computational cost and facilitates the efficient extraction of compact, high-value datasets from excessively large surveys.

SummaryThe Ecuadorian Andes are a complex region characterized by accreted oceanic terranes driven by the ongoing subduction of the oceanic Nazca plate beneath South America. Present-day tectonics in Ecuador are linked to the down-going plate geometry featuring the subduction of the aseismic, oceanic Carnegie Ridge, which is currently entering the trench. Using seismic tomography, we jointly invert arrival times of P and S waves from local and teleseismic earthquakes with surface wave dispersion curves to image the structure of the forearc and magmatic arc of the Ecuadorian Andes. Our dataset includes >100,000 travel-times recorded at 294 stations across Ecuador. Our images show the basement of the central forearc is composed of accreted oceanic terranes with high elastic wavespeeds. Inboard of the Carnegie Ridge, the westernmost forearc and coastal cordilleras display relatively low Vp and Vs and high Vp/Vs values, which we attribute to the increased hydration and fracturing of the overriding plate due to the subduction of the thick oceanic crust of the Carnegie Ridge. We additionally image across-arc differences in magmatic architecture. The frontal volcanic arc overlies accreted terranes and is characterized by low velocities and high Vp/Vs indicative of partial melt reservoirs which are limited to the upper crust. In contrast, the main arc displays regions of partial melt across a wider range of depths. The Subandean zone of Ecuador has two active volcanoes built on continental crust suggesting the arc is expanding eastward. The mid to lower crust does not show indications of being modified from the magmatic process. We infer that the slab is in the process of flattening as a consequence of early-stage subduction of the buoyant Carnegie Ridge.

SummaryLocating the sources of observed disturbances in potential-field data is a challenging problem due to the non-unique nature of the inverse problem. The Euler deconvolution method was created to solve this issue, particularly for idealized sources (such as spheres and planar vertical dykes). Euler deconvolution has become widely used in potential-field methods due, in large part, to its low computational cost and ease of implementation into software. However, it is widely known that Euler deconvolution suffers from some shortcomings: 1) non-uniqueness of the solution with respect to the depth of the source and the structural index (a parameter that represents the idealised shape of the source); 2) sensitivity to short-wavelength noise in the data derivatives which are used as inputs for the method. Here, we present a new method called Euler inversion which is a reformulation of the inverse problem of Euler’s homogeneity equation as an implicit mathematical model rather than a parametric one. Euler inversion is a constrained, non-linear inverse problem capable of estimating both the model parameters (location of the source and constant base level) and the predicted data (potential field and its derivatives). We show that Euler inversion is less sensitive than Euler deconvolution to short-wavelength noise and to the presence of interfering sources in the data window. By also estimating the predicted data, Euler inversion is also able to estimate the best integer structural index to be used for inversion. Our results show that the estimated structural index minimizes the data misfit and coincides with those of the simulated sources. Furthermore, most matrices involved in the method are either sparse or diagonal, making Euler inversion computationally efficient. Tests on synthetic data and a real aeromagnetic dataset from Rio de Janeiro, Brazil, demonstrate the effectiveness of Euler inversion to delineate sources with variable geometries and correctly estimate their depths.

SummarySolving the wave equation is essential to seismic imaging and inversion. The numerical solution of the Helmholtz equation, fundamental to this process, often encounters significant computational and memory challenges. We propose an innovative frequency-domain scattered wavefield modeling method employing neural operators adaptable to diverse seismic velocities. The source location and frequency information are embedded within the input background wavefield, enhancing the neural operator’s ability to process source configurations effectively. In addition, we utilize a single reference frequency strategy, which enables scaling from larger-domain forward modeling to higher-frequency scenarios, thereby improving our method’s accuracy and generalization capabilities for larger-domain applications. Several tests on the OpenFWI datasets and realistic velocity models validate the accuracy and efficacy of our method as a surrogate model, demonstrating its potential to address the computational and memory limitations of numerical methods.

SummaryThe response of well water temperature to earthquakes is crucial for understanding subsurface seismic fluid dynamics. However, recent studies have primarily focused on observations at a single depth and have employed single-aquifer models, which may lead to controversies when explaining fluid flow. In this study, we develop single- and double-aquifer models to estimate well-water temperature variations at different depths in response to changes in pore pressure, permeability, and aquifer recharge temperature. The results indicate that variations in aquifer pore pressure and permeability result in significant differences in vertical flow velocity and temperature changes at various depths. When the borehole bottom is impermeable, for a single aquifer, temperature variation is maximal above the aquifer and variable at the aquifer depth, but nearly zero below the aquifer; for two aquifers, different pore pressure and permeability changes in each aquifer produce distinct temperature variation patterns, with minimal temperature change below the lower aquifer. If the borehole bottom is permeable, temperature variation becomes obvious below the lower aquifer. When cold or hot water from the aquifers flows into the borehole, significant temperature perturbations remain confined within a few metres of the aquifer within one day. Finally, a field case study investigates the co-seismic water temperature responses at three depths in the Chuan No. 03 well, triggered by the 2008 Mw 7.9 Wenchuan earthquake. The double-aquifer model effectively explains the complex co-seismic temperature fluctuations at different depths. Observation at a single depth risk missing crucial information, and multi-depth temperature observation is a promising approach for interpreting and monitoring groundwater responses to earthquakes.

SummaryExtreme value statistics (EVS) is commonly used to model rare, extreme events such as natural disasters. This study proposes a method that integrates EVS and Bayesian estimation to enable the early forecasting of aftershock-induced ground shaking. The method is applied to continuous seismograms recorded immediately after a large earthquake. The proposed method is based on several key assumptions: the Gutenberg–Richter (G–R) and Omori–Utsu laws, as well as the proportionality between earthquake magnitude and the logarithmic maximum amplitude. Based on these assumptions, two metrics were computed at each seismic station: the exceedance probability of the maximum amplitude (EPMA) and the number of threshold value exceedances (EPNUM). While EPMA follows a long-tailed Fréchet distribution, with uncertainty spanning at least an order of magnitude, EPNUM follows a short-tailed Poisson distribution, with uncertainty typically varying by a factor of two. The performance of the proposed method was evaluated across three different types of aftershock sequences in Japan. The practical forecasting capability was demonstrated within 1 hour of the mainshock and was effective up to 7 days. Compared to conventional methods that rely on incomplete earthquake catalogues, the proposed approach demonstrated faster and more robust results. While the median forecast of maximum amplitude tended to be overestimated, possibly due to the potential nonlinear relationship between magnitude and logarithmic maximum amplitude, the forecast for the number of felt earthquakes did not show such bias. Because the proposed method is based on single-station processing, it can be applied in regions without a dense seismograph network or real-time earthquake monitoring system, as long as continuous ground motion data is available at the target site.

SummaryThe interpretation of subsurface resistivity structures in volcanic areas remains challenging and requires the selection of the most plausible configuration from various geological features that affect the resistivity measurements. A comprehensive physical study of the rocks in the target area is essential for accurate interpretation. A magnetotelluric (MT) survey conducted in the southeast flank of Mt. Ontake volcano detected several kilometres of low resistivity below 30 Ωm. However, the interpretation of resistivity anomalies remains to be verified; borehole data have been used to resolve this problem. In this study, geological, fracture, temperature, and resistivity structures obtained from borehole investigations were analysed using rock physical methods, and core samples were subjected to physical property measurements and microscopic observations. The geology observed in the borehole was sedimentary rock with a porosity of less than a few percent, except for the surface volcanic breccia and some intrusive granites. The borehole wall was poorly fractured. The groundwater temperature in the borehole was aligned with the standard geothermal gradient, 40 °C at a depth of 800 m. These results indicate the absence of a hydrothermal reservoir at the site. In contrast, the core samples contained pyrite-filled microfractures. All the samples had poor porosity and contained a small amount of clay minerals but showed uniformly lower resistivity values than expected for similar porosity. Some samples exhibited robust induced polarization effects. The pyrite contents of the samples were low. A pseudo-high-frequency conductivity model test using electrostatic field analysis on a numerical model of microfractured rock properties reproduced the low resistivity observed in the borehole investigation and MT survey. The specifically low resistivity values observed were appropriately reproduced. The findings of this study indicate that pyrite filling the microfractures can be a major contributor to the low resistivity below 1 Ωm in the study area and suggest that it is necessary to consider the network of microscale veins formed by conductive sulfide minerals, mainly pyrite, for the interpretation of subsurface resistivity in geothermal areas.

SummaryThe intricate architecture of plant root systems is crucial for nutrient and water uptake, significantly influencing plant growth and productivity. Induced polarization (IP) is a promising non-destructive technique for analyzing plant roots in their natural conditions. This study introduces a novel theoretical and numerical model to explain the significant low-frequency polarization of plant root cells observed in previous experiments. Our approach addresses the limitations of existing models by incorporating geometric constraints and internal mechanisms of cell polarization, particularly focusing on interfacial polarization across the cell membrane. Through comprehensive simulations, we investigate various geometries and boundary conditions, demonstrating that densely packed root cells exhibit significant polarization signals within a measurable frequency range due to coupling effects. Our findings align with experimental observations, indicating that the peak frequency is highly sensitive to cell arrangement and membrane properties, while the maximum phase shift remains consistent. This model provides a robust framework for interpreting polarization signals in root systems, offering potential applications for in-situ characterization of plant roots and enhancing the understanding of root dynamics under different environmental conditions.

SummaryIn the previous papers of this series, we have developed an in-depth analysis of the low-frequency complex conductivity response of volcanic (extrusive) rocks. We showed that the alteration of these rocks plays a key-role in determining their induced polarization properties, especially regarding the formation of smectite in response to the thermo-activated alteration of the volcanic glasses. We also considered the effects associated with the presence of magnetite and pyrite. In the present paper, we look at the induced polarization properties of igneous rocks like granites and granitoids. Usually, the alteration path of these rocks leads to the formation of kaolinite, a clay mineral with a much lower Cation Exchange Capacity (CEC) than smectite. Thirty-three core samples from 3 sites in France are saturated with NaCl solutions at 3 salinities (pore water conductivity of 0.1, 1.0, and 10 S m−1, 25°C) and their complex conductivity spectra are measured in the frequency range 0.01 Hz-45 kHz. As observed for volcanic rocks, the surface conductivity, normalized chargeability, and quadrature conductivity depend strongly on the CEC of the rock, which is independently measured with the cobalt-hexamine method. The (intrinsic) formation factor follows an Archie's type relationship with the connected porosity with a porosity (cementation) exponent of m = 1.70 ± 0.02, much smaller than for volcanic extrusive rocks. Like for volcanic rocks, a dynamic Stern layer model can be used to illustrate the behavior associated with the clay-minerals (mostly kaolinite). A field investigation is conducted in the Vosges (France) using a deep time-domain induced polarization survey reaching at a depth of investigation ∼400 m. We show how the electrical conductivity and the normalized chargeability can be used to image the water content and CEC of the granitic substratum. The conductivity of granite is found to be dominated by surface conductivity rather than by bulk conductivity and therefore Archie's law cannot be used as a conductivity equation to interpret field data as commonly done in ElectroMagnetic (EM) surveys.

SummaryWe present a method for ambient noise cross-correlation modelling and source inversion, which accounts for spatio-spectral variability in noise source distributions. It is based on numerical wavefield simulations in 2-D acoustic media. The source power spectral density is parameterized by a sum of a small number of spatial source distributions, each with a corresponding frequency spectrum held fixed during the inversion. Algorithmically, this is an extension of our previous work which assumed spatially homogeneous source spectra. In this paper, we use it to study the impact of incorrectly estimating source spectra from observed data. This is done using synthetic tests involving sources with closely spaced frequency spectra. The tests demonstrate that when the spatial variability of sources is either partially or wholly unaccounted for, the recovery of true source locations is compromised.

SummaryWe determine detailed 3-D velocity and anisotropy models of the crust beneath the central segment of the Tanlu fault zone (TLFZ) in east China using hand-picked arrival times of direct waves (Pg, Sg), refracted waves (Pn, Sn), and the Moho reflected waves (PmP and SmS) of local earthquakes recorded at 55 portable seismic stations of our newly deployed TanluArray and 45 Chinese provincial stations. Our results show that the pattern of seismic velocity and anisotropy is in good agreement with surface geological and tectonic features in the study region. Along the TLFZ, obvious low-velocity (low-V) anomalies are visible in the lower crust in the north of the Cangni fault, whereas high-velocity (high-V) anomalies appear in the upper and middle crust. The 1668 Tancheng earthquake (M8.5) and small-to-moderate earthquakes occurred in the high-V zones but underlain by low-V anomalies. The low-V zones may reflect significant effects of fluids derived from hot and wet upwelling flow in the big mantle wedge. Under the Suqian seismic gap, low-V anomalies are clearly imaged in the middle crust, whereas high-V anomalies appear in the lower crust, suggesting the absence of fluids in the deep crust, which could explain why few earthquakes occurred there. The fast velocity direction (FVD) of P-wave azimuthal anisotropy changes with depth. The FVD is parallel with the TLFZ strike in the upper crust in and around the TLFZ, reflecting the strike-slip features of the fault zone. Our results provide new seismic constraints on the structure and tectonic evolution of the TLFZ.

SummaryThe Tjörnes Fracture Zone in North Iceland is one of two transform zones in Iceland capable of generating earthquakes of magnitude ∼7. More than 150 years have passed since the last two major earthquakes occurred on the Húsavík-Flatey fault, one of the two main transform structures within the Tjörnes Fracture Zone. Given the seismic hazard posed to Húsavík and adjacent coastal communities, accurately determining the slip rate and locking depth of this fault is crucial for a robust earthquake hazard and risk assessment. In this study, we significantly expand the existing GNSS dataset for the Tjörnes Fracture Zone by incorporating more than a decade of additional data and doubling the number of stations. This expansion not only improves the spatial coverage of the network, but also refines the station velocities. We present an updated interseismic velocity field for North Iceland and implement a backslip model with nine dislocation segments to describe the plate boundary deformation. Additionally, we include a point pressure source for the ongoing broad uplift signal in the study area. Our findings indicate a locking depth of $7.3\substack{+0.9 -0.7}$ km and an average slip rate of 6.9 ± 0.2 mm/yr for the Húsavík-Flatey fault. With our updated approach, we can narrow down model parameter constraints from previous studies and thereby provide an enhanced understanding of the earthquake potential of this region.

SummaryUsing the published paleoseismological trenching data for 16 faults in Central Italy, we compile a new database of surface faulting earthquakes, having a quite stationary temporal distribution since 6000 BCE. By applying a probabilistic aggregation method, we correlate the event ages from distinct trenches on each fault, to construct all possible individual fault rupture scenarios, consistent with geological constraints. These inferred fault time histories are the basis for both individual fault and regional seismic hazard evaluation. We found that the mean recurrence time of each fault goes from about 1 to 4 thousand years for individual faults, whereas the value at regional scale is close to 120 yrs. The small size of individual fault data samples does not allow us to infer straightforward information on the fault temporal behavior, but only to evaluate the reliability of a chosen occurrence model for each fault. Therefore, hazard assessment is carried out by including the uncertainties related to both ages and probability distribution of the inter-event times. We find that both these have a large impact on the probabilities of next rupture for individual faults: these depend on basic features of the temporal model and on the relation between the elapsed time and the mean interevent time. At a regional scale, we cannot exclude the simplest possible model, i.e. the poissonian behavior, that provides quite stable probabilities of future events, close to 27% in the next 50 years.

SummaryHigh-resolution seismic tomographic images of active fault zones are essential for linking physical processes and observations during large ruptures. The 1999 Izmit and Düzce earthquakes offer a valuable opportunity to study local fault geology, fault structure, and rupture characteristics. By analyzing seismic data from the aftershocks of these earthquakes and long-term seismic observations in the region, I computed P-wave velocity variations along ∼160 km long ruptured segments of the NAF at km-scale resolution.The large velocity perturbations (±10 per cent) at depths of less than 5 km are associated with known geological structures alongside the fault zone. Lower velocities are observed in the Çınarcık, Izmit, Adapazarı, and Düzce basins, while higher velocities are associated with the metamorphic and ultramafic rocks of the Armutlu-Almacık zone.A sharp velocity contrast (±10 per cent) is observed along the 1999 Izmit rupture zone at depths of less than 5 km, but this contrast diminishes at greater depths. The Izmit rupture zone does not hold a large-scale bi-material interface but instead exhibits heterogeneties on shorter length scales. The low-velocity fault zones, ∼3 km thick, detected only in a limited section of the rupture zone.The Almacık block, with its high-velocity core (6.0–6.6 km/s) extending to a depth of about 13 km, plays a crucial role in shaping the multiple branches of the NAF and in the partitioning of strain. The rupture arrests of both 1999 Izmit and Düzce earthquakes occurred in the transition of Düzce and Karadere faults, where a transtensional structure develops at depth of 8-12 km within the zone along the northern boundary of the Almacık block.There appears to be a link between the coseismic slip during the Izmit earthquake and P-wave velocity gradient along the fault zone. Nonetheless, this apparent correlation requires further evaluation through dynamic rupture simulations.

SummaryA problem of growing importance in earthquake forecasting is how to compare probabilistic forecasting models with deterministic physical simulations and extract physical insights from their differences. Here we compare the time-independent Uniform California Earthquake Rupture Forecast Version 3 with a long earthquake catalog simulated by the multi-cycle Rate-State Quake Simulator (RSQSim). Shaw et al. (2018) generated a million-year rupture catalog for California from RSQSim simulations based on UCERF3 fault geometries and slip rates and found that the shaking hazard from the synthetic catalog was in good agreement with the UCERF3 hazard maps. We take this model-to-model comparison to the more granular level of individual faults and ruptures. We map RSQSim ruptures from the Shaw18 catalog onto equivalent UCERF3 ruptures by maximizing the mapping efficiency and ensuring that every RSQSim realizations is associated with a unique UCERF3 rupture. The full UCERF3 logic tree is used to approximate the prior distributions of individual rupture rates and fault subsection participation rates as independent gamma distributions. We formally test the ontological null hypothesis (ONH) that the empirical RSQSim rupture counts are statistically consistent with the UCERF3 rate distributions, given the sampling uncertainty of the RSQSim catalog and the epistemic uncertainty of the UCERF3 model. Testing individual rupture rates provides little evidence either for or against the ONH, owing to the predominance of large ruptures with low recurrence rates. However, at the subsection level, the statistically significant discrepancies are much more common than expected under the ONH. We obtain a 25% failure rate at the 5% significance level and a 15% failure rate at 1% level. The false discovery rates estimated by q-value calculations are low, so we can be confident that the same subsections would likely fail if tested against an independent million-year catalog generated by the same RSQSim model. Bayesian recalibration of the UCERF3 priors using the empirical RSQSim rates yields Gamma posterior distributions that can be derived analytically. The results of testing and recalibration, taken together, quantify how well RSQSim rupture rates agree with, and differ from, the UCERF3 forecast rates. We find that some of the discrepancies can be attributed to the differences in slip rates that drive the models, whereas others are governed by the RSQSim fault dynamics absent from UCERF3.