Papers

Three-dimensional probabilistic inversion of time-domain airborne electromagnetic (AEM) data is limited by the cost of the forward solve. Even though one simulation takes only tens of seconds, a Bayesian inversion of a survey of millions of soundings requires of order $10^{10}$ forward evaluations. To address this, we develop a continually learning neural-operator surrogate of the three-dimensional AEM forward operator that replaces the solver inside the Bayesian inversion. We start from the point of view that regardless of what geological prior is specified, Maxwell's laws remain invariant. Secondly, we avoid the limitation of learning on a single prior by continual learning on consecutive priors, which means our surrogate becomes richer as it is applied in future case studies, either by the authors, or by the scientific community. We use a validity check built on ensemble disagreement to divert cases with measurements outside the training range to the solver. Driven by the surrogate, the identical Markov chain Monte Carlo sampler reproduces the full-solver posterior, and its credible intervals cover the truth within 2.6 percentage points. Applied to the 2013 Capricorn TEMPEST survey in Western Australia, the surrogate inverts over two million soundings in seconds, a computation infeasible for the solver. Testing the geological prior against the entire survey costs minutes. The framework delivers uncertainty-quantified conductivity imaging at survey scale, which we believe is essential to perform near real-time mineral-systems targeting with geophysics.

We propose GeoFormer, a Geometry-Aware Transformer architecture specifically designed for prestack seismic data. Unlike Vision Transformer, whose tokens are extracted from 2D patches and primarily encode visual patterns, GeoFormer is designed for prestack seismic data by explicitly incorporating acquisition geometry. Each seismic trace is represented by a 5D unit consisting of the waveform and four source-receiver coordinates, from which two geometric attributes are derived: offset and relative elevation, where relative elevation is the receiver elevation minus the source elevation. GeoFormer therefore performs trace-level tokenization, where each token combines the waveform with these geometric attributes. To exploit these geometric attributes, GeoFormer introduces three geometry injection mechanisms operating at different levels of the Transformer pipeline. At the token level, GeomMLP replaces the classification token with a per-trace geometric representation derived from offset and elevation. At the normalization level, GeomAdaLN replaces uniform layer normalization with geometry-conditioned feature modulation. At the attention level, GeomAttnBias injects a parameter-free physical prior that geometrically proximate traces should attend more to each other. We validate GeoFormer on first-arrival picking, a representative seismic processing task that relies heavily on acquisition geometry. Experiments on four field datasets demonstrate that GeoFormer outperforms Vision Transformer and two task-specific baselines. Crucially, GeoFormer maintains robust picking accuracy under strong noise because its trace-level tokens encode geometric attributes that provide a physical prior independent of waveform quality. Ablation studies confirm that GeomMLP, GeomAdaLN, and GeomAttnBias each contribute to GeoFormer's picking accuracy.

Modeling frequency-dependent shear-wave propagation in fluid-saturated porous rocks remains a challenge today because existing poroelastic theories generally rely a large number of constitutive parameters that are difficult to constrain from seismic observations. We propose a dynamically consistent higher-gradient continuum model that reproduces shear-wave dispersion using only three effective material parameters: the macroscopic shear velocity and two characteristic length scales governing higher-order elasticity and gradient inertia. We validate the model against laboratory measurements of glycerin-saturated tight sandstone over a wide range of effective pressures. The inversion reproduces the measured shear-wave dispersion, including positive, weak and negative dispersion regimes. We finally discuss the relationship between the proposed formulation and micropolar theories, the computational challenges and the potential of reduced-order generalized continuum models for seismic exploration, reservoir characterization, and emerging agroseismology applications.

Waveform inversion seeks to estimate the wave speed of a heterogeneous, inaccessible medium, from time-resolved measurements of the waves at user controlled sensors. We consider this inverse problem for acoustic waves and an active array of source/receiver sensors that emit probing signals and measure the generated pressure waves. The forward map, from the wave speed to the measurements, is nonlinear and oscillatory. The oscillations cause cycle skipping, the main impediment to using the standard, nonlinear least-squares data fitting formulation, known as full waveform inversion (FWI). A recently introduced alternative waveform inversion approach computes from the measurements an algebraic surrogate of the wave operator, a reduced order model (ROM) matrix, which is then used to estimate the wave speed. The mapping from the measurements to the ROM is nonlinear, but well understood. It is computed efficiently, in a non-iterative manner. The nonlinear mapping from the ROM to the wave speed is less understood, and its approximation involves time-consuming optimization. Our goal in this paper is to use a neural network to map the ROM matrix to a nearby one, that has a simpler and explicit dependence on the wave speed. This simplifies and reduces the computational cost of the ROM-based waveform inversion. We introduce the methodology, called ROMNet, and test it with numerical simulations, using two training data sets: The first set consists of random media with variations of the wave speed modeled by a superposition of Gaussians with random amplitudes and standard deviations. The second is the publicly available GeoFWI dataset introduced for benchmarking FWI using deep learning. We compare the performance of ROMNet with the direct ROM-based inversion and with two representative deep learning approaches to FWI: “Fourier-DeepONet" and “InversionNet".

