Geophysical Journal International

AbstractWe investigate the modeling of P-wave reflections on the mantle transition zone (MTZ) discontinuities (Pv410p* and Pv660p*) using ambient seismic noise generated by distant oceanic sources. Using ray theory and waveform simulations, we assess biases in arrival times and amplitude ratios when interpreting noise correlations as Green‘s functions. Our results show that source distribution and the b-caustic effect strongly influence signal recovery. Simulations based on realistic oceanic models (WAVEWATCH III) demonstrate that appropriate source conditions significantly reduce biases. This approach enables reliable imaging of the MTZ, particularly in regions like the greater Alpine area with favorable microseismic source distribution.

SummaryHigh-resolution detection of hidden geological faults is vital for city planning, earthquake disaster prevention, and large-scale engineering construction. This study deployed 229 short-period seismometers across a 10×30 km region within the Xianlin area of Nanjing. Of which, 199 formed a 2-D array, and 30 formed a linear array. Various methods were applied to detect hidden faults in the study area. Using ambient noise tomography, a three-dimensional (3D) S-wave velocity structure was obtained from the surface to a depth of 6.0 km, allowing the first locations of a hidden fault to be mapped via velocity anomalies. A linear array was subsequently deployed based on these early findings, and the horizontal-to-vertical spectral ratio (HVSR) method was applied to estimate bedrock depth and define shallow fault features in greater detail. Finally, a shallow seismic exploration was performed to verify the detection results of ambient noise tomography and HVSR analysis. The results indicate the presence of a hidden fault in the study area, which manifests as a distinctive area of alteration in the high- and low-velocity anomalies in the 3D S-wave velocity structure. Significant variation was identified in the sediment layer thickness in the shallow subsurface, as observed in the HVSR records. In addition, shallow seismic exploration defined important wave-group phase-axis discontinuities in areas with abrupt sedimentary thickness changes. Thus, the hidden fault identified in this study is a normal fault with a nearly north-dipping direction, dip angle of approximately 60°, and fault displacement of approximately 30 m. By linking these results with previous data, it is possible to suggest that such hidden faults are part of the Mufushan–Jiaoshan Fault. Future urban designs and buildings must thoroughly consider the seismic dangers in this region and apply suitable mitigation strategies.

SummaryShipborne gravity anomaly data exhibit multi-vintage characteristics due to their extended temporal coverage. Currently, the measurement accuracy of gravimeters and the processing methods for shipborne gravity anomaly data have been significantly improved and refined. At this stage, the influence of temporal error on the processing of shipborne gravity anomaly data has become an issue that cannot be neglected. We propose a joint reprocessing method for multi-vintage shipborne gravity anomaly data considering temporal error effects. Firstly, the gross error of the shipborne gravity anomaly data is eliminated and filtered. When compensating for the survey line error, the time variable is added to the error equation in order to retain the temporal information in the observed value. The corrected shipborne gravity anomaly data by this method is closer to the real gravity field information. We applied this method to the real shipborne gravity anomaly data in the Philippine Sea. The results showed that the standard deviation of the discrepancy at the intersection points of the survey lines was reduced from the initial 13.46 mGal to 4.30 mGal. The shipborne gravity anomaly data processed after considering the temporal error effects conforms more closely to the actual gravity field information.

SummaryPassive surface wave method is increasingly being applied to urban subsurface exploration due to its non-invasiveness, low cost, and high efficiency. However, its imaging quality is often influenced by limited data acquisition time and the heterogeneous distribution of seismic ambient fields in complex urban environments. To extract coherent surface wave signals for seismic imaging in such challenging setting, we developed a multi-stage urban ambient noise deep clustering framework based on a convolutional autoencoder and deep embedded clustering algorithm. The initial clustering characterizes the distribution patterns of urban noise sources, which informs a secondary, finer clustering to select noise sources optimized for urban seismic imaging. Real-world experiment on the urban train noise field demonstrates our urban noise cluster framework effectively identifies and elucidates the temporal evolution patterns of moving train sources. Compared to traditional data selection methods, our approach yields superior dispersion measurements and significantly attenuates artifacts from the fundamental mode. Furthermore, by employing mode-specific clustering, we successfully capture the refined first overtone, enhancing the accuracy and depth resolution of seismic imaging. This study presents a new perspective to analyzing and selecting complex noise sources, significantly advancing seismic imaging and monitoring in alignment with emerging Artificial Intelligence trends.

