Geophysical Journal International

SummaryFirst-arrival traveltime tomography (FATT) is a widely used method for characterizing near-surface velocity structures in geotechnical engineering and resource exploration. We introduce an improved version of FATT, named first-arrival adjoint tomography (FAAT), which involves the joint inversion of first-arrival absolute traveltimes and differential traveltimes. Unlike absolute traveltimes, differential traveltimes, derived from common sources or receivers, offer heightened sensitivity to fine-scale structures near the receivers or sources, respectively. This dual sensitivity makes FAAT particularly effective in imaging highly heterogeneous media. However, the proximity of rays associated with differential traveltimes can lead to the instability of the inversion. To overcome this challenge, we simultaneously incorporate absolute traveltimes and differential traveltimes to update the velocity model. This approach improves the stability of the inversion process, leading to improved resolution of inverted results. We specifically employ the fast-sweeping method to solve the factored eikonal equation, providing robust solutions in models with complex geological structures. Furthermore, we address the inverse problem by computing the gradient of data misfit using the efficient adjoint-state method. Through numerical testing, we validate the effectiveness of FAAT in comparison to that using only absolute or differential traveltimes. Finally, we apply the proposed FAAT method to near-surface site characterization at Bishan-AMK Park in Singapore. Compared with FATT and validated against borehole data, FAAT demonstrates its ability to establish more reliable velocity models, revealing finer details and substantially improving geological interpretation.

SUMMARYNumerical models are a powerful tool for understanding and predicting the impact of landslide-generated tsunamis. We use here the HySEA code which incorporates a multilayer structure and non-hydrostatic pressure to simulate the tsunami generated by a potential submarine landslide located offshore of the Mayotte island. The island is surrounded by a lagoon and steep slopes weakened by the 2018 seismo-volcanic crisis. The influence of the input parameters and of the model assumptions is shown to change by a factor 2 the predicted maximum water free surface elevation, velocity and maximum inundation depth. This demonstrates the need of using numerical models for building local scaling laws to relate tsunami and landslide properties. Our results highlight the necessity of incorporating high-resolution bathymetry, in-depth variations through multilayer modeling and relevant landslide rheology to accurately predict tsunami impact. In case of strong topography variations as in Mayotte, using 4 layers seems to be a good compromise between accuracy and computational cost. Accounting for these effects would enable to refine hazard maps by identifying safe and high-risk coastal zones and to improve wave arrival time estimates, thus reducing tsunami-related risks in regions like Mayotte.

SummaryUnder compression in the brittle regime, rocks fail as fracture interaction, propagation and coalescence produce a process zone of high strains that forms the nucleus of the fault that eventually traverses the rock core. Early analyses suggested that this process zone extends through relatively intact rock close to the peak stress, and that the deformation that develops earlier in loading does not guide the location of this process zone. Here, we assess the predictability of the location of the process zone in triaxial compression experiments on Darley Dale sandstone. We develop supervised machine learning models to predict whether a particular location in the rock core will eventually host high shear and dilative strain, and thus eventually reside within the process zone. We calculate the local incremental strain tensors throughout the rock cores using digital volume correlation of X-ray tomograms captured during the experiments. We find that the models produce accuracies of 0.73 on average, indicating that it is possible to predict the location of the process zone using observations of the strain field preceding macroscopic failure. The average accuracy increases to 0.84 when the data are restricted to later in the experiment because the strain field becomes more similar to the final strain field that hosts the process zone. The models primarily rely on the shear strain field to predict the location of the process zone, likely because it is more spatially persistent than the dilative strain. The local neighborhood of the strain field is helpful for predicting whether a location will eventually host the process zone only when the data are restricted to one experiment. When the data includes several experiments, the models primarily rely on global statistics of the strain fields. The varying correlation lengths of the dilative and shear strain field of different experiments help explain this result.