Nonlinear shear softening in soils and rocks is commonly described using empirical modulus-reduction laws. In this work, we investigate nonlinear softening within the Eringen–Mindlin micromorphic framework by interpreting its additional constitutive parameters in terms of coupling, internal relaxation, and micro-inertia. The micromorphic continuum is viewed as a dynamically enriched effective medium in which the macroscopic deformation interacts with an internal microstructural deformation field. This interaction naturally produces frequency-dependent effective elastic properties and provides a physical connection between macroscopic wave propagation and unresolved microstructural dynamics. We show that nonlinear shear softening can be described by a progressive reduction of the coupling between macroscopic and internal deformation. As this coupling weakens with increasing shear strain, the effective shear modulus decreases. In the low-frequency limit, the resulting constitutive relation takes a hyperbolic form closely related to the classical Hardin–Drnevich law. The model is validated against laboratory measurements of shear-modulus reduction in sand and reproduces the observed nonlinear behavior with only a small number of physically interpretable parameters. The proposed framework provides a new description of nonlinear softening and frequency-dependent behavior without requiring explicit knowledge of the underlying microstructural geometry. Classical empirical modulus-reduction curves emerge as a macroscopic consequence of the interaction between macroscopic and internal deformation.

Modeling frequency-dependent seismic wave propagation in porous media remains a major challenge because classical poroelastic theories generally require a large number of poorly constrained constitutive parameters. We investigate seismic wave propagation within a higher-order gradient constitutive framework as an effective dynamic homogenization model for poroelastic media. The proposed formulation extends classical strain-gradient elasticity by incorporating an additional curvature contribution associated with rotational deformation, allowing independent descriptions of longitudinal and transverse wave dispersion within a unified continuum framework. A general dispersion relation is derived together with its macroscopic and microscopic velocity limits, providing a direct physical interpretation of the characteristic inertial and elastic length scales introduced by the model. The framework is validated against laboratory measurements of frequency-dependent P- and S-wave velocities in tight sandstones under different confining pressures. The inversion demonstrates that a common inertial characteristic length and a common strain-gradient characteristic length govern both wave types, while an additional curvature characteristic length is required only for shear waves. Comparison with classical Biot theory shows that the proposed framework reproduces the observed laboratory dispersion using only a small number of effective constitutive parameters, providing a reduced constitutive description of poroelastic wave propagation. These results suggest that higher-order gradient elasticity offers a physically interpretable and computationally efficient alternative for modeling frequency-dependent seismic wave propagation, with potential applications to seismic imaging, inversion, and natural hydrogen reservoir characterization.

Recent JWST observations challenge the traditional 'cosmic shoreline' from both sides, revealing thick volatile atmospheres on the hottest close-in 'lava worlds,' where irradiation should drive the most extreme escape, and bare rocky surfaces on cooler terrestrial planets around M dwarfs, where atmospheres would be expected to survive. Using a coupled atmosphere-interior evolution model, we show that atmosphere retention is governed not by a single escape boundary but by two: a hot, outgassing-regulated 'cosmic sandbar' and a cooler, escape-regulated 'cosmic shoreline,' separated by an 'airless valley' that may mark a graveyard of stripped sub-Neptune cores. The sandbar arises because long-lived magma oceans, sustained further by tidal heating from secular eccentricity excitation in multi-planet systems, keep most volatiles dissolved and expose only a small atmospheric reservoir to escape, whereas cooler planets solidify, sequestering volatiles in the deep solid mantle while overexposing the rest to loss. This two-regime structure recasts the single cosmic shoreline as two boundaries set by distinct physics: outgassing and escape. We provide time-evolving fits for both boundaries across G, K, and M stellar types as a function of volatile inventory, planetary mass, age, and tidal heating. Lava worlds with thick atmospheres are unlikely around stars cooler than K-type unless sustained by extreme tidal and/or other interior heating. This framework links atmosphere survival from USP lava worlds to habitable zone planets, informing target selection and interpretation for TRAPPIST-1 and JWST DDT characterization.