SummaryBoth short-term coseismic off-fault damage and long-term fault growth during interseismic periods have been suggested to contribute to the formation and evolution of fault damage zones. Most previous numerical models focus on simulating either off-fault damage in a single earthquake or off-fault plasticity in seismic cycles ignoring changes of elastic moduli. Here we developed a new method to simulate the damage evolution of fault zones and dynamic earthquake cycles together in a 2D anti-plane model. We assume fault slip is governed by the laboratory-derived rate-and-state friction law while the constitutive response of adjacent off-fault material is controlled by a simplified version of the Lyakhovsky-Ben-Zion continuum brittle damage model. This study aims to present this newly developed modeling framework which opens a window to simulate the co-evolution of earthquakes and fault damage zones. We also demonstrate one example application of the modeling framework. The example simulation generates coseismic velocity drop as evidenced by seismological observations and a long-term shallow slip deficit. In addition, the coseismic slip near the surface is smaller due to off-fault inelastic deformation and results in a larger coseismic slip deficit. Here we refer to off-fault damage as both rigidity reduction and inelastic deformation of the off-fault medium. We find off-fault damage in our example simulation mainly occurs during earthquakes and concentrates at shallow depths as a flower structure, in which a distributed damage area surrounds a localized, highly damaged inner core. With the experimentally based logarithmic healing law, coseismic off-fault rigidity reduction cannot heal fully and permanently accumulates over multiple seismic cycles. The fault zone width and rigidity eventually saturate at long cumulative slip, reaching a mature state without further change.

SummaryGreen’s function expressions for seismic interferometry in acoustic and elastic media have been extensively studied and applied across a wide range of applications, including surface-wave tomography and generating virtual shot gathers. However, analogous expressions for coupled acoustic-elastic media systems remain absent, despite their importance for analysing cross-correlation wavefields from ocean-bottom nodal and seismometer recordings and other seismic problems in marine settings. To address this issue, we derive convolution- and correlation-type reciprocity relations for physically coupled acoustic-elastic media by combining Rayleigh’s and Rayleigh-Betti reciprocity theorems, incorporating the constitutive equations governing coupling at the acoustic-elastic interface, and applying time-reversal invariance principles for an arbitrary 3-D inhomogeneous, lossless medium. The derived relationships show that the acoustic and elastic Green’s functions between any two observation points in the medium can be expressed as integrals of cross-correlations of wavefield observations at those locations, generated by sources distributed over an arbitrarily shaped closed surface enclosing the two observation points. When the Earth’s free surface coincides with the enclosing surface, integral evaluation is required only over the remaining portion of the closed surface. If the sources are mutually uncorrelated ambient sources, the Green’s function representation simplifies to a direct cross-correlation of wavefield observations at the two points, generated by a specific ambient source distribution on the closed surface. However, in practical scenarios, the ideal source distribution necessary to retrieve Green’s functions is rarely realized, for example, due to non-uniform illumination. To address these challenges, we represent the ambient cross-correlations as self-consistent observations and introduce a cross-correlation modelling methodology that accounts for practical limitations in source distribution for coupled acoustic-elastic media scenarios. We illustrate the theory by modelling ambient cross-correlation wavefields for a deep-water scenario.