AbstractActive tectonic movements and surface deformation are observed in the eastern Tibetan Plateau. Understanding variations in crustal thickness and the deep Moho interface is crucial for elucidating the expansion of the Tibetan Plateau. This study utilizes InSAR to derive vertical surface deformation and applies loading corrections based on Green's function method. Additionally, satellite gravity data are used to separate hydrological and tectonic signals to infer changes in the Moho interface. Our results indicate that the regionally averaged loading effects, estimated using localized Green's functions (LGF), account for approximately 16.5 per cent of the InSAR-derived vertical displacement field. This contribution exhibits significant spatial variability, exceeding 100 per cent in regions with strong hydrological activity. The loading calculation is highly sensitive to Earth model: the relative difference between the load displacement obtained using the local Green's function and that obtained with an average Green's function reaches 48.2 per cent. After applying loading corrections, a more accurate Moho uplift rate of –8.2 ± 3.1 mm/a is obtained. The findings support the conclusion that the Moho interface rises in the southern region of the study area, with crustal thinning, while the Moho surface sinks in the northern region, with crustal thickening.

AbstractIn recent years, offshore tsunami observation networks equipped with ocean bottom pressure gauges (OBPGs), such as S-net, DONET, and N-net, have been deployed around Japan, enabling real-time collection of high-quality tsunami data near the source. These networks make it possible to estimate the spatiotemporal variation of the tsunami wavefield using a data assimilation approach, and to predict coastal tsunamis from the initial or current tsunami wavefield. This study proposes a novel tsunami data assimilation method that uses physics-informed neural networks (PINNs) to estimate tsunami wavefields from the observed OBPG data. The neural network was optimised by minimising the sum of the data loss, which quantifies discrepancies from the tsunami data, and the physical loss, which quantifies the satisfaction of the linear long wave equation. This was performed to ensure that the estimated results are consistent with both the observed data and the physics of tsunami propagation, even when there are limited observational data and significant noise. We first validated the effectiveness of the proposed method using synthetic S-net OBPG data from the 2011 Tohoku-oki earthquake (Mw 9.0) tsunami. The results confirmed that by using both data and physical constraints in the PINN optimisation, the PINN could adequately assimilate the spatiotemporal distribution of the tsunami wavefield from OBPG data, even for predictions outside the network coverage area. The predicted tsunami waveforms at the coastal stations, computed from the estimated initial wavefield, showed good agreement with the actual waveforms. Next, we conducted an experiment using actual S-net OBPG data from the 2016 Fukushima-oki earthquake (Mw 6.9) tsunami. The initial tsunami source estimated by PINN was in good agreement with other studies based on waveform inversion, although the maximum source amplitude and maximum coastal tsunami heights were underestimated. We also conducted an experiment using N-net OBPG data from the 2024 Hyuganada earthquake (Mw 7.0) tsunami. The PINN could accurately estimate the initial tsunami source, even though the tsunami source of this event was located outside the N-net coverage area. Finally, we have shown that incorporating tsunami observations over time into the iterative optimisation of the PINN model allows for accurate and efficient tsunami data assimilation.

AbstractUnderstanding processes in the Critical Zone requires reliable information about the vadose-zone aquifer, its geometry, and spatial variability. Commonly, such information is obtained from boreholes, yet large areas might render their application prohibitively expensive. Additionally, limited geological a-priori information might bias the interpretation due to lateral geological changes smaller than the borehole sampling scale. The transient electromagnetic method (TEM) has emerged in the last decades as a well-suited method to efficiently investigate the subsurface, as required for many hydrogeological applications. The interpretation of TEM measurements relies mainly on deterministic inversions, offering only a limited insight on the uncertainty of the subsurface model. Uncertainty quantification, however, is essential for integrating TEM results into hydrogeological models. Hence, we propose a combined approach using both deterministic and stochastic inversion of TEM soundings to investigate the uncertainty of shallow (< 40 m) aquifers. Current stochastic approaches for TEM data rely on Markov chain Monte Carlo algorithms, which have to be run from scratch for each individual sounding. Alternatively, machine learning approaches, such as Bayesian Evidential Learning (BEL), can be much faster because they do not require retraining for every new data set. Hence, we investigate, in particular, the application of a single, common prior model space instead of multiple, individual prior model spaces to directly estimate the uncertainty of multiple TEM soundings. To this end, we combine forward modelling routines with the stochastic inversion approach BEL1D and assess our approach using both field data and numerical experiments.