We present a variational–adjoint interpretation of quantum mechanics in which the interaction between forward and adjoint wavefunctions defines a general sensitivity kernel of the system. This Fréchet interaction density quantifies how perturbations in the wavefunction (or system parameters) influence a chosen observable. The familiar Born probability density appears as a special case when the adjoint wavefunction is chosen as the complex conjugate of the forward wavefunction, for which the interaction density becomes real and non-negative. Within this framework, probability is a particular positive-definite form of sensitivity described by the Fréchet interaction density. In general, the variational–adjoint interpretation also produces Fréchet sensitivity kernels associated with other quantum observables, including momentum, energy, and spin. This suggests that the Born probability density belongs to a broader class of Fréchet sensitivity kernels associated with quantum observables. The proposed interpretation also provides a connection with time-symmetric interpretations of quantum mechanics and possible future applications in quantum control and quantum metrology.

Accurate traveltime computation for the eikonal equation is essential in seismic applications such as tomography and migration. The fast sweeping method (FSM) is widely used because of its unconditional stability and computational efficiency. We develop a 3D fast sweeping solver for vertical transverse isotropic (VTI) media that combines multiplicative factorization with a six-tetrahedron pyramidal stencil. The factorization removes the point-source singularity, while the stencil improves local accuracy. In the unfactored formulation, the six-tetrahedron finite-difference scheme requires solving only quadratic equations. After factorization, however, the local update for the perturbation factor becomes quartic, and the update systems on oblique stencil faces are substantially more complicated than those on non-oblique faces. To resolve these difficulties, we solve the update equation using Ferrari's method with a robust root-selection strategy and derive complete update formulas for all oblique-face configurations. For horizontally constrained faces, the characteristic constraint reduces to a linear relation by exploiting the structure of the VTI Hamiltonian, so the local system still reduces to a quartic equation. For mixed horizontal-vertical constrained faces, we design a bisection-based iterative solver. Numerical examples show that the proposed method effectively suppresses source-related errors and improves traveltime accuracy.

The Bandt-Pompe permutation entropy framework, alongside the complexity-entropy causality plane, has become a standard tool for characterizing the dynamical properties of time series. However, observational noise distorts ordinal pattern probability distributions in ways that can systematically misplace time series within the causality plane, compromising dynamical classification. This effect is particularly relevant for geophysical signals, which are typically poorly and irregularly sampled, and have a low signal-to-noise level. In this work, we characterize the distortions on ordinal pattern statistics using the formalism of majorization. We provide theoretical results and propose corrective strategies that restore discriminability under realistic measurement conditions. To achieve this, we introduce methodology that allows the characterization of noisy dynamical series and further allows the quantification of observational noise without the need of a fitting procedure. Finally, we illustrate our methodology by analyzing paleomagnetic records to determine if the geological evolution of the Earth dipole is better described by a stochastic or chaotic system.

Deep atmospheric convection governs South Asian monsoon variability, yet attempting to learn its latent world model directly from zero-inflated, heavy-tailed precipitation yields suboptimal predictive representations. Continuous atmospheric proxies, such as outgoing longwave radiation (OLR), express this convective organization far more coherently. We address this mismatch with cross-variable proxy learning: M-JEPA, a multiscale Monsoon Joint-Embedding Predictive Architecture, is pretrained on five continuous proxy fields over Lagrangian patches tracking moving convective systems—without rainfall supervision at any point. The resulting frozen representation is transferred to daily precipitation forecasts through a shared decoder trunk featuring parallel probabilistic and deterministic branches. Because rainfall is strictly unobserved during pretraining, downstream skill directly measures the predictive information captured in the latent rollout. A frozen-backbone probing framework with two controls (an identical architecture trained on rainfall alone, and a randomly initialized backbone) attributes the transfer specifically to proxy pretraining: direct rainfall training exhibits $36\%$ higher CRPS error ($7.52$ vs.\ $5.54$\,mm/day). Against the 51-member operational ECMWF ensemble, the transferred model attains a statistically resolved CRPS advantage ($6.81$ vs.\ $6.89$\,mm/day) and higher Brier skill ($+0.05$ vs.\ $-0.04$) using $15.4$M parameters on a single consumer GPU, concentrated at heavy-rain thresholds and fine spatial scales, while the ensemble retains an advantage in neighborhood skill and deterministic references on point metrics. The result provides a competitive monsoon precipitation forecast grounded in intraseasonal dynamics and a diagnostic framework for evaluating transferred atmospheric representations.

