Geophysical Journal International

SummaryWe revisit the budget of 20th century true polar wander (∼1°/Myr in the direction of 70°W) using a state-of-the-art adjoint-based reconstruction of mantle convective flow and predictions of ongoing glacial isostatic adjustment that adopt two independent models of Pleistocene ice history. Both calculations are based on a mantle viscosity profile that simultaneously fits a suite of data sets related to glacial isostatic adjustment (Fennoscandian Relaxation Spectrum, post-glacial decay times) and a set of present-day observations associated with mantle convection (long-wavelength gravity-anomalies, plate motions, excess ellipticity of the core-mantle boundary). Our predictions reconcile both the magnitude and direction of the observed true polar wander rate, with convection and glacial isostatic adjustment contributing signals that are 25-30% and ∼75% of the observed rate, respectively. The former assumes that large-scale seismic velocity heterogeneities are purely thermal in origin, and we argue that our estimate of the convection signal likely represents an upper bound due to the neglect of hypothesized compositional variations within the large low shear velocity provinces in the deep mantle.

AbstractThe earthquake size distribution is well described by the Gutenberg Richter Law, controlled by the b-value parameter. In recent decades, a great variety of methods for estimating the b-value have been proposed by the scientific community, despite the simplicity of this relationship. All these methods underlie the different views of individual modelers and, therefore, often generate inconsistent results. In this study, we perform a seismological experiment in which we compare different, commonly adopted, methodologies, to estimate the completeness magnitude and the b-value, for seismicity in Central Italy. The inter-method differences are on average equal to 0.4 and 0.3, for Mc and b, respectively, but reach much larger values, especially during more intense seismic activity. This shows that epistemic uncertainty in the b-value plays a more crucial role than intra-method uncertainties, opening new perspectives in the interpretation of discrepant, single studies.

SummaryGeological Carbon Storage (GCS) is one of the most viable climate-change mitigating net-negative CO2-emission technologies for large-scale CO2 sequestration. However, subsurface complexities and reservoir heterogeneity demand a systematic approach to uncertainty quantification to ensure both containment and conformance, as well as to optimize operations. As a step toward a Digital Twin for monitoring and control of underground storage, we introduce a new machine-learning-based data-assimilation framework validated on realistic numerical simulations. The proposed Digital Shadow combines Simulation-Based Inference (SBI) with a novel neural adaptation of a recently developed nonlinear ensemble filtering technique. To characterize the posterior distribution of CO2 plume states (saturation and pressure) conditioned on multimodal time-lapse data, consisting of imaged surface seismic and well-log data, a generic recursive scheme is employed, where neural networks are trained on simulated ensembles for the time-advanced state and observations. Once trained, the Digital Shadow infers the state as time-lapse field data become available. Unlike ensemble Kalman filtering, corrections to predicted states are computed via a learned nonlinear prior-to-posterior mapping that supports non-Gaussian statistics and nonlinear models for the dynamics and observations. Training and inference are facilitated by the combined use of conditional invertible neural networks and bespoke physics-based summary statistics. Starting with a probabilistic permeability model derived from a baseline seismic survey, the Digital Shadow is validated against unseen simulated ground-truth time-lapse data. Results show that injection-site-specific uncertainty in permeability can be incorporated into state uncertainty, and the highest reconstruction quality is achieved when conditioning on both seismic and wellbore data. Despite incomplete permeability knowledge, the Digital Shadow accurately tracks the subsurface state throughout a realistic CO2 injection project. This work establishes the first proof-of-concept for an uncertainty-aware, scalable Digital Shadow, laying the foundation for a Digital Twin to optimize underground storage operations.