SUMMARYThe Fennoscandian earthquake catalogue (FENCAT) assembles data on the natural seismicity in Fennoscandia, Northern Europe. We present an updated and standardized version of the catalogue originally published in the early 1990s. New instrumental data are recorded by the seismic networks of Denmark, Estonia, Finland, Norway, and Sweden, and analyzed by the Geological Survey of Denmark and Greenland, the Geological Survey of Estonia, the University of Helsinki in Finland, the University of Bergen and the NORSAR research foundation in Norway, and Uppsala University in Sweden. The updated catalogue provides the available earthquake parameters in a brief, user-friendly version: origin time, source coordinates, focal depth, macroseismic data (maximum intensity and radius of the area of perceptibility), up to three observed magnitudes, seismic moment estimate, and a standardized moment-related magnitude, mW(HEL), for each event. The standardized magnitude is defined in this paper and its relation to other magnitude scales is provided. Suspected non-earthquakes (e.g. frost events, explosions, human-induced events) have been removed. The standardized event magnitudes range from mW(HEL) -1.0 to 6.2. To enable the usage of earthquake data in a large variety of seismological, geological and earthquake engineering investigations, the data are not truncated at the low magnitude end.The updated catalogue, FENCAT (2021), contains about 23 000 earthquakes for the period 1467–2021 in an area bounded by 54–75°N latitudes and 0–45°E longitudes. The completeness and quality of the earthquake solutions is best within the areal coverage of the above-mentioned networks.

AbstractWhile modern thermal convection in rocky planets is controlled by a slow solid-state creep flow, the earliest stages of terrestrial planets likely experienced turbulent flow during which their silicate envelope was fully molten, usually called magma ocean. The main parameter separating the two regimes is the Prandtl number (Pr), which is so high for mantle convection to be usually assumed infinite, whereas magma oceans are characterized by Pr on the order of 1. We compared the results of isoviscous convection simulations performed with three codes: (GAIA, TLBM, StreamV). These codes are based on different numerical formulations and were used for modeling convection with Pr ranging from 1 to 1000, while exploring different convection intensity by varying the Rayleigh number (Ra) from 104 to 106. GAIA (Generic Automaton for planetary Interior Analysis) is a Finite Volume fluid flow and energy solver for the Navier-Stokes equations across arbitrary geometries. TLBM (Thermal Lattice Boltzmann Method) solves the mesocale momentum and energy distribution densities for colliding particles on a discrete lattice. StreamV is a Eulerian-Lagrangian Finite Volume code that solves the Navier-Stokes equations under the Boussinesq approximation. The codes are compared over 24 different simulation setups, analogue to the classical Blankenbach infinite Pr benchmark (Blankenbach et al., 1989a), but extending it to finite Pr and to two types of boundary conditions, free-slip and no-slip. We show that the results of the three codes are generally in good agreement, and discuss differences. Finite Pr solutions show a much richer dynamics varying from stable steady-state solutions, to oscillatory and chaotic ones, and converging to infinite Prandtl number solution for increasing values of Pr for larger Ra: Pr ≥ 100 is sufficient for Ra = 105 but Pr ≥ 1000 is required for Ra = 106. Our results offer a robust set of solutions useful for testing future finite Prandtl number convection codes.

SUMMARYOur viewpoint highlights the challenges faced by women in the Induced Polarization (IP) community (and elsewhere), particularly the persistent gender imbalance in scientific workshops. This underrepresentation in leadership roles, presentations, and discussions reflects broader systemic biases in academia. By sharing the experience of the 7th IP workshop, where the organizing team made deliberate efforts to prioritize gender and diversity in organizing committees, recruitment, and session formats, we demonstrate how intentional actions can create a more inclusive, gender-balanced environment. This approach is crucial not only for the IP community but for all research communities. Emphasizing diversity and inclusion fosters a welcoming atmosphere that encourages participation from all individuals, regardless their background. In turn, this enriches the research process by incorporating diverse perspectives, driving innovation, and improving scientific outcomes. We aim to inspire other academic communities to actively promote diversity and inclusivity, recognizing that such efforts not only benefit underrepresented groups but elevate the success and relevance of science as a whole.