This paper presents a scalable Bayesian framework for the inference of modern-day sea-level change and surface mass redistribution. To motivate this approach, we first review and quantify the performance of some standard methods for the analysis of satellite gravity and ocean altimetry observations. We find that these methods substantially underestimate uncertainties and are subject to systematic biases, deficiencies that stem from their incomplete treatment of sea-level physics and from uncertainty estimates based solely on the propagation of observational noise. To address these limitations, our approach combines three advances. First, we use recent developments in adjoint sea-level theory to embed the full physics into the forward and inverse modelling. Second, we formulate the Bayesian inverse problem in an infinite-dimensional setting, thereby avoiding discretisation artefacts and the underestimation of uncertainties inherent in truncated model spaces. Finally, our computational methods render such inversions tractable at full observational resolution while supporting joint model spaces and multiple data types. By employing a matrix-free approach with iterative solvers and randomised low-rank decompositions – implemented in the linked open-source libraries pygeoinf and pyslfp – basic calculations are possible on a single laptop, with the most intensive tasks parallelising trivially across available cores. We demonstrate the methodology through a series of synthetic experiments, culminating in a joint inversion of satellite gravity and ocean altimetry data that decomposes regionally averaged sea-level change into steric and manometric components with quantified uncertainties, including the degeneracies that remain.

Urban wildfire resilience depends on interactions among built form, vegetation condition, and extreme fire weather, yet city-scale risk models often overlook whether predictive skill transfers across neighborhoods. We developed a spatially validated GeoAI workflow for the January 2025 Palisades Fire, linking 12,081 CAL FIRE damage inspections to pre-fire Sentinel-2 vegetation indices, Landsat surface temperature, LANDFIRE fuels, terrain, and OpenStreetMap buildings and roads. Among 9,883 inspected residential structures, 5,566 were destroyed. Random cross-validation yielded ROC-AUC 0.92 for the integrated XGBoost model, but 1 km spatial block validation reduced performance to 0.75; logistic regression performed similarly and was better calibrated. Building count within 100 m was the strongest predictor, with destruction odds increasing 4.12-fold per standard deviation. Vegetation moisture and greenness showed opposing conditional associations: NDMI at 100-300 m was protective (OR 0.52), whereas NDVI at 30-100 m was positively associated with destruction after accounting for moisture (OR 1.74). Predictive information was concentrated at the 100-300 m neighborhood scale. A separate post-fire track mapped burn severity and vegetation recovery without leakage. The results support neighborhood-scale susceptibility screening, moisture-aware vegetation management, and spatial block validation as a minimum standard for single-event urban wildfire modeling.

Quantifying the interactions between neighboring projects will be essential as carbon storage operations expand in scale. However, many existing basin-scale studies are deterministic or rely on hypothetical scenarios. In addition, simulations that do not fully account for neighboring projects may suffice for initial permitting. In this study, we develop a unified basin-scale model for the Broom Creek Formation in North Dakota. The model includes a total of five (existing and planned) projects, involving 18~injection wells with a maximum potential injection rate of 30~MTPA. The unified model, constructed by integrating publicly available geological and project data, contains $44\times10^6$ cells and requires over two days to run. To enable faster simulations, we present a multilevel treatment that combines four-level nested local grid refinement, constructed through the solution of a constraint-satisfaction problem, with optimized power averaging for permeability in coarse regions. After demonstrating the accuracy of the coarse model, we use it to assess uncertainty in key quantities of interest (QoIs), including total CO$_2$ injected and plume area after 50~years, across a wide range of geological realizations and model parameters. In addition, we develop a deep learning surrogate model, which is used in global sensitivity analyses to identify dominant contributors to QoI variance. Comparisons between the unified model and standalone project simulations are presented to quantify inter-project interference. This can be substantial in some cases, e.g., for one project, the median total CO$_2$ injected decreases from 271~MT (standalone) to 192~MT (unified model).

We present a multi-parameter simulation-based inference framework for joint Bayesian recovery of subsurface velocity and acoustic impedance from seismic data. A score-based diffusion model is conditioned on two complementary Common Image Gathers (CIGs): an inverse-scattering CIG encoding reflectivity amplitude and an anti-ISIC CIG encoding kinematic velocity errors. The model simultaneously samples the posterior distributions of both parameters. Training labels are deliberately decoupled to prevent the model from exploiting the Gardner relationship: velocity targets are lightly smoothed to match the long-wavelength content of the anti-ISIC CIG, while impedance targets retain the unsmoothed ground truth. On the Compass benchmark the model achieves velocity SSIM of 0.967 (RMSE 0.050 km/s) and impedance SSIM 0.867 (RMSE 0.279 km/s g/cm^3), with velocity quality confirmed by CIG focusing.