SummaryThe Ecuadorian forearc, formed by the accretion of oceanic plateaus, island arcs, subduction of an aseismic ridge, records a history of long-lived subduction. The modern system includes subduction of the Carnegie Ridge and seamounts, young forearc coastal ranges, and translation of a forearc sliver from oblique subduction of the Nazca Plate beneath South America. The margin has experienced large megathrust earthquakes and exhibits slow-slip events and earthquake swarms. We present results from joint tomographic inversion of local earthquakes for 3D velocity structure and earthquake location. Our joint inversion uses seismic arrival-time data from local earthquakes recorded by permanent stations and dense seismic temporary networks deployed near the coast after the 2016 Mw 7.8 Pedernales megathrust rupture and across the entire northern forearc into the foothills of the Andes in 2021-2022. Our results show that seismicity distribution and megathrust rupture are controlled by inherited and modern structures in the upper plate forearc and subducting Nazca Plate. Forearc sedimentary basins observed as low-velocities (Vp < 5.8 km/s, Vs < 3.2 km/s) are dissected by forearc basement highs observed as fast velocities (Vp 6.6-7.2 km/s, Vs. 3.6-4.0 km/s). Localized deep depocenters adjacent to basement highs preserve older sedimentary sections beneath younger forearc deposits. Differences in velocity allow discrimination between oceanic plateau basement associated with the Piñón terrane beneath the forearc and accreted island arc terranes along the eastern forearc boundary with the Andes. Along the coast, basement velocities are consistent with a hydrated upper plate. We observe an apparent transient in Vp/Vs (higher to lower) in the upper plate after the 2016 megathrust rupture, representing a transient flux of fluids from the subducting slab into the upper plate triggered by the earthquake. We observe variable thickness of the subducting Nazca plate from ∼10 km north of the Carnegie Ridge reaching 20-25 km where the Carnegie Ridge subducts beneath the forearc. Lateral velocity variations in the subducting plate indicate heterogeneity along strike and dip associated with magmatic evolution of the ridge. High-velocity domains at depth correlate with seamounts and subducted relief along the Carnegie Ridge. A low-velocity zone marks the boundary between the subducting and overriding plates. The downdip termination of the Pedernales megathrust rupture coincides with structure of the Carnegie Ridge and along strike changes in the plate interface. The downdip edge of the rupture occurs where the low-velocity zone is absent, and the subducting Carnegie Ridge intersects the overlying mantle wedge. Earthquakes located with the joint inversion focus into tight clusters controlled by relief at the top of the subducting slab and basement structure in the overriding plate. Along the coast, seismicity shallows from south to north across the east-west striking Canandé Fault. South of the fault, seismicity locates predominantly within the subducting plate and plate interface. To the north, seismicity concentrates within the plate interface and upper plate. The northward shift in hypocenter depths and an offset in the eastern limit of thick subducting Nazca plate across the Canandé fault marks a significant transition in the forearc across the fault.

SummaryCharacterizing the hydraulic and geomechanical behavior of crystalline rocks is of importance for a wide range of geological and engineering applications. Geophysical methods in general and seismic techniques in particular are extensively used for these purposes due to their cost-effective and non-invasive nature. In this study, we combine legacy seismic observations to analyze the seismic attenuation and velocity characteristics in macroscopically intact regions of the granodiorite hosting the underground Grimsel Test Site in the central Swiss Alps across a wide frequency range. By focusing on data from the intact rock volumes we aim to assess the importance of viscoelastic effects in the crystalline host rock. Our results show consistent frequency-dependent characteristics of the seismic velocity and attenuation. We illustrate that it is possible to fit a microcrack-related wave-induced fluid flow (WIFF) model to the data over the entire frequency spectrum under examination extending from the Hertz to the Megahertz range. Utilizing complementary pressure-dependent ultrasonic measurements, we infer microcrack properties that validate the key parameters of the proposed WIFF model. These findings deepen our understanding of dispersion and attenuation mechanisms at the microscopic scale in crystalline environments, which is critical for a coherent analysis and integration of data from different seismic techniques as well as for the identification of dispersion and attenuation effects related to macroscale heterogeneities, such as fractures and faults.

SummaryThe South China block (SCB) formed after the collision between the Yangtze craton and the Cathaysia block along the Jiangnan orogenic belt at 980–820 Ma. Afterward, intense tectonic activity occurred in the SCB in the early Paleozoic and early Mesozoic. In the Mid–Late Jurassic and Early Cretaceous, the SCB experienced vigorous magmatic activity, which resulted in assemblage of mineralogenetic materials in the Youjiang basin (YB), Southeast coastal metallogenic belt (SCMB) and Wuyishan metallogenic belt (WMB). The mechanisms involved in the formation of various types of metallic ore deposits have attracted considerable attention. However, the crustal and uppermost mantle dynamics of the metallogenic mechanisms are still controversial. To address this issue, we conducted seismic tomography to image the velocity and azimuthal anisotropy of the crust and uppermost mantle beneath the SCB. In this study, an eikonal equation-based traveltime tomography method was used to invert a total of 143,473 high-quality P-wave first arrivals, which were obtained by manually picking the seismic waveforms of 3615 regional earthquakes recorded by 892 broadband seismic stations. After the inversion, we derived high-resolution images, in which we identified a strong low-velocity anomaly and weak azimuthal anisotropy in the uppermost mantle of the northern YB. Below the SCMB, a low-velocity body extends from the uppermost mantle to the bottom of the crust; the azimuthal anisotropy of the uppermost mantle is weak and does not exhibit a consistent fast-velocity direction (FVD). These characteristics can be attributed to the upwelling of hot materials and crustal partial melting. For both the northern YB and SCMB, the low-velocity anomaly is probably related to hot property of the uppermost mantle and weak azimuthal anisotropy may be due to the nearly vertical α-axis of olivine. These features indicate the upwelling of hot materials beneath the YB. The upwelling of hot materials carried metal elements from deep mantle to shallow crust, resulting in metal deposits. Beneath the WMB, the lower crust and uppermost mantle show high-velocity anomalies and moderate strong azimuthal anisotropy with a consistent NE–SW-oriented FVD. The high-velocity anomalies reflect cold and rigid properties of the lower crust and the uppermost mantle beneath the WMB; consistent FVD of azimuthal anisotropy may indicate ancient fossil anisotropy. These features suggest ancient continental relicts of the Cathaysia block under the WMB.

