Geophysical Journal International

SummarySeismic swarms are known to occur in regions with complex deformation and multiple fault systems. The identification of the affected faults, the evolution of the seismicity and the rupture characteristics are key to understand seismotectonics and seismic hazard of such areas. We here address the Atarfe-Santa Fe earthquake swarm with > 5000 events, including six magnitude 4+ earthquakes, recorded in 2021 in the Granada Basin area (S-Spain). We use continuous data from a dense local recording network and apply deep learning models to pick and associate seismic phase arrivals and construct an automatic event catalogue. A double difference relocation of 3196 earthquakes reveals the seismotectonic fine-structure of the swarm. We identify planar, southwest-dipping structures among the relocated hypocentres, consistent with the NW-SE trending, high-angle normal faults present in this sector of the Granada basin. Earthquakes concentrate between 4-7 km depth along three different lineaments. The distinctive pattern of three equidistant, parallel segments allows an association of the swarm with the Ermita los Tres Juanes, Atarfe and Pinos Puente normal faults. These faults outcrop at the upward extrapolation of the swarm, forming an arrangement of three structures that mimic the geometry of the relocated seismicity at depth. In the course of the swarm, the seismicity jumped from the Ermita los Tres Juanes to the Atarfe fault, in midst of a rapid succession of three magnitude 4 earthquakes within 20 minutes, then migrated laterally along both faults, and later migrated basinwards to the Pinos Puente fault, which produced fewer and smaller events. We estimate apparent source time functions for five earthquakes (Mw 4.1 to 4.4) through the deconvolution of empirical Green's functions from the records, suggesting rupture propagation towards NW, N and E directions. An isochrone back projection of apparent source time functions suggests km-scale ruptures with simple slip distributions, showing lateral and updip components of rupture. Our results shed light on the complexity of this seismic swarm in terms of the fault network involved, the propagation of seismicity across the faults and the variable directions of individual ruptures.

SummarySeismic tomography is routinely used to image the Earth’s interior using seismic data. However, in practice, data limitations lead to discretised inversions or the use of regularisations, which complicates tomographic model interpretations. In contrast, Backus-Gilbert inference methods make it possible to infer properties of the true Earth, providing useful insights into the internal structure of our planet. Two related branches of inference methods have been developed – the Subtractive Optimally Localized Averages (SOLA) method and Deterministic Linear Inference (DLI) approaches – each with their own advantages and limitations. In this contribution, we show how the two branches can be combined to derive a new framework for inference, which we refer to as SOLA-DLI. SOLA-DLI retains the advantages of both branches: it enables us to interpret results through the target kernels, rather than the imperfect resolving kernels, while also using the resolving kernels to inform us on trade-offs between physical parameters. We therefore highlight the importance and benefits of a more careful consideration of the target kernels. This also allows us to build families of models, rather than just constraining properties, using these inference methods. We illustrate the advantages of SOLA-DLI using three case studies, assuming error-free data at present. In the first, we illustrate how properties such as different local averages and gradients can be obtained, including associated bounds on these properties and resolution information. Our second case study shows how resolution analysis and trade-offs between physical parameters can be analysed using SOLA-DLI, even when no data values or errors are available. Using our final case study, we demonstrate that SOLA-DLI can be utilised to obtain bounds on the coefficients of basis function expansions, which leads to discretised models with specific advantages compared to classical least-squares solutions. Future work will focus on including data errors in the same framework. This publication is accompanied by a SOLA-DLI software package that allows the interested reader to reproduce our results and to utilise the method for their own research.

