Model simulations of past climates are increasingly found to compare well with proxy data at a global scale, but regional discrepancies remain. A persistent issue in modeling past greenhouse climates has been the temperature difference between equatorial and (sub-)polar regions, which is typically much larger in simulations than proxy data suggest. Particularly in the Eocene, multiple temperature proxies suggest extreme warmth in the southwest Pacific Ocean, where model simulations consistently suggest temperate conditions. Here we present new global ocean model simulations at 0.1° horizontal resolution for the middle-late Eocene. The eddies in the high-resolution model affect poleward heat transport and local time-mean flow in critical regions compared to the non-eddying flow in the standard low-resolution simulations. As a result, the high-resolution simulations produce higher surface temperatures near Antarctica and lower surface temperatures near the equator compared to the low-resolution simulations, leading to better correspondence with proxy reconstructions. Crucially, the high-resolution simulations are also much more consistent with biogeographic patterns in endemic-Antarctic and low-latitude-derived plankton, and thus resolve the long-standing discrepancy of warm subpolar ocean temperatures and isolating polar gyre circulation. The results imply that strongly eddying model simulations are required to reconcile discrepancies between regional proxy data and models, and demonstrate the importance of accurate regional paleobathymetry for proxy-model comparisons.

We conduct a high-resolution teleseismic receiver function investigation of the subducting plate interface within the Alaskan forearc beneath Kodiak Island using data collected as part of the Alaska Amphibious Community Seismic Experiment in 2019. The Kodiak node array consisted of 398 nodal geophones deployed at ~200-300 m spacing on northeastern Kodiak Island within the southern asperity of the 1964 Mw9.2 Great Alaska earthquake. Receiver function images at frequencies of 1.2 and 2.4 Hz show a coherent, slightly dipping velocity increase at ~30-40 km depth consistent with the expected slab Moho. In contrast to studies within the northern asperity of the 1964 rupture, we find no evidence for a prominent low-velocity layer above the slab Moho thick enough to be resolved by upgoing P-to-S conversions. These results support evidence from seismicity and geodetic strain suggesting that the 1964 rupture connected northern (Kenai) and southern (Kodiak) asperities with different plate interface properties.

Although the Brune source model describes earthquake moment release as a single pulse, it is widely used in studies of complex earthquakes with multiple episodes of high moment release (i.e., multiple subevents). In this study, we investigate how corner frequency estimates of earthquakes with multiple subevents are biased if they are based on the Brune source model. By assuming complex sources as a sum of multiple Brune sources, we analyze 1,640 source time functions (STFs) of Mw 5.5-8.0 earthquakes in the SCARDEC catalog to estimate the corner frequencies, onset times, and seismic moments of subevents. We identify more subevents for strike-slip earthquakes than dip-slip earthquakes, and the number of resolvable subevents increases with magnitude. We find that earthquake corner frequency correlates best with the corner frequency of the subevent with the highest moment release (i.e., the largest subsevent). This suggests that, when the Brune model is used, the estimated corner frequency and therefore the stress drop of a complex earthquake is determined primarily by the largest subevent rather than the total rupture area. Our results imply that the stress variation of asperities, rather than the average stress change of the whole fault, contributes to the large variance of stress drop estimates.

The fastest projected rates of sea level rise appear in models which include “the marine ice cliff instability (MICI),” a hypothesized but mostly unobserved process defined by rapid, brittle failure of terminal ice cliffs that outpaces viscous relaxation and ice-shelf formation. The response of Crane Glacier to the Larsen B Ice Shelf collapse has been invoked as observational evidence of this process, but limited data coverage in space and time has limited our ability to meaningfully refine our understanding of cliff failure processes using that event. Updates to Crane’s subglacial topography show that much of its terminus retreat occurred in floating, not grounded ice, but its retreat, arrest, and regrowth over the last decade indicate brittle (not viscous) processes dominated during the 2 years following ice shelf collapse. If retreat occurred by cliff failure, maximum cliff heights would have been 111 m, consistent with process models that incorporate damaged ice.