SummarySeismic signals generated by near-surface explosions, with sources including industrial accidents and terrorism, are often analysed to assist post-detonation forensic characterisation efforts such as estimating explosive yield. Explosively generated seismic displacements are a function of, amongst other factors: the source-to-receiver distance, the explosive yield, the height-of-burst or depth-of-burial of the source and the geological material at the detonation site. Recent experiments in the United States, focusing on ground motion recordings at distances of <15 km from explosive trials, have resulted in empirical models for predicting P-wave displacements generated by explosions in and above hard rock (granite, limestone), dry alluvium, and water. To extend these models to include sources within and above saturated sediments we conducted eight explosions at Foulness, Essex, UK, where ∼150 m thicknesses of alluvium and clay overlie chalk. These shots, named the Foulness Seismoacoustic Coupling Trials (FSCT), had charge masses of 10 and 100 kg TNT equivalent and were emplaced between 2.3 m below and 1.4 m above the ground surface. Initial P-wave displacements, recorded between 150 and 7000 m from the explosions, exhibit amplitude variations as a function of distance that depart from a single power-law decay relationship. The layered geology at Foulness causes the propagation path that generates the initial P-wave to change as the distance from the source increases, with each path exhibiting different amplitude decay rates as a function of distance. At distances up to 300 m from the source the first arrival is associated with direct propagation through the upper sediments, while beyond 1000 m the initial P-waves are refracted returns from deeper structure. At intermediate distances constructive interference occurs between P-waves propagating through the upper sediments and those returning from velocity-depth gradients at depths between 100 and 300 m. This generates an increase in displacement amplitude, with a maximum at ∼800 m from the source. Numerical waveform modelling indicates that observations of the amplitude variations is in part the consequence of high P- to S-wave velocity ratios within the upper 150 m of saturated sediment, resulting in temporal separation of the P- and S-arrivals. We extend a recently developed empirical model formulation to allow for such distance-dependent amplitude variations. Changes in explosive height-of-burst within and above the saturated sediments at Foulness result in large P-wave amplitude variations. FSCT surface explosions exhibit P-wave displacement amplitudes that are a factor of 22 smaller than coupled explosions at depth, compared to factors of 2.3 and 7.6 reported for dry alluvium and granite respectively.

SummaryIntermittent fluid injection aims at inferring and steering hydraulic transmissivity and has become an integral part of reservoir stimulation techniques. Modeling the poroelastic response of such pumping operations poses new challenges with respect to the hydromechanical coupling. This is because when a fluid pressure perturbation is introduced in the pore space of a deformable porous rock, it will induce a stress perturbation in the solid phase and this is accompanied by pore boundary motion. Within the limits of quasi-static linear poroelasticity, we analyze the macroscopic signatures of pore boundary motion during injection, i.e., when the rock frame is mechanically loaded, and after injection stop, i.e., when pore boundaries tend to relax back into equilibrium. We show that there is a pumping sequence that allows to harness the energy associated with pore boundary motion accumulated during the frame-loading cycle. Our results foster the need to distinguish how pressure diffusion in poroelastic solids proceeds: either fluid transport is of compressible or incompressible nature and the respective diffusion constant depends on undrained or drained poroelastic moduli.