SummaryWe present a novel model-based seismic data redatuming method based upon time-domain data-assimilated wavefield (DAW) reconstruction. Seismic redatuming refers to the positioning of sources and/or receivers from the acquisition surface to a virtual datum level along with the computation of the dataset that would have been recorded by the virtual acquisition. Physically, DAWs approximate the true wavefields by least-squares time-reversal extrapolation of the recorded data in a background subsurface model. To this end, DAWs are first formulated as the sum of the background wavefields and an approximation of the scattered wavefields by the unknown contrasts between the true model and the background model. We estimate the scattered fields via the estimation of their volume scattering sources by fitting in a least-squares sense the differences between the recorded data and the simulated data in the background model, namely the restriction at receivers of the scattered fields. These underdetermined linear source problems involve two main steps: a multidimensional deconvolution of the scattered data by the data-domain Hessian, followed by the propagation of the deconvolved scattered data backward in time. Once the source problems have been solved, DAWs follow by solving the wave equation in the background model using the sum of the experimental sources and the scattering sources as extended sources. Finally, the redatumed datasets, whose accuracy depends on the precision of the background medium between the acquisition surface and the datum, are readily obtained by sampling the DAWs at arbitrary datum levels. The primary challenge in computing accurate DAWs lies in the multidimensional deconvolution of the scattered data, which requires solving a data-domain normal equation with preconditioned Krylov-subspace iterative methods, where each iteration requires one forward and one backward simulation. We demonstrate the accuracy of the method and discuss its computational cost with several benchmarks representing various experimental environments (onshore with weathered layers, shallow and deep offshore).

SummaryDistributed acoustic sensing (DAS) has emerged as a potential solution to the sparse instrumentation issue in the world’s oceans. DAS involves repurposing fibre optic cables into dense receivers. The spatial undersampling limits our understanding of fundamental oceanic processes, like ocean dynamics. We use long-term DAS recordings from Svalbard, Norway, over two roughly perpendicular fibre segments to analyse ocean surface gravity wave (OSGW) signals and gain additional insight into their dynamics. This fibre layout allows estimation of the angle of arrival for OSGW generated under different weather conditions, while the long-term recording allows one to study seasonal variations. By investigating different wind directions, we observe two sets of OSGW arrivals: swells generated by distant storms and waves generated by the local winds. The swells consistently originate from the south-west, whereas the wind-forced OSGW follows, more or less, the local wind direction. Moreover, we conduct a detailed analysis of the recorded swell waves by computing their origin time, great circle propagation distance, group velocity, incidence angle, location, and interference pattern. This yields important data that can be used to characterise local and distant Atlantic storms. Only one receiver system has been employed, and more receivers are needed to validate the results obtained here and gain additional insight into the OSGW signals recorded on DAS systems.

SummaryInversion of geophysical data usually exhibits strong non-uniqueness, arising from sparse data coverage, limited number of measurements, inherent nonlinearity of governing physical laws, noise, and other factors. Methods based on Monte Carlo sampling are commonly used to explore the posterior model distributions, but these approaches are computationally demanding. Variational inference (VI) provides an alternative by transforming a high-dimensional sampling problem into an optimization problem, thereby significantly reducing the computational time. However, conventional VI methods, which typically use simple distribution families, like Gaussians, to approximate the posterior, may lack flexibility necessary to capture the complexity of the posterior distributions. Normalizing flows (NFs), a type of deep generative models, address this limitation by transforming a simple initial distribution into a highly complex target distribution through a sequence of invertible and differentiable transformations. In this study, we develop an NF-based VI method and apply it to electromagnetic (EM) data. This approach allows for explicit integration of prior knowledge and reference models into the inversion process. Both synthetic tests and field applications on EM data demonstrate that NF-based inversion effectively recovers the posterior model distribution in a more efficient manner, while providing excellent data fitting performance. Unlike many other machine learning algorithms, NFs do not require a training set, making it highly transferable across various inversion problems with minimal adjustments. The proposed NF-based method offers a more robust and computationally efficient solution to uncertainty quantification and shows great potential to be extended to solve 3-D geophysical Bayesian inversions, a major challenge that the geophysical community has faced for decades.