Characterizing the large M4.7+ seismic events during the 2018 Kilauea eruption is important to understanding the complex subsurface deformation at the Kilauea summit. The first 12 events (May 17 - May 26) are associated with the explosive eruptions and the remaining 50 events (May 29 - August 02) onwards are accompanied by large-scale caldera collapses. It is challenging to resolve their locations and mechanisms because of the shallow source depths, complex velocity structure, and significant non double-couple components. We show the necessity to combine multiple geophysical data including broadband seismometers, accelerometers and infrasound sensors to resolve different aspects of the seismic source. The seismic moment tensor solutions using near-field summit stations show the early events are highly isotropic. Infrasound data and particle motion analysis identify the Halema’uma’u reservoir as the inflation source. For the later collapse events, two independent moment tensor inversions using local and global stations consistently show that asymmetric slips occur on inward-dipping normal faults along the northwest corner of the caldera. Infrasound simulation results suggest there may be inflation during the collapse events although not resolvable seismically. Our findings show that the summit events are characterized by both inflation and asymmetric slip, which are consistent with geodetic data. Based on the location of the slip and microseismicity, the caldera may have failed in a ‘see-saw’ manner: small continuous slips in the form of microseismicity on the southeast corner of the caldera, compensated by large slips on the northwest during the large collapse events.

Geodetic altimeters provide unique observations of the river surface longitudinal profile due to their long repeat periods and densely spaced ground tracks. This information is valuable for calibrating hydraulic model parameters, and thus for producing reliable simulations of water level for flood forecasting and river management, particularly in poorly instrumented catchments. In this study, we present an efficient calibration approach for hydraulic models based on a steady-state hydraulic solver and CryoSat-2 observations. In order to ensure that only coherent forcing/observation pairs are considered in the calibration, we first propose an outlier filtering approach for CryoSat-2 observations in data-scarce regions using simulated runoff produced by a hydrologic model. In the hydraulic calibration, a steady-state solver computes the WSE profile along the river for selected discharges corresponding to the days of CryoSat-2 overpass. In synthetic calibration experiments, the global search algorithm generally recovers the true parameter values in portions of the river where observations are available, illustrating the benefit of dense spatial sampling from geodetic altimetry. The most sensitive parameters are the bed elevations. In calibration experiments with real CryoSat-2 data, validation performance against both Sentinel-3 WSE and in-situ records is similar to previous studies, with RMSD ranging from 0.43 to 1.14 m against Sentinel-3 and 0.60 to 0.73 against in-situ WSE observations. Performance remains similar when transferring parameters to a one-dimensional hydrodynamic model. Because the approach is computationally efficient, model parameters can be inverted at high spatial resolution to fully exploit the information contained in geodetic CryoSat-2 altimetry.

On 15 January 2022, the submarine volcano on the southwest Pacific island of Tonga violently erupted. Thus far, the ionospheric oscillation features caused by the volcanic eruption have not been identified. Here, observations from the Super Dual Auroral Radar Network (SuperDARN) radars and digisondes \change{are}{were} employed to analyze ionospheric oscillations in the Northern Hemisphere caused by the volcanic eruption in Tonga. Due to the magnetic field conjugate effect, the ionospheric oscillations were observed much earlier than the arrival of surface air pressure waves, and the maximum negative line-of-sight (LOS) velocity of the ionospheric oscillations exceeded 100 m/s in the F layer. After the surface air pressure waves arrived, the maximum LOS velocity in the E layer approached 150 m/s. A maximum upward displacement of 100 km was observed in the ionosphere. This work provides a new perspective for understanding the strong ionospheric oscillation caused by geological hazards observed on Earth.

