Geophysical Journal International

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.

SummaryAdvance detection is a pivotal technology in tunnel construction, designed to precisely invert the velocity distribution of geological formations, thereby ensuring both construction safety and operational efficiency. The unscented hybrid simulated annealing (UHSA) algorithm is a global optimization technique that has been successfully applied to the advance detection of tunnels. However, the UHSA exhibits a slow convergence rate under complex geological conditions and is prone to getting trapped in local optima. To address this issue, we propose an improved method based on the whale optimization algorithm (WOA), referred to as MEWOA. MEWOA incorporates a population-initialization method based on the tent map and reverse learning and integrates a nonlinear convergence factor, a frequency-fluctuation mechanism, and the low-soaring strategy of the red-tailed hawk algorithm (RTH). These enhancements notably improve the accuracy and convergence speed of the algorithm. We conducted a qualitative analysis of the enhancement mechanisms in MEWOA using two functions and performed quantitative experiments on four tunnel models. The experimental results demonstrate that MEWOA outperforms UHSA, achieving higher accuracy in solving inversion problems with multiple velocity models.

SummaryAs probabilistic tsunami hazard analysis (PTHA) focuses more on assessments for localized, populous regions, techniques are needed to identify a subsample of representative earthquake ruptures to make the computational requirements for producing high-resolution hazard maps tractable. Moreover, the greatest epistemic uncertainty in seismic PTHA is related to source characterization, which is often poorly defined and subjective. We address these two salient issues by applying streamlined earthquake rupture forecasts (ERFs), based on combinatorial optimization methods, to an unsupervised machine learning workflow for identifying representative ruptures. ERFs determine the optimal distribution of a millennia-scale sample of earthquakes by inverting the observed slip rate on major faults. We use two previously developed combinatorial optimization ERFs, integer programming and greedy sequential, to produce the optimal location of ruptures with seismic moments sampled from a regional Gutenberg-Richter magnitude-frequency distribution. These ruptures in turn are used to calculate peak nearshore tsunami amplitude, using computationally efficient tsunami Green's functions. An unsupervised machine learning workflow is then used to identify a small sub-sample of the earthquakes input to ERFs for onshore PTHA analysis. We eliminate epistemic uncertainty related to source distribution under traditional PTHA analysis; in its place, a quantifiable, less subjective, and generally smaller uncertainty related to the input to ERFs is included. The Nankai subduction zone is used as a test case, where previous ERFs have been conducted. Results indicate that the locations of representative earthquakes are sensitive to choice of magnitude-area relation and to whether a minimum cumulative stress objective is imposed on the fault. In general, incorporating ERFs into PTHA provide a physically self-consistent method to incorporate fault slip information in determining representative earthquakes for onshore PTHA, eliminating a major source of epistemic uncertainty.

SummaryCurrent efforts to correctly categorize natural events from suspected explosion sources with data that is collected by ground- or space-based sensors presents historical challenges that remain unaddressed by the Event Categorization Matrix (ECM) model. Smaller historical events (lower-yield explosions) may have data available from fewer measurement techniques than are available today, and therefore, a historical event record can lack a complete set of discriminants. The covariance structures can also differ between such observations of event (source-type) categories. Both obstacles are problematic for the classic Event Categorization Matrix model. Our work addresses this gap and presents a Bayesian update to the previous Event Categorization Matrix model, termed the Bayesian Event Categorization Matrix model, which can be trained on partial observations and does not rely on a pooled covariance structure. We further augment the Event Categorization Matrix model with Bayesian Decision Theory so that false negative or false positive rates of an event categorization can be reduced in an intuitive manner. To demonstrate improved categorization rates for the Bayesian Event Categorization Matrix model, we compare an array of Bayesian and classic models with multiple performance metrics using Monte Carlo experiments. We use both synthetic and real data. Our Bayesian models show consistent gains in overall accuracy and lower false negative rates relative to the classic Event Categorization Matrix model. We propose future avenues to improve Bayesian Event Categorization Matrix models’ decision making and predictive capability.

SummaryThe response of porous rocks to fluid flow and external loads is critical to a wide range of geophysical and geotechnical problems, and is described by the long-established theory of poroelasticity. Despite the wealth of literature existing on the subject, little investigation has been devoted to the analysis of the state of stress and pore pressure regimes within porous media under the action of gravity force. Here, we present analytical solutions for the effective stress and pore pressure in gravitationally loaded porous rocks, in both drained and undrained conditions, and compare them with results from Finite Element numerical models. We also apply our models to a ground deformation episode observed in the Campi Flegrei caldera, Italy. We find that the numerical results accurately reproduce our analytical solutions and show how accounting for gravity-induced stress and pore pressure regimes is critical to model stress and deformation in poroelastic media accurately. Specifically, we highlight how failing to assign realistic initial conditions, or neglecting gravity altogether, may lead to misleading results and interpretations of geophysical observables, such as ground deformation and rock failure.

SummaryPrecursory spurious arrivals, commonly observed in ambient noise correlations, are generated by near-field noise sources. We have developed an inversion method to evaluate the noise source distribution based on the precursory waves. This method is applied to the BASIN experiment in Los Angeles, revealing that the noise sources show coherent patterns with features like faults and structure boundaries. Our spectral analysis indicates that the energy of the noise source in generating the precursory signal is predominant at higher frequencies, suggesting a shallow source(<200m). We conclude that near-field noise is primarily produced by scattering from geological structures with significant velocity contrasts, such as faults at shallow depths. This method offers a new way to map faults using ambient noise correlations.

SummaryThere are currently few studies in the literature on source and background seismic noise power distribution in Italy. In this research, the seismic noise recorded by 233 broadband (BB) seismic stations of the Italian Seismic Network (ISN), operating continuously between 2015 - 2018, was investigated. Starting from the average Power Spectral Density (PSD) calculated for each selected station, the seismic noise powers in four subsets from 0.025 to 30 Hz frequency bands were analysed. Using the Inverse Distance Weighted (IDW) interpolation method, the background noise distribution in the entire Italian territory was observed, producing noise interpolation maps for each of the four frequency bands. Furthermore, regional seismic noise models and local anomalies for each subset of the frequency bands were found by applying 2D moving average spatial filtering techniques. In addition, in this research, a first attempt of linear regression analysis is performed to discover possible relationships between seismic noise powers and geographical parameters (elevation site and minimal station-coast distance) and the role of meteorological parameters such as rainfall. The large dataset obtained allowed for the assessment of the main characteristics and sources of seismic noise at all sites. One can thus know the distribution of noise levels in the Italian territory and in particular study their main sources related to natural and anthropic ambient vibrations. To improve the analysis, a comparison between the seismic noise maps and the completeness magnitude map was carried out, which showed the effectiveness of the national seismic network.

SummaryThe Seoul Mega City (SMC), one of the most densely populated regions in the world, has experienced an increase in seismic activity, raising concerns about the potential presence of concealed faults that could trigger clustered earthquakes in the Seongbuk area (SB), located in the central part of the SMC. This study aims to identify such a fault through the interpretation of terrestrially measured gravity field. By employing dip-curvature analysis and the first vertical derivative calculation of the gravity field, supported by S-wave velocity models, reflection seismic profiles, geological data, and borehole data analyses, we provide compelling evidence for the existence of a deeply buried NS-trending fault (DF0) associated with the Dongducheon Fault system. This fault is likely responsible for the clustered seismic activities observed in the SB. The confirmation of DF0’s presence highlights the critical need for further geophysical investigations to better understand seismic risks in the SMC, which has a history of significant seismic events.