SummaryNumerical simulations of earthquakes and seismic wave propagation require accurate material models of the solid Earth. In contrast to purely elastic rheology, poroelasticity accounts for pore fluid pressure and fluid flow in porous media. Poroelastic effects can alter both the seismic wave field and the dynamic rupture characteristics of earthquakes. For example, the presence of fluids may affect cascading multi-fault ruptures, potentially leading to larger-than-expected earthquakes. However, incorporating poroelastic coupling into the elastodynamic wave equations increases the computational complexity of numerical simulations compared to elastic or viscoelastic material models, as the underlying partial differential equations become stiff. In this study, we use a Discontinuous Galerkin solver with Arbitrary High-Order DERivative time stepping (ADER-DG) of the poroelastic wave equations implemented in the open-source software SeisSol to simulate 3D complex seismic wave propagation and 3D dynamic rupture in poroelastic media. We verify our approach for double-couple point sources using independent methods including a semi-analytical solution and a finite-difference scheme and a homogeneous full-space and a poroelastic layer-over-half-space model, respectively. In a realistic carbon capture and storage (CCS) reservoir scenario at the Sleipner site in the Utsira Formation, Norway, we model 3D wave propagation through poroelastic sandstone layers separated by impermeable shale. Our results show a sudden change in the pressure field across material interfaces, which manifests as a discontinuity when viewed at the length scale of the dominant wavelengths of S- or fast P-waves. Accurately resolving the resulting steep pressure gradient dramatically increases the computational demands, requiring high-resolution modeling. We show that the Gassmann elastic equivalent model yields almost identical results to the fully poroelastic model when focusing solely on solid particle velocities. We extend this approach using suitable numerical fluxes to 3D dynamic rupture simulations in complex fault systems, presenting the first 3D scenarios that combine poroelastic media with geometrically complex, multi-fault rupture dynamics and tetrahedral meshes. Our findings reveal that, in contrast to modeling wave propagation only, poroelastic materials significantly alter rupture characteristics compared to using elastic equivalent media since the elastic equivalent fails to capture the evolution of pore pressure. Particularly in fault branching scenarios, the Biot coefficient plays a key role in either promoting or inhibiting fault activation. In some cases, ruptures are diverted to secondary faults, while in others, poroelastic effects induce rupture arrest. In a fault zone dynamic rupture model, we find poroelasticity aiding pulse-like rupture. A healing front is induced by the reduced pore pressure due to reflected waves from the boundaries of the poroelastic damage zone. Our results highlight that poroelastic effects are important for realistic simulations of seismic waves and earthquake rupture dynamics. In particular, our poroelastic simulations may offer new insights on the complexity of multi-fault rupture dynamics, fault-to-fault interaction and seismic wave propagation in realistic models of the Earth’s subsurface.

SummaryIn this work we apply a dedicated 4D full waveform inversion workflow to short offset streamer data from the Sleipner CO2 storage field in the North Sea. We consider a baseline dataset acquired in 1994 and a monitor dataset acquired in 2008. Accessing to only short offset data raises significant difficulties for full waveform inversion. In this case the penetration of diving waves, which controls the depth where quantitative updates of the velocity can be expected, do not reach the zone of interest where the CO2 is injected. For this reason, we propose to combine an efficient time-lapse full waveform inversion strategy, which we call simultaneous, with a reflection oriented full waveform inversion workflow. The latter has been introduced in the literature as a way to circumvent short-offset limitation and increase the ability of full waveform inversion to update the macro-velocity model at depth by exploiting the reflection paths, using a prior step of impedance reconstruction. We first illustrate the interest of this combined strategy on a 2D synthetic model inspired from the Sleipner area. Then we apply it to the Sleipner field data, using as baseline model the one we present in a companion paper, where our reflection oriented workflow is presented. Our combined approach yields reliable estimates of the changes due to the CO2 injection, characterized by velocity reductions of up to 400 m.s−1 and strong impedance contrasts at depths of 800-1000 m, which consistent with previous FWI studies. Furthermore, the spatial distribution of CO2 changes aligns with conventional seismic time-migration results from earlier studies, following a north-south migration trend.