With the utilization of a novel synergistic approach, we constrain the vertical distribution of water vapor on Mars with measurements from nadir-pointing instruments. Water vapor column abundances were retrieved simultaneously with PFS (sensing the thermal infrared range) and SPICAM (sensing the near-infrared range) on Mars Express, yielding distinct yet complementary sensitivity to different parts of the atmospheric column. We show that by exploiting a spectral synergy retrieval approach, we obtain more accurate water vapor column abundances compared to when only one instrument is used, providing a new and highly robust reference climatology from Mars Express. We present a composite global dataset covering all seasons and latitudes, assembled from co-located observations sampled from seven Martian years. The synergy also offers a way to study the vertical partitioning of water, which has remained out of the scope of nadir observations made by single instruments covering a single spectral interval. Special attention is given to the north polar region, with extra focus on the sublimation of the seasonal polar cap during the late spring and summer seasons. Column abundances from the Mars Climate Database were found to be significantly higher than synergistically retrieved values, especially in the summer Northern Hemisphere. Deviances between synergy and model in both magnitude and meridional variation of the vertical confinement were also discovered, suggesting that certain aspects of the transport and dynamics of water vapor are not fully captured by current models.

As a consequence of a diminishing sea ice cover in the Arctic, activity is on the rise. The position of the sea ice edge, which is generally taken to define the extent of the ice cover, changes in response to dynamic and thermodynamic processes. Forecasts for sea ice expansion due to an advancing ice edge will provide information that can be of significance for operations in polar regions. However, the value of this information depends on the quality of the forecasts. Here, we present methods for examining the quality of forecasted sea ice expansion and the geographic location where the largest expansion are expected from the forecast results. The algorithm is simple to implement, and an examination of two years of model results and accompanying observations demonstrates the usefulness of the analysis.

El Niño-Southern Oscillation (ENSO) shows a large diversity of events that is modulated by climate variability and change. The representation of this diversity in climate models limits our ability to predict their impact on ecosystems and human livelihood. Here, we use multiple observational datasets to provide a probabilistic description of historical variations in event location and intensity, and to benchmark models, before examining future system trajectories. We find robust decadal variations in event intensities and locations in century-long observational datasets, which are associated with perturbations in equatorial wind-stress and thermocline depth, as well as extra-tropical anomalies in the North and South Pacific. Some climate models are capable of simulating such decadal variability in ENSO diversity, and the associated large-scale patterns. Projections of ENSO diversity in future climate change scenarios strongly depend on the magnitude of decadal variations, and the ability of climate models to reproduce them realistically over the 21st century.

The diurnal-eastward propagating tide with zonal wavenumber 3 (DE3) has gained significant attention due to its ability to preferentially propagate to the ionosphere and thermosphere (IT) from the tropical troposphere, thus effectively coupling these atmospheric regions. In this work, we demonstrate the existence of a pronounced zonal wavenumber 4 (WN4) structure in the low-latitude ionosphere during May 27 - June 5, 2020 using concurrent in-situ total ion number density measurements from the Scintillation Observations and Response of The Ionosphere to Electrodynamics (SORTIE) and the Ionospheric Connection Explorer (ICON) satellites. Temperature observations from the Thermosphere Ionosphere Mesosphere Energetics Dynamics Sounding of the Atmosphere using Broadband Emission Radiometry (TIMED/SABER) instrument near 105 km and output from the Specified Dynamics Whole Atmosphere Community Climate Model with thermosphere and ionosphere eXtension (SD/WACCM-X) demonstrate that this global-scale ionospheric WN4 structure is due to DE3 propagating through the lower thermosphere.

Tropical instability waves (TIWs) in the equatorial eastern Pacific (EEP) exhibit 25–40-day westward-propagating fluctuations with seasonal and inter-annual variations, which are stronger during the July–December and La Niña periods. They likely transfer their energy northward by forming barotropic Rossby waves (BTRWs). Long-term near-bottom current measurements at 10.5°N and 131.3°W during 2004–2013 revealed a spectral peak at 25–40 days, where significant coherences were found with satellite-measured sea surface height in a wide region of EEP with maxima around 5°N. Simulated deep currents from a data-assimilated ocean model agree with the observed near-bottom currents, and both currents vary seasonally and interannually, consistent with the typical characteristics of TIW. Further analyses using 25–40-day bandpass-filtered barotropic velocity data from the model revealed that they reasonably satisfied the theoretical dispersion relation of TIW-induced BTRW (BTRWTIW). This study provides the first observational evidence showing BTRWTIW propagating northward above 10°N in the northeastern Pacific.