Simulating the trajectory of surface drifters in the ocean matters for search and rescue, pollution tracking, oil and chemical spill response, and marine risk analysis. Accurate prediction remains difficult, as widely acknowledged in the literature, and also illustrated by the “Forecasting Floats in Turbulence” challenge issued by the US Defense Advanced Research Projects Agency (DARPA) in 2021, and which ultimately led to this paper. The main source of error usually comes from uncertain ocean currents, while errors in wind forcing and object drift properties are often smaller [Dagestad and Röhrs, 2019]. Here, we use an open one-year dataset of Sofar Spotter trajectories together with several ocean and atmospheric hindcast products to test data-driven drift models at scale. We compare three approaches: i) a standard (baseline) drifter trajectory simulation based on one ocean model and one atmospheric model, ii) linear regression (LR) models that fuse all available predictors, and iii) neural networks (NN) using similar inputs. A simple LR model that combines all predictors performs equally well as the NN. Because LR is simpler, cheaper, and more robust, we retain it as the preferred approach. In 2-day trajectory prediction, this improves the Liu-Weisberg skill score by around 40% relative to the baseline. These findings apply to hindcast mode; applying this methodology for forecast mode remains for future work.

Data-driven methods for earthquake and tsunami early warning rely on large ensembles of rupture scenarios and their resulting waveforms, but generating such datasets with repeated high-fidelity seismic and tsunami simulations is prohibitively expensive. We exploit the linear time-invariant structure of both dynamics to precompute elastic Green's functions and acoustic-gravity adjoint responses, reducing the source-to-waveform map to two consecutive convolution operators. We evaluate these convolutions with a distributed, FFT-accelerated GPU pipeline that partitions the large seafloor grid across GPUs and directly generates the final observation waveforms. We demonstrate the scalability of this pipeline for the Cascadia Subduction Zone with 963 subfaults, 2,416,530 seafloor grid points, 64 observation locations, and 256 timesteps, requiring 9.45 TiB of aggregate GPU memory. On 64 GB200 GPUs within one NVL72 domain, the pipeline generates waveforms in 24 ms per rupture once the response operators are resident, enabling large rupture ensembles to be evaluated within minutes.

The French Massif Central and the Eifel Volcanic Field represent two of the enigmatic features in France and Germany, respectively. Both areas have been affected by Cenozoic volcanism whose origin is still debated. Using 3-D numerical geodynamic modelling of recently developed solid-mechanical constitutive laws including viscoelastic-viscoplastic behaviour, we explored a new hypothesis to highlight the contribution of localised deformational (inelastic) heating towards the evolution of localised volcanism as being due to far-field-induced compressive stresses from the Africa-Eurasia convergence. Our simulations show that variations in crustal thicknesses can help localise inelastic strain near the Moho of relatively thicker crust, which leads to enhanced heating. Localised viscoplastic deformation above the brittle-ductile transition in the crust also formed conjugate deformation bands, resulting in localised changes in the topography which can be compared to present-day topography. The first-order picture that emerges from our simulations shows an association of elevated topography and tectonic structures that localise deformation, with a geothermal anomaly sufficient to explain the elevated heat flow and cause a potential localised zone for partial melting in hydrated or CO\textsubscript{2}-bearing mantle below the Moho.

Flood-monitoring sensors are often placed according to coverage, access, or expected inundation. However, the value of a measurement depends on the prediction or decision it is intended to inform. Using tRIBS-Urban simulations and a neural-network surrogate of the August 2014 metropolitan Detroit flood, we examine how this learning target changes single-sensor placement. Across 2,576 candidate locations, we compare parameter-oriented optimal experimental design (PO-OED), which values expected information gain (EIG) about model parameters, with goal-oriented optimal experimental design (GO-OED), which values EIG about specified flood predictions. We also examine how parameter EIG evolves during the event, and illustrate that parameter learning translates unevenly into reductions in predictive uncertainty across locations and lead times. Under GO-OED, point-depth targets favor nearby locations, whereas regional-average and regional maximum-depth targets can favor nonlocal locations. Weighted multi-point objectives retain similar broad spatial patterns, although their computed max-EIG locations differ. Public geospatial data further provide illustrative feasibility and contextual classifications for deployment screening. These results show how monitoring objectives shape sensor placement in optimal experimental design, and motivates an objective-first workflow that defines the intended prediction and priorities, applies field-verified restrictions, and ranks locations by EIG.