SummaryDetermining the precise pattern of crustal deformation enhances our comprehension of crustal deformation traits and the significant earthquakes. By incorporating 21 additional continuous GNSS stations along with existing ones, we generated an updated GNSS velocity field for the northeastern Tibetan Plateau. Using the back-slip dislocation model, we calculated the average slip rates of three major active faults: the Haiyuan fault, the Liupanshan fault, and the Helanshan fault. Our findings indicated that the regional crustal movement does not conform to the equilibrium principles typically associated with the triple junction-like tectonics. This suggests the existence of a newly active tectonic belt within the Longxi block. Consequently, we proposed a revised block model that incorporates a right-lateral shear zone within the Longxi block to account for the observed crustal deformation in the northeastern Tibetan Plateau. Our study indicates that the right-lateral shear zone significantly contributes to the northeastward expansion of the Tibetan Plateau, accounting for approximately 82 per cent of strain accumulation, while the remaining 18 per cent accumulates along the Liupanshan fault. The revised block model emphasizes the pivotal role of the Haiyuan fault and the right-lateral shear belt as the key tectonic factors shaping the crustal deformation pattern. Our result enables a comprehensive understanding of both the spatial variations observed in the GNSS velocity field and the spatial distribution of significant earthquakes in the region.

SummaryFull waveform inversion (FWI) has enjoyed increased attention the past decade, becoming the state of the art for estimating parameters influencing wave propagation in a medium. However, only a few recent emerging efforts have attempted to tackle the challenge of uncertainty quantification in FWI. In this study, we suggest joining FWI with the Bayesian approach, where we provide a post-processing step with an advantageous starting point defined by the global minimum stemming from a deterministic FWI algorithm. Then, using the local ensemble transform Kalman filter (LETKF), we obtain the uncertainty as a follow-up step to the FWI procedure. Within a probabilistic Bayesian inversion framework, the LETKF uses local seismic data to update sets of variables in the subsurface domain. Seismic data for each shot and receiver in the time-domain is in this way matched with subsurface layers, and assimilated in a sequential manner. The methodology is showcased on a realistic model of the Gullfaks field in the North Sea, where we study effects of various seismic acquisition design set-ups, algorithm and model parameter settings. We investigate how these acquisition designs and parameters influence the uncertainty reduction and bias of the inversion results. We highlight the importance of studying statistical performance metrics to ensure a balance between bias and underestimation of uncertainty.

SummaryModel-based information about the global water cycle, in particular the redistribution of terrestrial water masses, is highly relevant for the understanding of Earth system dynamics. In many geodetic applications, hydrological model results play an important role by augmenting observations with a higher spatio-temporal resolution and gapless coverage. Here we demonstrate the feasibility of the high-resolution, open-source hydrological model OS LISFLOOD to simulate terrestrial water storage (TWS) variations with a spatial sampling of up to about 5 km (0.05○). Validation against data from satellite gravimetry reveals that the choice of the maximum soil depth has a significant impact on long-term trends in TWS, mainly in the deepest soil layer. We find that refining the soil depth definition effectively reduces spurious TWS trends, while preserving accuracy in modeled river discharge. Using the modified model set-up, we show that in many regions TWS from OS LISFLOOD fits better to observations than TWS from the Land Surface Discharge Model (LSDM) routinely operated at the GFZ and used in geodetic applications worldwide. The advantage of the high spatial resolution of the OS LISFLOOD implementation is shown by comparing vertical surface displacements to GNSS observations in a global network of stations. The data set presented here is the first application of OS LISFLOOD to generate quasi-global (regions south of 60○S excluded) daily 0.05○ TWS fields for a 23-year period (2000–2022).

SummaryAnalysis of the 48-year of satellite laser ranging (SLR) data of multiple satellites shows that long-term variations in the Earth's dynamical oblateness represented by the second-degree zonal harmonic J2 is best characterized by the superposition of a quadratic trend, 10.5-year and 18.6-year variations. These variations result from climate-related mass changes, tides, and core flow-induced variations at the core-mantle boundary (CMB). We determined that the global ocean's response to the lunar attraction at the 18.6-year period is near equilibrium, with an amplitude of 0.4735 ± 0.008 cm and an error of ∼11% relative to the modeled amplitude (0.4224 cm), and ∼2±3 degrees of phase lag. The 18.6-year Love Number was found to be 0.01375-i0.00553 with an error of 2% for both the real and imaginary parts and a phase correction π for the imaginary part of the IERS 2010 anelasticity model. The nominal frequency-independent anelasticity Love number, k₂, was determined to be 0.3022 ± 0.0001 for the 18.6-year period, based on a reference frequency of 200 seconds and α = 0.1514 for mantle anelasticity for mantle anelasticity. This study also reveals a significant gravitational signal (3.06×10−11) in J2 obstructs the Earth's mantle anelastic response to the 18.6-year tidal forces, reflecting in phase shift of π the imaginary part of the IERS 2010 Love number. This signal can be characterized by a positive Love number of 0.01106 in the modeling of the variation in J2 coupling with the 18.6-year tide. This signal is possibly produced by the core dynamics, which creates a gravitational signal in J2 with an amplitude of 3.36×10−11 at the decadal time scale and could account for ∼70% of the observed 10.5-year variation.