The 39.8 ka Campanian Ignimbrite (CI) is the largest caldera-forming eruption of the Campi Flegrei during the Quaternary, which had a global-scale impact on the environment and human populations. The cooling following the eruption and the several effects of it strongly affected the paleoenvironment and the migration of hominids in Europe. The volume of the eruption is necessary to constrain the climate model of this area in the past. However, despite a large number of studies, the Dense Rock Equivalent (DRE) volume estimates range from 60 to 300 km. Here we present a review of the previous volume evaluations and a new calculation of the volume of the ignimbrite. This estimate is constrained by an isopach map that reconstructs the paleo-topography during the eruption. The preserved total bulk extra-caldera volume of the ignimbrite is estimated at 61.5 km ± 5.5 km. The total PDC deposit volume is then corrected for erosion, ash elutriation, the intracaldera deposit volume and the volume of tephra deposited in the sea. The total final volume estimate of the eruption ranges from 164.9 km – 247.7 km DRE. This value corresponds to a mass of 4.30 - 6.46 x 10 kg, a magnitude (M) of 7.7 and a VEI of 7. This M makes the CI the largest-magnitude Quaternary eruption in the Mediterranean area. The new detailed estimation of CI eruption physical parameters confirms this event as capable of having significantly affected human activity and the environment on a large scale at the time of the eruption.

A new automated method to retrieve charge layer polarity from flashes, named Chargepol, is presented in this paper. Using data from the NASA Lightning Mapping Array (LMA) deployed during the RELAMPAGO field campaign in Cordoba, Argentina, from November 2018 to April 2019, this method estimates the polarity of vertical charge distributions and their altitudes and thicknesses (or vertical depth) using the very-high frequency (VHF) source emissions detected by LMAs. When this method is applied to LMA data for extended periods of time, it is capable of inferring a storm’s bulk electrical charge structure throughout its life cycle. This method reliably predicted the polarity of charge within which lightning flashes propagated and was validated in comparison to methods that require manual assignment of polarities via visual inspection of VHF lightning sources. Examples of normal and anomalous charge structures retrieved using Chargepol for storms in Central Argentina during RELAMPAGO are presented for the first time. Application of Chargepol to five months of LMA data in Central Argentina and several locations in the United States allowed for the characterization of the charge structure in these regions and for a reliable comparison using the same methodology. About 13.3% of Cordoba thunderstorms were defined by an anomalous charge structure, slightly higher than in Oklahoma (12.5%) and West Texas (11.1%), higher than Alabama (7.3%), and considerably lower than in Colorado (82.6%). Some of the Cordoba anomalous thunderstorms presented enhanced low-level positive charge, a feature rarely if ever observed in Colorado thunderstorms.

Snowpack accumulation in forested watersheds depends on the amount of snow intercepted in the canopy and its partitioning into sublimation, unloading, and melt. A lack of canopy snow measurements limits our ability to evaluate models that simulate canopy processes and predict snowpack and water supply. Here, we tested whether monitoring changes in wind-induced tree sway can enable snow interception detection and estimation of canopy snow water equivalent (SWE). We monitored hourly tree sway across six years based on 12 Hz accelerometer observations on two subalpine conifer trees in Colorado. We developed an approach to distinguish changes in sway frequency due to thermal effects on tree rigidity versus intercepted snow mass. Over 60% of days with canopy snow had a sway signal in the range of possible thermal effects. However, when tree sway decreased outside the range of thermal effects, canopy snow was present 93-95% of the time, as confirmed with classifications of PhenoCam imagery. Using sway tests, we converted significant changes in sway to canopy SWE, which was correlated with total snowstorm amounts from a nearby SNOTEL site (Spearman r=0.72 to 0.80, p<0.001). Greater canopy SWE was associated with storm temperatures between -7 C and 0 C and wind speeds less than 4 m/s. Lower canopy SWE prevailed in storms with lower temperatures and higher wind speeds. We conclude that monitoring tree sway is a viable approach for quantifying canopy SWE, but challenges remain in converting changes in sway to mass and further distinguishing thermal and mass effects on tree sway.