Video recordings can reveal how rapidly fault displacement develops at the Earth's surface, but camera motion and recording artifacts can obscure the ground signal. We analyze a secondary copy of security-camera footage that captured surface displacement during the 28 July 2026 Kumamoto earthquake; the native recording was unavailable and could not be recovered. Two independent image-tracking methods were used. Optical flow follows identifiable image features, whereas normalized cross-correlation (NCC) template matching follows fixed image patches by their similarity. Both measured target-region motion relative to spatially separated reference regions while correcting motion shared by the recording. Image displacement was calibrated to the magnitude of the field-measured offset vector: 1.05 m right-lateral and 0.90 m east-side-up, or 1.383 m in total. We characterize the principal rise by the time required for displacement to progress from 20% to 80% of the selected final level. Across prespecified endpoint choices, optical flow gives 0.866-0.901 s and NCC gives 0.910-0.928 s. These durations correspond to average rates of 0.920-0.958 and 0.895-0.912 m/s, respectively. Checks using independently published tracking windows reproduce the displacement scale, although exact timing is more sensitive in spatially restricted tests. The record also shows an early apparent peak and decline followed by renewed apparent horizontal displacement. Because that later motion may represent either continued ground displacement or the geometry of the secondary recording, neither the permanent endpoint nor physical overshoot can be determined. The most robust conclusion is that the central part of the surface displacement developed in approximately 0.9 s at an average rate near 0.9 m/s.

Periodic surface-wave forcing can reorganise turbulent flows through coherent spectral response, synchronization, intermittency, and changes in short-term predictability, so its dynamical effect need not vary monotonically with forcing frequency. We reanalyse laboratory acoustic Doppler velocimetry records obtained at constant discharge under four conditions (0, 0.5, 0.67, and 1~Hz) using one common KNN–GMAE forecast-error-growth protocol. A distance-weighted $k$-nearest-neighbour predictor generates out-of-sample forecasts over multiple horizons, and the slope of the early quasi-linear region of $\ln(\mathrm{GMAE})$ versus physical forecast time is reported as a finite-horizon forecast-error-growth rate, $\lFEG$. For the full 120-s records, $\lFEG$ is 5.68, 0.85, 3.39, and 4.11~s$^{-1}$ for 0, 0.5, 0.67, and 1~Hz, respectively. The ordering 0~Hz $>$ 1~Hz $>$ 0.67~Hz $>$ 0.5~Hz is preserved in all seven nearby parameter configurations, indicating that the comparative result is not an artefact of a single KNN setting. A 40-s sliding-window analysis with a 10-s step reveals substantial temporal structure. The no-wave condition remains predominantly high and the 0.5-Hz condition predominantly low, whereas the 1-Hz record has weaker local support for a single exponential-growth regime: only 3 of 9 windows satisfy the adopted early-fit criterion $\RFEG\geq0.90$, compared with 8/9, 7/9, and 7/9 for 0, 0.5, and 0.67~Hz.

We used 3-D thermomechanical modelling to investigate conditions during subduction-zone initiation and early thermal development with focus on the Lesser Antilles. Our model imposes a convergence velocity of 2 cm per year and incorporates heating caused by irreversible deformation of mantle and crustal rocks, using elasticity, creep, and non-associative plastic flow laws. Our results show that deformational heating before slab development is unexpectedly strong. After several million years, buckling and heating due to irreversible deformation create distinctive patterns of topography and surface heat flow that resemble present-day observations, despite the slab and subduction interface being incompletely developed. Within the Caribbean plate, plate buckling produces a high topographic ridge underlain by a large positive thermal anomaly of approximately 200 K, centred just below the Moho. The conductive thermal boundary layer transporting this heat to the surface thins from about 100 km to 10 km beneath the topographic maximum, allowing the ridge to rise above sea level. This thermal structure suggests the potential initiation of a volcanic arc approximately 180 km from the inter-plate contact. A hot zone at 30-50 km depth has pressures consistent with those inferred from Lesser Antilles primitive magmas and represents the most plausible location for partial melting of Caribbean mantle if volatiles are present. The thick Caribbean crust, approximately 20-25 km, is also heated sufficiently for possible silicic melt generation. The inferred lithospheric thickness of 50-100 km aligns with tomography studies. Thus, subduction thermal structure is strongly influenced by several million years of initiation processes.

Timely, paddock-scale estimates of plant-available water (PAW) can support crop management in water-limited grain systems. We present the Sydney Soil Water-Energy Balance (SWEB) model, a scalable, physically consistent modelling framework that integrates satellite and climate data with soil properties to estimate crop evapotranspiration and root-zone soil moisture (RZSM) at daily, 30 m resolution. SWEB was validated nationally against different SM monitoring networks across Australia, which demonstrated robust performance across diverse grain-growing environments (correlation coefficients generally ranged from 0.70 to 0.85 for most regions). In future applications, mid-season PAW nowcasts could be combined with water-use-efficiency approaches to estimate water-limited yield potential and inform responsive management decisions, such as nitrogen fertiliser top-up recommendations.