SummaryTo clarify the 3-D crustal and upper mantle structure of the Bolivian orocline in the Central Andes, we conduct azimuthal anisotropy tomography using newly measured teleseismic fundamental mode Rayleigh-wave phase and amplitude data at periods of 25-110 s. Our tomography shows that the subducting Nazca slab is imaged as a high-velocity zone beneath the study region, except for areas where the Nazca ridge is subducting. Azimuthal anisotropy in the subducting slab generally exhibits trench-parallel fast-velocity directions (FVDs), but it becomes complex in and around the subducting Nazca ridge. Low-velocity anomalies with trench-normal FVDs exist in the mantle wedge beneath active arc volcanoes and backarc regions beneath Altiplano. A significant high-velocity zone with relative weak anisotropy exists in the crust of the overriding plate above the Peruvian flat slab in the study region. In contrast, low-velocity anomalies with trench-parallel FVDs are revealed in the crust beneath Altiplano. Furthermore, a high-velocity zone with depth-varying FVDs appears beneath Eastern Cordillera and its surrounding regions, which may indicate the westward underthrusting cratonic lithosphere. These tomographic features well capture the primary 3-D structure of the middle-lower crust and upper mantle beneath the Bolivian orocline, which results from the subduction of an oceanic lithosphere and the delamination and underthrust of a continental lithosphere, leading to the second-highest plateau on Earth.

SummaryIn this study, we applied the “in-situ Vp/Vs method” to monitor variations of seismic velocity ratio (Vp/Vs) within swarms, providing insights into eruption processes. This method, particularly effective in volcanic regions, estimates Vp/Vs by comparing P- and S-wave arrival times of closely located earthquake pairs, reducing errors from unknown crustal velocity variations and is well-suited for detecting rapid changes associated with volcanic swarms. Our study focused on seismic swarms on the Reykjanes Peninsula, south-west Iceland where, swarms have been frequent since 2017 and led up to eruptions in 2021, 2022, and 2023. We analyzed the entire period (2017–2023) as well as the 2021 swarm separately using data from over 40,000 seismic events recorded by the REYKJANET network. We observed significant decrease in the Vp/Vs ratio before major pre-eruption swarms, compared to the background Vp/Vs value of 1.78. From the 2020 swarm, we observed a lower Vp/Vs of 1.72, but the lowest estimated value was 1.70, associated with the 2021 pre-eruption swarm that preceded Fagradalsfjall's first eruption after 7000 years. Reduced Vp/Vs ratios were also noted before the 2022 and 2023 eruptions, suggesting supercritical fluids in the crust during these stages. We also introduce the concept of “change points” to interpret Vp/Vs variations along the dyke. Change points denote specific locations or times of significant Vp/Vs shifts, potentially indicating subsurface changes such as fluid influx or new fracturing from magma intrusion. Identifying these points allows us to pinpoint key moments when the system undergoes substantial changes, offering insights into eruption timing and location. Focusing on 2021 pre-eruption swarm, interestingly the spatial change point found in a location very close to the eruption site. Temporal analysis identified two main change points: the first corresponding with initial activity in the northern dyke and the second with a shift to the southern segment, ultimately leading to eruption. These points mark stages in magma progression, with each showing an initial rapid Vp/Vs drop that could indicate CO₂-rich fluid infiltration, followed by an increase as magma enters. The in-situ Vp/Vs method's sensitivity to changes in seismic properties makes it a powerful tool for real-time volcanic monitoring. By detecting critical Vp/Vs changes with minimal computational demand, this method has potential for integration with online seismic networks, providing an effective early warning system for volcanic hazards.

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.