Open cracks and cavities play important roles in fluid transport. Underground water penetration induces microcrack activity, which can lead to rock failure and earthquake. Fluids in cracks can affect earthquake generation mechanisms through physical and physicochemical effects. Methods for characterizing the crack shape and water saturation of underground rock are needed for many scientific and industrial applications. The ability to estimate the status of cracks by using readily observable data such as elastic-wave velocities would be beneficial. We have demonstrated a laboratory method for estimating the crack status inside a cylindrical rock sample based on a vertically cracked transversely isotropic solid model by using measured P- and S-wave velocities and porosity derived from strain data. During injection of water to induce failure of a stressed rock sample, the crack aspect ratio changed from 1/400 to 1/160 and the degree of water saturation increased from 0 to 0.6. This laboratory-derived method can be applied to well-planned observations in field experiments. The in situ monitoring of cracks in rock is useful for industrial and scientific applications such as the sequestration of carbon dioxide and other waste, induced seismicity, and measuring the regional stress field.

An analysis of concurrent extreme events of continental precipitation and Integrated Water Vapour Transport (IVT) is crucial to our understanding of the role of the major global mechanisms of atmospheric moisture transport, including that of the landfalling Atmospheric Rivers (ARs) in extratropical regions. For this purpose, gridded data on CPC precipitation and ERA-5 IVT at a spatial resolution of 0.5 º were used, covering the period from Winter 1980/1981 to Autumn 2017. For each season, and for each point with more than 400 non-dry days, several copula models were fitted to model the joint distribution function of the two variables. At each of the analysed points, the best copula model was used to estimate the probability of a concurrent extreme. At the same time, within the sample of observed concurrent extremes, the proportion of days with landfalling ARs was calculated for the whole period and for two 15-year subperiods, one earlier period and one more recent (warmer) period. Three metrics based on copulas were used to analyse carefully the influence of IVT on extreme precipitation in the main regions of occurrence of AR landfall. The results show that the probability of occurrence of concurrent extremes is strongly conditioned by the dynamic component of the IVT, the wind. The occurrence of landfalling ARs accounts for most of the concurrent extreme days of IVT and continental precipitation, with percentages of concurrent extreme days close to 90% in some seasons in almost all the known regions of maximum occurrence of landfalling ARs, and with percentages greater than 75% downwind of AR landfall regions. This coincidence was lower in tropical regions, and in monsoonal areas in particular, with percentages of less than 50%. With a few exceptions, the role of landfalling ARs as drivers of concurrent extremes of IVT and continental precipitation tends to show a decrease in recent (warmer) periods. For almost all the landfalling AR regions with high or very high probabilities of achieving a concurrent extreme, there is a general trend towards a lower influence of IVT on extreme continental precipitation in recent (warmer) periods.