SummaryTeleseismic P-wave coda autocorrelation has been increasingly applied to subsurface structure characterization, given its potential to infer velocities. However, the inversion of coda autocorrelation data has not been extensively investigated regarding data processing (stacking and move-out correction), inversion approaches (Monte Carlo or metaheuristic), model parameterization, and applicability. Here, we propose an inversion method for teleseismic P-wave coda autocorrelation based on particle swarm optimization and a treatment of uncertainty. This inversion method utilizes the arrival time information of reflected (or converted) waves contained in the binned stack waveforms, demonstrating promising model adaptability and robustness. Synthetic data tests show that this method accurately inverted various geological models without prior information, such as the number of crustal layers, surface sedimentary layers, and low-velocity zones within the crust. The method was successfully applied to the QSPA station near the South Pole, revealing an ice sheet thickness of approximately 2900 m, with a 340 m thick low shear-wave velocity ice layer at the base, likely containing up to 15% water. Beneath the ice sheet, we infer a 400 m thick subglacial sediment layer. The uncertainties of the thickness of the low shear wave velocity ice and the sedimentary layer are 150 m and 10 m, respectively. These findings and the potential of the proposed method open up new directions for glacier dynamics research in the region. Additionally, we apply the method to the BOSA station near Kimberley, South Africa, which confirms clear Moho and intracrustal interfaces, consistent with receiver functions and deep seismic reflection data results. This study improves the inversion algorithm for teleseismic P-wave coda autocorrelation and expands its application scenarios.

SummaryGeophysical simulations for complex subsurface structures and material distributions require the evaluation of partial differential equations by means of numerical methods. However, the mentioned high complexity often yields computationally very costly simulations, especially for electromagnetic (EM) and seismic methods. When used in the context of parameter estimation or inversion studies, this aspect severely limits the number of simulations that are affordable. However, especially for structured model analysis methods, such as global sensitivity analyses or inversions, often thousands to millions of forward simulation runs are required. To address this challenge, we propose utilizing a physics-based machine learning method, namely the non-intrusive reduced basis method, aiming at constructing low-dimensional surrogate models to significantly reduce the computational cost associated with the numerical forward model while preserving the physical principles. We demonstrate the effectiveness and benefits of the surrogate models using broadband Magnetotelluric (MT) responses of a 2-D model that mimics a conceptual volcano-hosted geothermal system. Next to being a first such application, we also show how ML reduced basis method can be adapted to consistently treat complex-valued variables – an aspect that has been overlooked in previous studies. Additionally, reducing computation time by several orders of magnitude through the surrogate enables us to perform a global sensitivity analysis for MT applications. Despite additional insights, such an analysis has been normally deemed infeasible given the high computational burden. The methods developed here are presented in a generalized form, making this approach feasible for other electromagnetic techniques with a low-dimensional parameter space.

SummaryThe triple junction between the North American, Caribbean, and Cocos plates at the Guatemala-Mexico border is not well understood. It forms a broad region from around the active Tacaná volcano up to the Guatemala City graben. Tacaná is the westernmost active volcano of the Central American volcanic arc and is located at the intersection of four major active faults: the Polochic, Motagua, Jalpatagua, and Tonalá faults. Using seismicity around the Tacaná volcano, we show that there is moderate to low tectonic seismic activity between the Guatemala City graben and the Tacaná volcano, possibly related to the ancient extremes of the Motagua and Jalpatagua faults. Therefore, we speculate that the triple junction would be located onshore, around the Tacaná volcano.We located earthquakes around the Tacaná volcano between January 2017 and October 2018, a period that includes the large Mw8.2 Tehuantepec (Chiapas) earthquake of 8 September 2017, located ∼190 km away. We identified four distinct types of seismicity, interpreted as having tectonic, hydrothermal, intermediate depth magmatic, and deep magmatic origins. The tectonic seismicity occurred at depths between ∼5 km and ∼30 km b.s.l., and may be associated with three faults around the Tacaná volcanic complex. These faults are oriented in NE-SW, aligned with the four Tacaná volcanic edifices; NW-SE, consistent with the Jalpatagua fault; and approximately EW, corresponding to the Motagua fault. The hydrothermal seismicity is observed at shallow depths, from the subsurface to about 2 km b.s.l., predominantly in the western sector of the Tacaná summit, and partially beneath the San Antonio volcano, an area known for intense hydrothermal activity. This seismicity is spatially related to the shallow portions of the same three faults described above. The intermediate depth magmatic seismicity is detected at depths between 5 and 12 km b.s.l. and is interpreted to be related to the presence of a shallow magma chamber beneath the Tacaná volcanic complex. Finally, the deep magmatic seismicity is located in the eastern part of the Tacaná, at depths ranging from 15 km to about 22 km b.s.l. This seismicity is interpreted to be due to a vertical dike intrusion that connects a deep magma reservoir located between 30 km and 40 km depth, to the hypothesized shallower magma chamber associated with the intermediate depth seismicity.