Cross-correlation analysis was applied to long-term onshore broadband records from April 2004 to March 2021 to detect and relocate shallow very low frequency earthquakes (VLFEs) southeast off the Kii Peninsula, along the Nankai Trough, Japan. We then determined the moment rate functions of detected shallow VLFEs using the Monte Carlo-based simulated annealing method. According to this new comprehensive catalog, shallow VLFEs are widespread beneath the accretionary prism toe, but shallow VLFEs with large cumulative moments are localized around the western edge of the paleo-Zenisu ridge, which is subducted beneath southeast off the Kii Peninsula. Our results from long-term shallow VLFE catalog are well consistent with previous studies in this region, suggesting that heterogeneous structures and stress conditions due to the subducted paleo-Zenisu ridge promote the occurrence of shallow slow earthquakes. The relocated shallow VLFE epicenters illustrated three major episodes characterized by a similar activity area and five minor episodes characterized by different areas. The three major episodes exhibited slow frontal migration with different initiation locations, directions, and speeds, as well as several rapid reverse migrations. Episodes of minor activity were distributed in different locations within part of the area of major activity. Different patterns of shallow VLFE migration could reflect temporal changes in the pore-fluid distribution or stress conditions of the plate boundary.

Crustal-stored magma reservoirs contain exsolved volatiles which accumulate in the reservoir roof, exerting a buoyancy force on the crust. This produces surface uplift and sudden loss of volatiles through eruption results in syn-eruptive subsidence. Here, we present three-dimensional, visco-elasto-plastic, numerical modeling results which quantify the ground deformation arising from the growth and release of a volatile reservoir. Deformation is independent of crustal thermal distribution and volatile reservoir shape, but is a function of volatile volume, density and depth and crustal rigidity. We present a scaling law for the volatiles’ contribution to syn-eruptive subsidence and show this contributes ~20% of the observed subsidence associated with the 2015 Calbuco eruption. Our results highlight the key role that volatile-driven buoyancy can have in volcano deformation, show a new link between syn-eruptive degassing and deflation, and highlight that shallow gas accumulation and release may have a major impact on ground deformation of volcanoes.

Temporal changes in the subsurface seismic velocity structure reflect the physical processes that modulate the properties of the media through which seismic waves propagate. These processes, such as healing of the surface damage zone and deep crustal deformation, are described by similar functional forms and operate on similar timescales, making it difficult to determine which process drives the observed changes. We examine earthquake-induced velocity changes using the measured lag-time time series τ(t) of the repeating earthquake sequences since the 2004 Mw 9.2 Sumatra and 2005 Mw 8.6 Nias earthquakes. The S-coda velocity changes (δVS, equivalent to -τS) recover steadily during the 2005−2015 period. The Rayleigh wave velocity changes (δVLR, or -τLR) undergo transient recovery, followed by a strong δVLR reduction in late 2007. δVS recovery is most likely driven by deep processes, whereas the temporal breaks in δVLR recovery in 2007 mostly reflect surface damage and healing induced by the strong ground motions of the 2004 Mw 9.2, 2005 Mw 8.6, 2007 Mw 8.4 and Mw 7.9 Bengkulu, and 2008 Mw 7.3 Simeulue earthquakes. The observed differences between the temporal variations in δVS and δVLR can distinguish deep processes from healing of the surface damage zone.

Faults in carbonate rocks show both seismic and aseismic deformation processes, leading to a wide range of slip velocities. We deformed two centimeter-scale cores of Carrara marble at 25°C, under in-situ conditions of stress of 2-3 km depth, and imaged the nucleation and growth of creeping faults using dynamic synchrotron X-ray microtomography with micrometer spatial resolution. The first sample was under a constant confinement of 30 MPa and no pore fluid. The second sample was under a confinement in the range 35-23 MPa, with 10 MPa pore fluid pressure. We increased the axial stress by steps until creep deformation occurred and imaged deformation in 4D during creep. The samples deformed with a steady-state strain rate when the differential stress was constant, a process called creep. However, for both samples, we also observed transient events that include the acceleration of creep, i.e., creep bursts, phenomena similar to slow slip events that occur in continental active faults. During these transient creep events, strain rates increase and correlate in time with strain localization and the development of system-spanning fault networks. In both samples, the acceleration of opening and shearing of microfractures accommodated creep bursts. Using high-resolution time-lapse X-ray micro-tomography imaging, and digital image correlation, during triaxial deformation allowed quantifying creep in laboratory faults at sub-grain spatial resolution, and demonstrates that transient creep events (creep bursts) correlate with the nucleation and growth of faults.