Oceanic quantities of interest (QoIs), e.g., ocean heat content or transports, are often inaccessible to direct observation, due to the high cost of instrument deployment and logistical challenges. Therefore, oceanographers seek proxies for undersampled or unobserved QoIs. Conventionally, proxy potential is assessed via statistical correlations, which measure covariability without establishing causality. This paper introduces an alternative method: quantifying dynamical proxy potential. Using an adjoint model, this method unambiguously identifies the physical origins of covariability. A North Atlantic case study illustrates our method within the ECCO (Estimating the Circulation and Climate of the Ocean) state estimation framework. We find that wind forcing along the eastern and northern boundaries of the Atlantic drives a basin-wide response in North Atlantic circulation and temperature. Due to these large-scale teleconnections, a single subsurface temperature observation in the Irminger Sea informs heat transport across the remote Iceland-Scotland ridge (ISR), with a dynamical proxy potential of 19%. Dynamical proxy potential allows two equivalent interpretations: Irminger Sea subsurface temperature (i) shares 19% of its adjustment physics with ISR heat transport; (ii) reduces the uncertainty in ISR heat transport by 19% (independent of the measured temperature value), if the Irminger Sea observation is added without noise to the ECCO state estimate. With its two interpretations, dynamical proxy potential is simultaneously rooted in (i) ocean dynamics and (ii) uncertainty quantification and optimal observing system design, the latter being an emerging branch in computational science. The new method may therefore foster dynamics-based, quantitative ocean observing system design in the coming years.

Subsurface sequestration of carbon dioxide (CO2) requires long-term monitoring of the injected CO2 plume to prevent CO2 leakage along the wellbore or across the caprock. Accurate knowledge of the location and movement of the injected CO2 is crucial for risk management at a geological CO2-storage complex. Conventional methods for locating/assessing the injected CO2 plume in the subsurface assume a geophysical model, which is specific and may not be applicable to all types of CO2-injection reservoirs and scenarios. We developed an unsupervised-learning-based visualization of the subsurface CO2 plume that adapts and scales based on the data without requiring an assumption of the geophysical model. The data-processing workflow was applied to the cross-well tomography data from the SECARB Cranfield carbon geo-sequestration project. A multi-level clustering approach was developed to account for data imbalance due to the absence of CO2 in the large portion of the imaged reservoir. The first level of clustering differentiated CO2-bearing regions from the non-CO2 bearing regions and achieved a silhouette score of 0.85, a Calinski-Harabasz index of 160666, and a Davies-Bouldin index of 0.43, which are indicative of high quality, reliable clustering. The second level of clustering further differentiated the CO2-bearing regions into regions containing low, medium, and high CO2 content. Overall, the multi-level clustering achieved a silhouette score, Calinski-Harabasz index, and Davies-Bouldin index of 0.74, 59656, and 0.32, which confirm the high quality and reliability of the newly proposed unsupervised-learning-based visualization. Three distinct clustering techniques, namely k-means, mean-shift, and agglomerative, generated similar visualizations. In terms of the adjusted Rand index, the similarity of clusters identified by the three distinct clustering techniques is around 0.98, which indicates the robustness of the cluster labels assigned to various regions of the CO2-injection reservoir. Further, we find certain geophysical signatures, such as Fourier transform and wavelet transform, to be highly relevant and informative indicators of the spatial distribution of CO2 content.

An ocean iodine cycling model is presented, which predicts upper ocean iodine speciation. The model comprises a three-layer advective and diffusive ocean circulation model of the upper ocean, and an iodine cycling model embedded within this circulation. The two primary reservoirs of iodine are represented, iodide and iodate. Iodate is reduced to iodide in the mixed layer in association with primary production, linked by an iodine to carbon (I:C) ratio. A satisfactory model fit with observations cannot be obtained with a globally constant I:C ratio, and the best fit is obtained when the I:C ratio is dependent on sea surface temperature, increasing at low temperatures. Comparisons with observed iodide distributions show that the best model fit is obtained when oxidation of iodide back to iodate is associated with mixed layer nitrification. Sensitivity tests, where model parameters and processes are perturbed, reveal that primary productivity, mixed layer depth, oxidation, advection, surface fresh water flux and the I:C ratio all have a role in determining surface iodide concentrations, and the timescale of iodide in the mixed layer is sufficiently long for non-local processes to be important. Comparisons of the modelled iodide surface field with parameterisations by other authors shows good agreement in regions where observations exist, but significant differences in regions without observations. This raises the question of whether the existing parameterisations are capturing the full range of processes involved in determining surface iodide, and shows the urgent need for observations in regions where there are currently none.