SummaryThermal properties such as thermal conductivity (TC), thermal diffusivity (TD), and specific heat capacity (SHC) are essential for understanding and modelling the subsurface thermal field. In sedimentary basins, these parameters play a key role in characterizing the present-day thermal state and predicting its evolution, for example, in response to future geo-energy utilizations. Given the wide range of potential geo-energy utilizations and the frequent lack of sufficient sample material, many studies have focused on developing accurate prediction approaches. Machine learning (ML) offers promising non-linear statistical methods to enhance the mapping of interrelations between standard geophysical well-log readings and thermal rock properties. In this study, we introduce an open-access tool for computing profiles of thermal rock properties from standard geophysical borehole logging data, building upon and extending previous petrophysical studies. The tool employs various machine-learning approaches trained on large, physically modelled synthetic datasets that account for mineralogical and porosity variability across major sedimentary rock groups (clastic rocks, carbonates, and evaporates). It establishes functional relationships between thermal properties and different combinations of standard well-log data, including sonic velocity, neutron porosity, bulk density, and the gamma-ray index. We trained four different models including linear regression, AdaBoost, Random Forest, and XGBoost using 80 per cent of the synthetic group data for model development, including training and hyperparameter tuning through cross-validation. The remaining 20 per cent was held out as an independent test set for statistical validation, feature recognition, and input variable importance analysis. A total of 15 input log combinations (including all one, two, three, and four well-log configurations) were evaluated across four machine learning models (linear regression, AdaBoost, Random Forest, and XGBoost), resulting in 180 trained models. The model's predictive accuracy and reliability were further evaluated against independent laboratory drill-core measurements reported in recent studies. Our results indicate that the best-performing predictive models vary depending on the available log-combinations. However, XGBoost frequently outperforms other models in sedimentary rocks. When at least two well logs are provided as input variables, the best-performing models predict thermal conductivity with an uncertainty below 10 per cent relative to borehole validation data (with laboratory-measured thermal conductivity). In most tested model cases and for most input log combinations, predictive errors for thermal conductivity range between 10 and 30 per cent at the (point measurement) sample scale (cm to half a meter). However, when averaged over geological formations or borehole intervals (tens to thousands of meters), the accuracy of thermal-conductivity predictions improves significantly, with uncertainties of the interval mean conductivity dropping below 5 per cent for large intervals. For specific heat capacity, prediction accuracy for the best-performing models at the measurement scale is typically better than 5 per cent. Thermal diffusivity reflects a larger variation, accumulating the uncertainties from conductivity and heat capacity. The presented log-based Python prediction tool provides an automated means to compute thermal parameters using the most suitable ML model for given well-log inputs, facilitating enhanced thermal characterization in sedimentary settings. This has practical relevance for geothermal or hydrocarbon exploration, or subsurface storage projects.

SummaryThis paper investigates the seismic activity and velocity structure in the Three Gorges Reservoir (TGR) region using high-quality travel time data from an extensive seismic observation network. The primary goal is to understand the relationship between the three-dimensional velocity structure and seismicity within the reservoir area. We employed advanced inversion techniques to develop detailed 3-D models of the P- and S-wave velocities and analyzed the focal mechanisms of significant seismic events. Our results reveal that there are substantial lateral variations in the upper crustal velocity structure, with high-velocity zones in the northeastern region of Badong and lower velocities in the Zigui Basin (ZGB). The sedimentary layers in the ZGB are 6–8 km thick, and low S-wave velocity anomalies extend from this depth and are correlated with the Triassic formations. The seismic activity patterns show that the earthquakes in the Badong region were concentrated along three east–west trending belts within the core of an anticline. These patterns suggest that the geological structures and fluid infiltration significantly influence the seismicity. In particular, the M5.1 Badong earthquake occurred at the boundary of a high-velocity zone and was associated with a seismic belt extending from shallow to deeper depths. The results of this study highlight the complex interactions between rock heterogeneity, fault dynamics, and fluid effects, providing a comprehensive analysis of reservoir-induced seismicity. This work provides a better understanding of the physical mechanisms driving seismic activity in large reservoir systems and provides insights relevant to seismic hazard assessment and reservoir management.

SummaryCharacterizing ore deposits or mining dumps in terms of mineral content and grain size remains a challenge. Since the 1950s the Induced Polarization (IP) method has been successfully applied in ore prospecting. However, reliably interpreting field survey data requires comprehensive laboratory studies to establish a link between the IP parameters from empirical or phenomenological models and the type and quantity of ore minerals. In this study, we use numerical electrical networks to replicate the complex electrical resistivity spectra observed in experiments on sand-pyrite-water mixtures. A network consists of a 3D assembly of resistors, representing the saturated pore space, and leaky capacitors simulating the electrical behaviour of ore minerals. A sophisticated fitting procedure enables the precise determination of resistor and capacitor parameters, ultimately leading to strong agreement between measured and synthetic IP spectra. The results obtained from the 3D network align well with the classical Pelton model, which is based on a simple equivalent circuit. Our findings indicate that the network's chargeability depends on the fraction of capacitors in the system (i.e. the number of capacitors divided by the number of capacitors and resistors), and that the Pelton time constant of the measured spectra is closely related to the resistor and capacitor parameters. We argue that a 3D approach offers a more realistic framework, paving the way for future studies on the effects of ore grain size distribution, and the spatial arrangement of ore grains.