Surface displacements, including co-seismic fault slip, have occasionally been captured by closed-circuit television (CCTV) cameras. Such video recordings may provide information unavailable from seismometers located away from the fault. We analyzed videos of surface displacement from three earthquakes—the 2018 Hualien earthquake, the 2025 Mandalay earthquake, and the 2026 Kumamoto earthquake—using a unified methodology. We removed the effects of strong ground motion by aligning each frame so that reference regions in the image remained fixed. This analysis enables the temporal evolution and, in particular, the duration of surface displacement to be estimated. Although the three earthquakes ranged from $M_\mathrm{w}$6.4 to $M_\mathrm{w}$7.7, the tracked motion lasted approximately 2 s in all three cases; in the Hualien case, the motion included compression of the ground rather than fault slip alone. Comparison with field surveys indicates that the maximum velocity was on the order of 1 m/s in each case. The 2026 Kumamoto earthquake, however, exhibited a more complex velocity history than the other two events. The three earthquakes involved different mechanisms of surface fault displacement, suggesting that these differences may also influence the complexity of the displacement time histories.

Even for tension-dominated fracture propagation, porous materials may deform plastically adjacent to the propagating fracture. As is common for porous materials, existing phase-field models typically employ a non-associative flow rule for plasticity, and a Helmholtz free energy based on strain and fluid pressure. This work revisits the hydromechanically coupled formulation of the phase-field model for fracture in poro-elastoplastic media by analyzing the strength surface and fracture driving force. Our analyses show that these common choices of flow rule and free energy will lead to a discontinuous strength surface across the tension-compression transition. A non-associative flow rule introduces a jump at the strength surface, while treating fluid pressure-rather than fluid content-as the independent variable in the Helmholtz free energy omits a coupling term from the phase-field driving force, also breaking continuity. Incorporating an associative Drucker-Prager flow rule and this omitted coupling term ensures a continuous strength surface and the accurate fracture driving force. The proposed model exhibits improved accuracy in hydromechanical responses when compared against the analytical solution of the Kristianovich-Geertsma-de Klerk hydraulic fracturing benchmark. Numerical simulations of hydraulic fracturing and biaxial compression in poro-elastoplastic media show that the model can reproduce both shear-dominated fractures induced by mechanical disturbance and tension-dominated fractures driven by fluid injection in saturated porous media.

Horizontal earthquake ground motions exhibit substantial changes in amplitude with changes in orientation. While the effects of style of faulting on central tendency measures of ground motion intensity such as RotD50 are well established in ground motion models, less attention has been given to how faulting mechanism influences the complete directional characteristics of ground motions. This study addresses this gap by conducting a comprehensive probabilistic characterization of ground motion directionality for strike-slip and reverse-slip earthquakes using 3,862 records with moment magnitudes greater than five. Two directionality parameters were examined across 58 oscillator periods ranging from 0.01 to 10 s. The directionality parameters examined consist of pseudo-acceleration response spectral ordinates at all non-redundant orientations normalized by either RotD100 or RotD50 intensities (denoted as {\eta} and {\nu}, respectively). Nonlinear regression was used to develop mechanism-specific models for the geometric mean, logarithmic standard deviations, and full probability distributions of these parameters. Results show that strike-slip earthquakes exhibit notably different directionality compared to reverse-slip earthquakes, particularly at longer periods. Strike-slip earthquakes are, on average, 1% to 44% more polarized than reverse-slip earthquakes at periods above 0.5 s, and have RotD100/RotD50 ratios that are 1% to 3.7% larger. The developed models enable engineers and engineering seismologists to obtain mechanism-specific estimates of pseudo-acceleration response spectral ordinates at specific orientations, for applications such as regional seismic risk assessments and orientation-specific probabilistic seismic hazard analyses.

Earthquake sequence forecasting requires models that can learn nonlinear history dependence while retaining robust statistical structure. We develop a scaling-law-informed neural marked point process, termed Fusion, that combines neural representations of catalog history with temporal features derived from the Epidemic-Type Aftershock Sequence model and magnitude information derived from the Gutenberg–Richter law. The model separates the magnitude cutoff applied to the input catalog from the fixed target-event threshold, allowing lower-magnitude earthquakes to inform forecasts without changing the target-event set. For the 2016–2017 Amatrice–Visso–Norcia sequence, Fusion achieves the highest target-event temporal likelihood when lower-magnitude events are retained, outperforming both ETAS and a purely neural point-process baseline. Event-wise and cumulative analyses show sustained timing gains through substantial portions of the Visso and Norcia sequences. Across five benchmark catalogs, catalog-specific neural training with a fixed ETAS prior yields the highest temporal likelihood at the minimum evaluated magnitude cutoff. Magnitude likelihood shows no consistent predictive gain beyond the Gutenberg–Richter-based ETAS reference, indicating that the additional information captured by Fusion is primarily temporal. These results show that lower-magnitude catalog histories and empirical scaling-law information complement neural sequence learning for target-event timing.