SummarySoutheast Asia, bordered by significant tectonic plates such as the Indo-Australian, Pacific, and Philippine Sea Plates, is distinguished by its frequent tectonic activity and complex geological structures, making it one of the most dynamically evolving regions worldwide. In this study, we introduce a novel 3D P-wave velocity model of the upper mantle and transition zone in Southeast Asia using regional seismic traveltime tomography based on first-arrival data from the International Seismological Center. We employ an adjoint-state tomography method with normal-vectors independence to accurately invert for 3D velocities using a 1D reference model. Synthetic tests confirm the reliability of our model in delineating features of the subduction zone and the surrounding region. Our inversion results highlight distinct subducted slabs within the subduction zone and a pervasive low-velocity zone beneath Sundaland, which may be associated with lithospheric thermal weakening. Additionally, a mushroom-shaped low-velocity anomaly attributed to the Hainan mantle plume is identified beneath Hainan Island. The low-velocity anomaly observed beneath the western part of the Java Sea may be attributed to the combined effects of Sunda-Java slab subduction, lower-mantle flow through the Sunda Strait, and the influence of the Hainan mantle plume. Notably, beneath the Andaman Sea, we observe an east-west elongation of the northern Sumatra slab, potentially linked to the clockwise rotational opening of the Andaman Sea. Additionally, three potential rifts are identified beneath the subducting Sumatra-Java slab: beneath the Toba Volcano, the Sunda Strait, and the eastern segment of Java Island. Extensive high-velocity anomalies beneath the Philippine Islands and the South China Sea suggest a double-sided subduction process involving the Proto-South China Sea slab.

SummaryInversion of a given geophysical dataset cannot be complete without assessing the resolution and uncertainties associated with the model obtained. However, model appraisal may still be a challenging task from both a theoretical and a computational point of view. To tackle the problems of model estimation and appraisal, we introduce the Subtractive Optimally Localized Averages (SOLA) method to the geophysical electromagnetic community, through the example of linear inversion of induced polarization (IP) data. SOLA is a variant of the Backus-Gilbert method: it is computationally more efficient but also allows one to specify directly the target local averages of the Earth’s properties to be estimated, including their uncertainties. SOLA offers great flexibility in the construction of averaging kernels, via the design of target kernels, and direct control over the propagation of data errors into the local-average estimates. With SOLA we obtain a collection of i) local averages of the ‘true’ Earth model, accompanied with their ii) averaging kernels and iii) uncertainties. We investigated the performance of SOLA for the 2–D tomographic inversion of a field IP data set. The obtained chargeability model compares well with previous studies, and, most importantly, its resolution (the spatial extent of the averaging kernels) and uncertainties can be interrogated. We conclude that SOLA is a promising approach for geophysical-electromagnetic linear(ised) tomographies. In the case of IP inversion, to construct chargeability models and evaluate their robustness.

SummaryFault systems have geometrically complex structures in nature, such as stepovers, bends, branches, and roughness. Many geological and geophysical studies have shown that the geometrical complexity of fault systems in nature decisively influences the initiation, arrest, and recurrence of seismic and aseismic events. However, a vast majority of models of slip dynamics are conducted on planar faults due to algorithmic limitations. We develop a 3D quasi-dynamic slip dynamics model to overcome this restriction. The calculation of the elastic response due to slip is a matrix-vector multiplication in boundary element method, which can be accelerated by using Hierarchical Matrices. The computational complexity is reduced from O(N2) to O(Nlog N), where N is the number of degrees of freedom used. We validate our code with a static crack analytical solution and the SEAS benchmark/validation exercise from Southern California Earthquake Center. We further employ this method on a realistic fault system with complex geometry that was reactivated during the 2023 Kahramanmaraş–Türkiye doublet earthquakes, generating slip sequences that closely match real observations.

SUMMARYThe spherical harmonic coefficient Level-2 products of the Gravity Recovery and Climate Experiment (GRACE) mission are affected by north-south stripe noise. Toward this end, we have developed a new filter named Variational Mode Decomposition spatial (VMDS) filter that transforms the Equivalent Water Height (EWH) map derived from GRACE level-2 product into a one-dimensional sequence, which is then filtered by using variational mode decomposition. This approach overcomes the limitations of the singular spectrum analysis spatial (SSAS) filter, which well performs in the medium-frequency band but omits the high-frequency NSS noise. We thus put the VMDS filter behind the SSAS filter to utilize the good performance of SSAS in the medium-frequency band and thus propose a combined filter termed SV. A closed-loop simulation demonstrates the better ability of SV to suppress NSS noise and preserve signal at the grid scale compared to the SSAS filter. In the real-world scenario, the SV solution achieves a noise level (46.68 mm of EWH) below that for SSAS and DDK7 solutions (53.52 and 53.68 mm of EWH, respectively) over the ocean at low latitudes. Moreover, the well-documented water level of Lake Victoria and the well-modeled coseismic gravity change of the Mw9.2 2004 Sumatra-Andaman earthquake demonstrate that the SV filter efficiently preserves localized mass evolutions while suppressing north-south stripe noise. Such short-wavelength signals usually miss in highly filtered spherical harmonics (e.g. DDK5 and DDK6) solutions or are significantly inconsistent for various mass concentration solutions.