SummaryThe most widely used method to derive global geomagnetic field models for historical and longer timescales has long been regularized non-linear least squares inversion. It is based on spherical harmonics for the spatial part and cubic B-splines for the temporal dynamics. Recently, different versions of Bayesian inversion have been applied for this purpose. Early literature on the traditionally used formalism states the inverse problem in a Bayesian setting and discusses uncertainty estimation via the posterior covariance, but this view was lost in subsequent studies in the geomagnetic community. Here we aim to provide both geomagnetic field modellers and users of such models with a comparative view of the methods to enable them to better evaluate strengths and weaknesses of different models. We first describe the connection between regularized least squares and Bayesian inversion in general form in a linear, one-dimensional setting. A fully Bayesian perspective allows interpreting the regularization term as a form of prior and offers new ways of comparing models from both approaches. We then discuss the particular case of geomagnetic field modelling. We find that in comparison to Bayesian modelling approaches the prior corresponding to the widely used regularization does not imply reasonable field properties and does not lead to meaningful uncertainty estimates.

AbstractThis study presents a refined interpretation of the Jeokjung-Chogye Basin (JCB), a confirmed meteorite impact structure in South Korea, by integrating high-resolution gravity data, microtremor measurements, and borehole information. A total of 1 700 gravity stations including 1 000 newly acquired in 2023, were used alongside horizontal-to-vertical spectral ratio (HVSR) analysis and well-log constraints to characterize subsurface structures. To isolate the impact-induced deformation from overlying sedimentary effects, gravity stripping was applied to remove the signal from post-impact ejecta deposits. The residual gravity field was analysed using dip-curvature mapping and Euler deconvolution, which revealed concentric ring structures with displaced centres. These asymmetries, corroborated by 3D forward gravity modelling using IGMAS+, suggest a northeast-to-southwest impact trajectory with an oblique incidence angle of approximately 45°, contrasting with earlier estimates of ~55° from east to west. The final 3D density model achieves a strong correlation with observed anomalies (R ≈ 0.95) and successfully resolves variations in the autochthonous and basement layers.

SummaryThe mantle transition zone (MTZ) plays an important role in the global material circulation, slab dynamics, and seismogenesis of deep earthquakes in subduction zones. Here we construct fine MTZ structures of both P and SH waves beneath the northeast Asia continental margin using improved grid-search waveform modelings, based on high-quality triplicated waveforms of three deep Kuril earthquakes recorded by the China Digital Seismograph Network (CDSN). We find a high-velocity anomaly (HVA, average 3.3 per cent δVp and 2.3 per cent δVs) with a thickness of 130-138 km in the lower MTZ. The HVA hosts a top interface with positive velocity contrasts (δVp: 4.3 per cent, δVs:3.2 per cent), while the 660-km discontinuity (660) shows reduced velocity contrasts (δVp: 3.6 per cent, δVs: 5.1 per cent) and negligible depressions of less than 10 km. The HVA we detected likely implies the thickened stagnant Pacific slab that may alter localized heat exchanges between the MTZ and lower mantle. The increased Vp/Vs ratio (∼1.85) indicates a water-rich state (∼0.42 wt per cent) inside the stagnant slab, evidencing the deep water transportation by the slab subduction. We infer that the interior localized dehydration of hydrous minerals within the stagnant slab may trigger the large outboard 1990 Mw 7.2 Sakhalin (Kuye) earthquake. Our results can provide more insight into slab dynamics and seismogenesis of deep earthquakes in northeast Asia.