The spectra of earthquake waveforms can provide important insight into rupture processes, but the analysis and interpretation of these spectra is rarely straightforward. Here we develop a Bayesian framework that embraces the inherent data and modeling uncertainties of spectral analysis to infer key source properties. The method uses a spectral ratio approach to correct the observed waveform spectra of nearby earthquakes for path and site attenuation. The objective then is to solve for a joint posterior probability distribution of three source parameters -- seismic moment, corner frequency, and high-frequency falloff rate -- for each earthquake in the sequence, as well as a measure of rupture directivity for target events with good azimuthal station coverage. While computationally intensive, this technique provides a quantitative understanding of parameter tradeoffs and uncertainties and allows one to impose physical constraints through prior distributions on all source parameters, which guide the inversion when data is limited. We demonstrate the method by analyzing in detail the source properties of 14 different target events of magnitude M5 in southern California that span a wide range of tectonic regimes and fault systems. These prominent earthquakes, while comparable in size, exhibit marked diversity in their source properties and directivity, with clear spatial patterns, depth-dependent trends, and a preference for unilateral directivity. These coherent spatial variations source properties suggest that regional differences in tectonic setting, hypocentral depth or fault zone characteristics may drive variability in rupture processes, with important implications for our understanding of earthquake physics and its relation to hazard.

The Tatun volcano group (TVG) proximate to the 7-million-population metropolis Taipei has long been a central concern in volcanic hazard in the Taiwan. While the TVG has been previously considered an extinct volcano, recent evidences suggested a much younger age of the last eruption event (~6000 years) and the possible existence of magma reservoir beneath the TVG. To examine and unveil the TVG magmatic plumbing system in detail, the local P-wave travel time data and the teleseismic waveform data from a new island-wide Formosa Array Project are combined for a 3D joint inversion. The new model reveals a pronounced magma reservoir with high fraction of melt (~16.5%) beneath the TVG at 5-25 km and better illuminate the deep magmatic structures that may suggest delamination processes in the northern Taiwan. Plain Language Summary The Tatun volcano group (TVG) is close to the metropolis Taipei and poses a direct threat to the 7 million population nearby. While the TVG was thought to be extinct in the past, recent studies found much younger eruption records (~6000 years) and some indications of alive magma chamber underneath. To better understand whether the magma chamber exists and what its geometry is, the Formosa Array project deployed 120 broadband stations in the northern Taiwan well covering the TVG area. Integrating the Formosa Array dataset, a new 3-D P-wave velocity model is constructed and shows a clear image of the magma chamber beneath the TVG at the depth of 5-25 km. Deep magmatic structures are also imaged and suggest the crust in the northern Taiwan was thickened and has possibly peeled off during the mountain building processes.

The sciences struggle to integrate across disciplines, coordinate across data generation and modeling activities, produce connected open data, and build strong networks to engage stakeholders within and beyond the scientific community. The American Geophysical Union (AGU) is divided into 25 sections intended to encompass the breadth of the geosciences. Here, we introduce a special collection of commentary articles spanning 19 AGU sections on challenges and opportunities associated with the use of ICON science principles. These principles focus on research intentionally designed to be Integrated, Coordinated, Open, and Networked (ICON) with the goal of maximizing mutual benefit (among stakeholders) and cross-system transferability of science outcomes. This article 1) summarizes the ICON principles; 2) discusses the crowdsourced approach to creating the collection; 3) explores insights from across the articles; and 4) proposes steps forward. There were common themes among the commentary articles, including broad agreement that the benefits of using ICON principles outweigh the costs, but that using ICON principles has important risks that need to be understood and mitigated. It was also clear that the ICON principles are not monolithic or static, but should instead be considered a heuristic tool that can and should be modified to meet changing needs. As a whole, the collection is intended as a resource for scientists pursuing ICON science and represents an important inflection point in which the geosciences community has come together to offer insights into ICON principles as a unified approach for improving how science is done across the geosciences and beyond.