SUMMARYWe propose a stable and efficient method for high-degree regional lithospheric magnetic field modeling based on spherical Slepian functions, achieving significant improvements in the computational efficiency of solving large linear systems by reducing both complexity and memory requirements. This method leverages the orthogonality of Slepian basis functions on regional domains R to improve the stability of regional modeling. The block-diagonal structure of the normal equation matrix and the sparse representation by Slepian functions are simultaneously exploited to achieve a two-stage matrix compression. Additionally, a Bayesian Information Criterion (BIC)-based strategy is introduced to determine the optimal truncation number for the Slepian basis functions, ensuring a balance between computational efficiency and data fidelity. Our method was validated by directly inverting synthetic regional lithospheric magnetic field data generated from a spherical harmonic model up to degree 1050. The modeling process showed a reduced memory requirement of approximately eight orders of magnitude. The high-degree model successfully reconstructed the desirable magnetic field at different heights, and the residual statistical analysis results showed that the spatial variation of the lithospheric magnetic field was accurately captured. In addition, the proposed method can be readily extended to gravity modeling and other applications that utilize spherical harmonic analysis. The matrix compression techniques adopted provide an ideal framework for parallel computing, showing a wide range of application potential.

SUMMARYSeismic ambient noise tomography has been a powerful tool for seismic imaging, but most existing approaches fail to accurately predict detailed correlation waveforms in the presence of spatially heterogeneous noise distributions. Full-waveform ambient noise inversion allows for high-resolution waveform-based inversion even in substantially heterogeneous noise fields. Unfortunately, the computational cost of this approach is constrained by the noise distribution, rather than the inversion domain: global-scale applications are viable, but smaller-scale applications require a modified approach. We present a general finite-domain full-waveform ambient noise inversion methodology, providing approximation mechanisms for treating out-of-domain propagation from distant noise sources. This makes the problem tractable on much smaller domains. In a numerical example, we demonstrate that this approach enables inversion for both structure within a chosen domain and approximate noise source distribution outside it.

AbstractThe British Isles' lithospheric structure, shaped by a dynamic geological history, remains incompletely understood, particularly regarding anelastic parameters, such as attenuation. In this study, we present a teleseismic attenuation model for the British Isles, using time-domain analysis of teleseismic P-wave data from 28 deep earthquakes. We constructed a 2D differential attenuation map (Δt*) that reveals significant regional variations. Our findings show a weak anticorrelation between Δt* and shear wave velocity at upper mantle depths, suggesting that variations in the lithosphere-asthenosphere system influence this pattern. A high-attenuation zone extends from Scotland across the Irish Sea to southwest England, potentially linked to mantle upwelling associated with the Iceland plume. This model provides new insights into the mantle dynamics beneath the British Isles, offering a crucial reference for future geophysical studies in the region.

SUMMARYInduced seismicity poses a significant challenge to the safe and sustainable development of Enhanced Geothermal Systems (EGS). This study explores the application of machine learning (ML) for forecasting cumulative seismic moment (CSM) of induced seismic events to evaluate reservoir stability in response to fluid injections. Using data from the Cooper Basin (Australia), the St1 Helsinki geothermal project (Finland), and a controlled laboratory injection experiment, we evaluate ML models that integrate catalog and operational features with various frameworks. Results indicate that feature-rich models outperform simpler ones in complex seismic environments like the Cooper Basin and laboratory cases, where seismicity is promoted by earthquake interaction and fault reactivation. However, in scenarios like St1 Helsinki, with minimal event clustering, additional features offer limited predictive benefits. While ML models are promising, several challenges impede reliable forecasting, including data scarcity from operational wells, the extrapolation demands of cumulative output (i.e. CSM), and the difficulty of predicting abrupt CSM increases for large seismic events. Enhancing model robustness requires synthetic data augmentation and improved feature selection capable of capturing diverse reservoir dynamics. These advancements may enable more accurate near real-time forecasts of problematic induced seismic events, informing operational decisions to mitigate seismic risks while maximizing energy extraction, and hence offering a pathway for broader adoption of ML in renewable energy development and management.