SummaryWe used broadband ocean bottom seismometer (BBOBS) data from the RHUM-RUM experiment to derive the compliance function and estimate the shear velocity (Vs) structure of the subsurface at several sites beneath the Indian Ocean. The primary objective is to map the geological features of poorly explored marine regions, utilizing the compliance function, a measure of seafloor deformation in response to infragravity pressure signals at low frequencies (0.003 to 0.04 Hz). Compliance is the transfer function between vertical displacement and pressure, which is most sensitive to subsurface shear velocities. Our analytical process involves several data processing steps, including the removal of glitches, filtering out seismic events, minimizing tilt effects, calibrating pressure gauges, searching over the frequency and coherence domains to determine the optimal data window, and performing depth-velocity inversion using Monte Carlo method, specifically the Metropolis-Hastings algorithm. We present the ’ComPy’ software, which automates these processing steps for seafloor compliance analysis. The data, recorded over 13 months in 2012-2013 over a large region stretching from La Reunion Island to the Central Indian Ridge (CIR) and the South-West Indian Ridge (SWIR) (water depths of 3 to 5 km), confirm the stability of the compliance function over time. Depth-velocity inversions of the derived compliance measurements, using the Metropolis-Hastings algorithm, illuminate the Vs structure of the oceanic crust down to 8 km. Low Vs anomalies in the crust at the SWIR are consistent with significant serpentinization of a crustal component of tectonically exhumed mantle-derived peridotites.

SummaryIn this study, we present a new algorithm and accompanying software for 3D gravity inversion of density structures. The algorithm combines the strengths of the Bayesian approach-which incorporates prior model information through a variogram model-with the advantages of the Tikhonov regularization framework to address the challenge of depth resolution. We also provide a detailed derivation of the procedure for calculating and fitting the 3D experimental variogram, which serves as a fundamental input to the algorithm. The software implementing the proposed algorithm was developed using the widely adopted computational programming language Matlab. To evaluate its effectiveness, we conducted four representative experiments, ranging from simple to complex scenarios. The synthetic results demonstrate that incorporating model covariance constraints yields a more localized and better-focused density distribution compared to results obtained without such constraints. Additionally, we tested the algorithm's robustness by introducing noise into the observation data. The results show that the proposed method is resistant to noise and maintains strong performance. Finally, we applied the algorithm and software to real field data and compared the results with those from previous studies. The comparison confirms that our method is capable of producing reliable, high-resolution 3D density models, with the added advantage of integrating prior information.

SummaryReliable earthquake forecasts depend on modellers making suitable choices regarding data selection and model implementation. We explore some of these choices in order to construct short-term forecasts for independent background events in Italy, using a Log-Gaussian Cox Process to describe the spatial intensity of seismicity as a function of physical spatial covariates and a random field. We explore correlations between a range of physical covariates available for Italy, and investigate how they might contribute to observed spatial seismicity patterns. We find that historic smoothed seismicity, strain rates and elevation are most useful in describing instrumentally-observed seismicity. We use point process intensity models constructed with combinations of these covariates with past seismicity to generate Bayesian earthquake forecasts, constructing simulated catalogue forecasts of future earthquake catalogues in Italy, discussing modelling choices we make along the way. Finally, we test the performance of these forecasts in a pseudo-prospective manner for the 2010-2021 period and using historical data to assess how well they might perform for extreme events. While the inclusion of historic events improves forecasts, models without historic seismicity also perform well, indicating a degree of stationarity in the process. We demonstrate that combinations of covariates can add predictive power over modern catalogue data alone, and provide a framework for future modelling.