The prediction of equilibrium beach profiles under tidal influence is of fundamental importance for sustainable coastal development, informing shoreline protection strategies and managing coastal ecosystems under changing environmental conditions. However, it remains challenging due to the highly nonlinear interactions among wave, tide, and sedimentary processes. Traditional empirical and numerical models often exhibit limited adaptability across diverse coastal environments, with especially pronounced limitations in beach systems where tidal processes are important . To improve data-driven prediction under these conditions, this study proposes MorphoGP, a unified category-specific Gaussian process framework for predicting equilibrium beach profiles (EBPs) under tidal influence. The framework first introduces a ContourCluster model based on contrastive learning to classify tide-influenced beach morphologies automatically. Within each morphological category, a specialized Gaussian process expert learns statistical associations between environmental descriptors including waves, tides, and sediments and the beach profile's shape. A Gating Net then integrates the outputs of all experts through a probabilistic weighting mechanism to produce the final prediction. Evaluated on data from over 180 beach profiles from tide-influenced coasts along the Chinese coast, MorphoGP achieves improved predictive performance compared with conventional and deep learning models, reducing the test RMSE by about 59.3% compared with the best baseline and achieving a final RMSE of 0.297 m. The proposed framework provides a physically informed, data-driven tool for equilibrium beach-profile prediction under tidal influence and coastal management, while stronger process-level physical coupling remains an important direction for future development.

Classical seismic data reconstruction relies on manually designed structural priors and iterative operators, whose coupled design space is far larger than manual trial and error can explore systematically. Deep-learning methods encode the reconstruction rules in learned weights rather than in an explicit operator that can be inspected and modified. We propose SeisEvo (Seismic Algorithm Evolution), which does not optimize a single reconstruction result but searches for the algorithm that produces it. Starting from a classical reconstruction algorithm, an LLM-driven multi-agent search modifies only the components that the user has opened for editing, without prescribing the mechanism to be discovered. Candidates that violate the physical constraints of the task are rejected outright, and the remaining ones are scored by execution. The output is neither an agent system nor a neural network, but a standalone white-box algorithm that requires no agent or neural network at inference time. For interpolation without added noise, the search discovered a residual-gated, phase-aligned dip-consistency projection; Evo-POCS improves the SNR over classic POCS by 3.49 dB on average across missing ratios from 30% to 70%. For simultaneous interpolation and denoising, it discovered a reliability-grouped singular-value shrinkage; Evo-MSSA improves the average reconstruction SNR by more than 7 dB over classic MSSA and by more than 3 dB over a stronger rank-reduction baseline. Both operators retain their gains on data not used during the search. To the best of our knowledge, this is the first study to formulate the design of a seismic reconstruction operator as a constrained, LLM-driven program evolution task. Agentic algorithm evolution can thus complement deep learning in discovering explicit, inspectable, and deployable seismic processing algorithms.

The complex-step method is traditionally derived from the Taylor expansion of an analytic function and is widely used as a numerical technique for derivative approximation. We present an alternative formulation based on the Cauchy–Riemann equations and show that the classical complex-step relation arises naturally from the harmonic structure of holomorphic functions. In particular, the complex-step method admits two complementary harmonic interpretations: as a Cauchy problem, in which the derivative is identified with the normal datum of the imaginary component on the real axis, and as a reconstruction problem in a strip, in which the finite imaginary perturbation provides the upper-boundary data. The latter formulation leads explicitly to the strip Poisson and conjugate Poisson kernels and their derivatives. A related harmonic reconstruction framework in the upper half-plane leads to the Poisson, conjugate Poisson, and Cauchy kernels as elementary reconstruction operators for harmonic and holomorphic functions. Extending this reconstruction from ordinary boundary functions to finite measures yields the classical Stieltjes transform and its inversion formula. The same measure-theoretic structure appears in spectral theory, where scalar matrix elements of the resolvent are Stieltjes transforms of the associated spectral measures. These results establish a common complex-analytic structure connecting complex-step differentiation, harmonic reconstruction, Stieltjes inversion, and spectral reconstruction, while distinguishing the boundary-value problems through which the corresponding information is recovered.