Supercritical CO2 (ScCO2) invades oilwell cement under geological CO2 sequestration conditions. With the penetration of ScCO2, cement structure prone to damage when the coupled effects of capillarity and carbonation were found. Microstructural evolution of oilwell cement samples was investigated by the CT scanning and the quantitative image-based analysis and show that ScCO2 with the high humid condition would penetrate much deeper than the dry ScCO2 because of the capillarity effects. Due to the deep saline condition in the sequestration formation, the penetration of ScCO2 was retarded by the salt deposition, comparing with the ultrapure water (UP water) conditions. For further assessment of this coupled mechanism, the permeability property and contact angle changes were proposed to analyse the interface region between ScCO2, saline/UP water and oilwell cement.

Jupiter possesses the strongest magnetic field of all planets in the solar system. Unique information about the dynamo process acting at Jupiter can be inferred by modelling and interpreting its field. Using the fluxgate magnetometer measurements acquired during the four years of the Juno mission, we derive an internal and secular magnetic field model in spherical harmonics. The static part is derived to degree 16 with a secular time variation to degree 8. We use properties of the power spectrum of the static field to infer the upper boundary of the dynamo convective region at 0.830±0.022 Jupiter radius. This confirms the role of the transition layer in the field generation inside Jupiter. The secular variation timescales indicate that advective effects dominate the dynamo and the secular variation structures estimated at the dynamo radius suggest that the complex flow involves non-zonal features.

We characterized the texture, composition, and seismic properties of the lithospheric mantle atop the Hawaiian plume by petrostructural analysis of 48 spinel-peridotite xenoliths from four localities in three Hawaiian islands. Coarse-porphyroclastic peridotites with variable degrees of recrystallization, recorded by growth of strain-free neoblasts onto the deformed microstructure, predominate. Full evolution of this process produced equigranular microstructures. Some peridotites have coarse-granular microstructures. Coarse-granular and coarse-porphyroclastic peridotites have strong orthorhombic or axial-[100] olivine crystal-preferred orientations (CPO). Recrystallization produced moderate dispersion and, locally, changed the olivine CPO towards axial-[010]. Enrichment in pyroxenes relative to model melting trends and pyroxenes with interstitial shapes and CPO uncorrelated with the olivine CPO suggest refertilization by reactive melt percolation. The unusual spatial distribution of the recrystallized fraction, Ti-enrichment, and REE-fractionation in recrystallized, equigranular, and coarse-granular peridotites support that these microstructures are produced by static recrystallization triggered by melt percolation. However, there is no simple relation between microstructure and chemical or modal composition. This, together with marked variations in mineral chemistry among samples, implies multiple spatially heterogeneous melt-rock reaction events. We interpret the coarse-porphyroclastic microstructures and CPO as representative of the original oceanic lithosphere fabric. Annealing changed the microstructure to coarse-granular, but did not modify significantly the olivine CPO. Recrystallization produced moderate dispersion of the CPO. “Normal” oceanic lithosphere seismic anisotropy patterns are therefore preserved. Yet Fe-enrichment, refertilization, and limited heating of the base of the lithosphere may reduce seismic velocities by up to 2%, partially explaining negative velocity anomalies imaged at lithospheric depths beneath Hawaii.

The goal of seismic inversion is to obtain subsurface properties from surface measurements. Seismic images have proven valuable, even crucial, for a variety of applications, including subsurface energy exploration, earthquake early warning, carbon capture and sequestration, estimating pathways of sub-surface contaminant transport, etc. These subsurface properties (such as wave speed, density, and elastic velocities) influence the transmission of seismic waves through the subsurface media, and well-understood physics models (so-called “forward models”) can be used to predict what surface measurements would be made, for any given subsurface configuration. Seismic inversion is the inverse problem: given actual surface measurements, infer what subsurface configuration would give rise to those measurements. Like most inverse problems, seismic imaging is ill-posed, meaning many different subsurface configurations can give rise to the same surface measurements. Iterative optimization algorithms for the inverse problem are typically very computationally expensive because they require many evaluations of the forward model, which is itself computationally expensive. A further challenge is the different sensitivity of subsurface properties to the seismic data; density for example is more difficult to accurately infer than is P-wave velocity. But recent advances in algorithms and computing provide an opportunity for remarkable progress in seismic inversion, and efficient solutions to previously infeasible problems have been obtained using data-driven approaches (such as the deep learning methods that were developed primarily for problems in computer vision). The excellent performance of learning-based methods arises from its ability to exploit large amounts of high-quality training data, without the need for hand-designed features. Unlike computer vision, however, seismic inversion is not a data-rich domain. There is a relatively small amount of field data in existence due to the high cost of acquisition, and as a result of its commercial value, a very limited amount is publicly available. To alleviate the data scarcity issue and improve model generalization, there has been growing interest in combining physics knowledge with machine learning for solving seismic inversion problems. This review will survey methods for incorporating physics knowledge with machine learning (primarily deep neural networks) to solve computational seismic inversion problems. We will provide a structured framework of the existing research in the seismic inversion community, and will identify technical challenges, insights, and trends.