SummaryThe rapid development of the global economy brings frequent artificial seismic events related to industrial activities. The identification of faults revealed by natural earthquakes, along with the temporal and spatial distribution of b-values derived from these seismic events, is influenced by their occurrence. Correcting identification of seismic events has significant influence on seismic risk assessment. In this paper, we present a tri-branch convolutional neural network (CNN) model, known as EQTypeNet, to automatically classify three classes of events (natural earthquakes, explosions, and collapses) in China, using waveform, spectrogram, and amplitude ratio features as inputs. The EQTypeNet has developed in two steps: 1. Binary classification: Mitigating data imbalance by supplementing non-natural events outside the test set area. In the classification of explosions and earthquakes, we obtained an individual station macroF1 of 0.99 and a classification accuracy of 98.7 per cent. In the classification of collapses and earthquakes, we obtained an individual station macroF1 of 0.99 and a classification accuracy of 99.5 per cent. 2. Ternary classification of natural earthquakes, explosions, and collapses: Enhancement of the model's classification in small areas by transfer learning. In the transfer learning of Shanxi and Northeast China, the F1-score of individual stations reached 0.98 and 0.97, with corresponding accuracies of 98.6 per cent and 96.8 per cent. The EQTypeNet applies to different magnitude ranges for seismic event discrimination and shows a satisfactory application performance in a large area. It will be expected to achieve higher performance in the ternary classification by further supplementing the training with suitable physical features.

SummaryIn February 2025, a strong seismic crisis occurred 35 km northeast of Santorini, an active volcano located in South Cyclades, Greece, a region of distributed extensional active faulting. The GNSS data shows an inflation of the volcano since August 2024. We model it with a magma source of 7.7 106 m3 located 3.1 km under the north-central floor of the caldera, near the inflation centre of 2011–2012. After 24 January 2025, the seismic activity, until then localised within the caldera, shifted offshore Santorini and increased with eight Mw ≥ 5 events and ground motion exceeding 5 mm at Syros at 110 km from the epicentres. The GNSS data is consistent with a model of dislocation involving a south-east dipping normal fault located between the Kolumbo submarine volcano and the Anydros islet, 18 km long, 12 km wide, with a tip at 7.5 km depth and ~3 m of slip. The ~17.5 1018 Nm corresponding seismic moment, much greater than the ~1.5 1018 Nm of the recorded earthquakes, reveals the probable occurrence of a slow-slip earthquake of equivalent magnitude 6.8. This event might have been triggered by a small dyke injected from the inferred source beneath Santorini. However, the subsidence recorded at Santorini and Anydros is incompatible with the hypothesis of a large dyke injected beneath Kolumbo-Anydros.

SummarySeasonal variations in Nepalese seismicity have been reported with varying degrees of confidence. We re-investigate these claims by analysing 20 years of Nepalese seismicity before the 2015 Gorkha earthquake, as detected by the Nepalese national network, and focusing on earthquakes located along the eastern and central sections of the Nepalese Main Himalayan Thrust. Using several declustering techniques, we find no statistically robust evidence of seasonal seismicity in the studied record, regardless of magnitude threshold above completeness. This suggests that previously reported seasonality may be restricted to the western section of the Nepalese orogeny, may be an artefact, or may indicate that nucleation times of earthquakes are longer than the year. We also investigate potential annual variations in the Gutenberg-Richter b-value, given its recent observed modulation by transient stressing. Additionally, we use large-scale mass redistribution estimated from the monthly gravity field retrieved from the Gravity Recovery And Climate Experiment and Follow-On (GRACE/-FO) missions, to resolve stress variations at depth induced by transient surface loads. We find that the mean annual b-value peaks when seasonal Coulomb stress rates reach their minimum value at the height of the summer rainy season. When considering the combined effect of tectonic and seasonal loading, this corresponds to a recurring period of stress reversal, when Coulomb stress momentarily decreases. This suggests that periodic clamping of the Main Himalayan Thrust reduces the likelihood of earthquakes growing to larger magnitudes in accordance with hierarchical rupture models. The susceptibility of b-value to stress variations of roughly 0.1 points.kPa−1 is consistent with recent estimates of b-value sensitivity to transient loading, although it remains high when compared to the stress-dependency associated with both static differential stress, and with long-term evolution during the seismic cycle. This discrepancy points to the large impact of stress transients on the dynamics of seismic rupture.