The Antarctic Ice Sheet (AIS) response to past warming consistent with the 1.5–2°C ‘safe limit’ of the United Nations Paris Agreement is currently not well known. Empirical evidence from the most recent comparable period, the Last Interglaciation, is sparse, and transient ice-sheet model simulations are few and inconsistent. Here we present new results from transiently-forced ice-sheet modelling experiments. We evaluate our results against near and far-field proxy reconstructions and find good agreement. Our simulations indicate that the AIS contributed approximately 4 m to global mean sea level, peaking at 126 ka BP, with ice lost primarily from the Amundsen but not Ross or Weddell Sea sectors. The AIS thinned in the area of the Wilkes Subglacial Basin but did not retreat. Continuing beyond present day our model predicts that the West Antarctic Ice Sheet may already be predisposed to collapse even in the absence of further environmental change.

Global groundwater volumes in the upper 2 km of the Earth’s continental crust – critical for water security – are well estimated. Beyond these depths, a vast body of largely saline and non-potable groundwater exists down to at least 10 km —a volume that has not yet been quantified reliably at the global scale. Here, we estimate the amount of groundwater present in the upper 10 km of the Earth’s continental crust by examining the distribution of sedimentary and cratonic rocks with depth and applying porosity-depth relationships. We demonstrate that groundwater in the 2-10 km zone (what we call ‘deep groundwater’) has a volume comparable to that of groundwater in the upper 2 km of the Earth’s crust. These new estimates make groundwater the largest continental reservoir of water, ahead of ice sheets, provide a basis to quantify geochemical cycles, and constrain the potential for large-scale isolation of waste fluids.

A volcanic eruption is usually preceded by seismic precursors, but their interpretation and use for forecasting the eruption onset time remain a challenge. Eruption processes in geysers are similar to volcanoes, but occur more frequently. Therefore, geysers are useful sites for testing new forecasting methods. We tested the application of Permutation Entropy (PE) as a robust method to assess the complexity in seismic recordings of the Strokkur geyser, Iceland. Strokkur features several minute-long eruptive cycles, enabling us to verify in 63 recorded cycles whether PE behaves consistently from one eruption to the next one. We performed synthetic tests to understand the effect of different parameter settings in the PE calculation. Our application to Strokkur shows a distinct, repeating PE pattern consistent with previously identified phases in the eruptive cycle. We find a systematic increase in PE within the last 15s before the eruption, indicating that an eruption will occur. We quantified the predictive power of PE, showing that PE performs better than seismic signal strength or quiescence when it comes to forecasting eruptions.

Seismic anisotropy produced by aligned olivine in oceanic lithosphere offers a window into mid-ocean ridge dynamics and the state of Earth’s upper mantle. Yet, interpreting anisotropy in the context of grain-scale deformation processes observed in laboratory and natural olivine samples has proven challenging due to the vast length scale differences. We bridge this observational gap by estimating the first in situ elastic tensor of oceanic lithosphere using compressional- and shear-wavespeed anisotropy observations, with fast azimuth parallel to the fossil-spreading direction. This observation is compared with a database of 123 petrofabrics from the literature to infer olivine crystallographic orientations and shear strain accumulated within the lithosphere. Findings indicate D-type olivine lattice-preferred orientation (LPO) with girdled [010] and [001] crystallographic axes and strain accumulation of 300–400%, challenging the prevailing assumption of A-type LPO. This LPO is consistent with olivine deformation during seafloor spreading via grain-size sensitive dislocation-accommodated grain boundary sliding (disGBS), rather than grain-size insensitive dislocation creep. Such deformation implies in situ grain sizes of 0.3–15 mm, smaller on average than steady-state predictions for pure olivine from laboratory calibrations, which may be attributed to grain-boundary pinning by secondary phases or an underestimate of the importance of disGBS by standard flow laws. This work demonstrates the ability to integrate laboratory- with seismological-scale observations of seismic anisotropy and provides new constraints on in situ grain size and associated rheology during near-ridge mantle deformation.