On June 7th, 2021 the Juno spacecraft visited Ganymede and provided the first in situ observations since Galileo’s last flyby in 2000. The measurements obtained along a one-dimensional trajectory can be brought into global context with the help of three-dimensional magnetospheric models. Here we apply the magnetohydrodynamic model of Duling et al. (2014) to conditions during the Juno flyby. In addition to the global distribution of plasma variables we provide mapping of Juno’s position along magnetic field lines, Juno’s distance from closed field lines and detailed information about the magnetic field’s topology. We find that Juno did not enter the closed field line region and that the boundary between open and closed field lines on the surface matches the poleward edges of the observed auroral ovals. To estimate the sensitivity of the model results, we carry out a parameter study with different upstream plasma conditions and other model parameters.

The North Anatolian Fault (NAF) has produced numerous major earthquakes. After decades of quiescence, the Mw 6.8 Elazig earthquake (January 24, 2020) has recently reminded us that the East Anatolian Fault (EAF) is also capable of producing significant earthquakes. To better estimate the seismic hazard associated with these two faults, we jointly invert Interferometric Synthetic Aperture Radar (InSAR) and GPS data to image the spatial distribution of interseismic coupling along the eastern part of both the North and East Anatolian Faults. We perform the inversion in a Bayesian framework, enabling to estimate uncertainties on both long-term relative plate motion and coupling. We find that coupling is high and deep (0-20 km) on the NAF and heterogeneous and superficial (0-5 km) on the EAF. Our model predicts that the Elazig earthquake released between 200 and 250 years of accumulated moment, suggesting a bi-centennial recurrence time.

Microseismicity, induced by the sliding of a glacier over its bed, can be used to characterize frictional properties of the ice-bed interface, which are a key parameter controlling ice stream flow. We use naturally occurring seismicity to monitor spatiotemporally varying bed properties at Rutford Ice Stream, West Antarctica. We locate 230000 micro-earthquakes with local magnitudes from –2.0 to –0.3 using 90 days of recordings from a 35-station seismic network located ~40 km upstream of the grounding line. Events exclusively occur near the ice-bed interface and indicate predominantly flow-parallel stick-slip. They mostly lie within a region of interpreted stiff till and along the likely stiffer part of mega-scale glacial landforms. Within these regions, micro-earthquakes occur in spatially (<100 m radius) and temporally (mostly 1-5 days activity) restricted event-clusters (up to 4000 events), which exhibit an increase, followed by a decrease, in event magnitude with time. This may indicate event triggering once activity is initiated. Although ocean tides modulate the surface ice flow velocity, we observe little periodic variation in overall event frequency over time and conclude that water content, bed topography and stiffness are the major factors controlling microseismicity. Based on variable rupture mechanisms and spatiotemporal characteristics, we suggest the event-clusters relate to three end-member types of bed deformation: (1) continuous creation and seismogenic destruction of small-scale bed-roughness, (2) ploughed clasts and (3) flow-oblique deformation during landform-formation or along bedrock outcrops. This indicates that multiple processes, simultaneously active during glacial sliding, can accommodate stick-slip behaviour and that the bed continuously reorganizes.

Short term aftershock incompleteness (STAI) can strongly bias any analysis built on the assumption that seismic catalogs have a complete record of events. Despite several attempts to tackle this issue, we are far from trusting any dataset in the immediate future of a large shock occurrence. Here we introduce RESTORE (REal catalogs STOchastic REplenishment), a Python toolbox implementing a stochastic gap-filling method, which automatically detects the STAI gaps and reconstructs the missing events in the space-time-magnitude domain. The algorithm is based on empirical earthquake properties and relies on a minimal number of assumptions about the data. Through a numerical test, we show that RESTORE returns an accurate estimation of the number of missed events and correctly reconstructs their magnitude, location and occurrence time. We also conduct a real-case test, by applying the algorithm to the Mw 6.2 Amatrice aftershocks sequence. The STAI-induced gaps are filled and missed earthquakes are restored in a way which is consistent with data. RESTORE, which is made freely available, is a powerful tool to tackle the STAI issue, and will hopefully help to implement more robust analyses for advancing operational earthquake forecasting and seismic hazard assessment.

Modelling the planetary heat transport of small bodies in the early Solar System allows us to understand the geological context of meteorite samples. Conductive cooling in planetesimals is controlled by thermal conductivity, heat capacity, and density, which are functions of temperature (T). We investigate if the incorporation of the T-dependence of thermal properties and the introduction of a non-linear term to the heat equation could result in different interpretations of the origin of different classes of meteorites. We have developed a finite difference code to perform numerical models of a conductively cooling planetesimal with T-dependent properties and find that including T-dependence produces considerable differences in thermal history, and in turn the estimated timing and depth of meteorite genesis. We interrogate the effects of varying the input parameters to this model and explore the non-linear T-dependence of conductivity with simple linear functions. Then we apply non-monotonic functions for conductivity, heat capacity and density fitted to published experimental data. For a representative calculation of a 250 km radius pallasite parent body, T-dependent properties delay the onset of core crystallisation and dynamo activity by ~40 Myr, approximately equivalent to increasing the planetary radius by 10 %, and extend core crystallisation by ~3 Myr. This affects the range of planetesimal radii and core sizes for the pallasite parent body that are compatible with paleomagnetic evidence. This approach can also be used to model the T-evolution of other differentiated minor planets and primitive meteorite parent bodies and constrain the formation of associated meteorite samples.

The November 2017 magnitude (MW) 5.5 Pohang, Korea, earthquake, induced by an Engineered 24 Geothermal Systems (EGS) project, caused one fatality and ~US$300M of economic 25 consequences. It has been proposed that a significant probability of such losses was predictable 26 beforehand, from the small earthquakes caused by well-stimulation, so the project should have 27 been suspended, implying that its developer was remiss for not doing so. This argument depends 28 on the low (~0.61) estimated b-value of this earthquake population. However, it is shown that 29 many of the magnitude determinations are inaccurate (underestimated) and the true b-value is 30 higher (1.12 for one subset). The probability of any earthquake as large as MW=5.5, predicted 31 beforehand by extrapolation, was thus much lower than has been claimed. This analysis 32 highlights the necessity of taking care over accuracy when reporting datasets like this, especially 33 in situations where such analyses might influence criminal trials of EGS developers.

We present a coupled finite element spectral boundary integral framework for modeling sequences of earthquakes and aseismic slip on a 2-D planar rate-and-state fault with off-fault visco-plastic response in the plane strain approximation. The model resolves both slow aseismic deformation and inertia effects during rapid slip. As an application of the method, we perform two sets of simulations with different choices of cohesion to explore the co-evolution of fault slip, bulk plasticity and local stress fields. The first set implements a relatively large value of the cohesion parameter, which results in limiting inelastic strain accumulation to dynamic rupture phases. The second set implements a smaller cohesion, allowing for plastic strain to accumulate in both seismic and aseismic phases. For the first model, our results indicate that the extent and distribution of plastic strain depend on the angle of maximum compressive principal stress. At larger angles, inelastic strain accumulates on the extensional side of a dynamically propagating rupture. At smaller angles, the extent of plasticity is limited to the compressional side of the domain. At smaller cohesion values, off-fault plasticity may occur during aseismic slip, which alters the nucleation characteristics and earthquake sequence pattern. Furthermore, our results at lower cohesion values indicate that plastic strain accumulation may occur in both the extensional and compressional sides of the off-fault bulk even at higher angles of maximum compression. This produces damage patterns that deviate from the traditional off-fault fan-like distribution observed in dynamic rupture simulations and emphasizes the significance of long-term deformation in interpreting observations.