From August 2018 to May 2019, Kīlauea’s summit exhibited unique, simultaneous, inflation and deflation, apparent in both GPS time series and cumulative InSAR displacement maps. This deformation pattern provides clear evidence that the Halema‘uma‘u (HMM) and South Caldera (SC) reservoirs are distinct. Post-collapse inflation of the East Rift Zone (ERZ), as captured by InSAR, indicates concurrent magma transfer from the summit reservoirs to the ERZ. We present a physics-based model that couples pressure-driven flow between these magma reservoirs to simulate time dependent summit deformation. We take a two-step approach to quantitatively constrain Kīlauea’s magmatic plumbing system. First, we jointly invert the InSAR displacement maps and GPS offsets for the location and geometry of the summit reservoirs, approximated as spheroidal chambers. We find that HMM reservoir has an aspect ratio of ~1.8 (prolate) and a depth of ~2.2 km (below surface). The SC reservoir has an aspect ratio of ~0.14 (oblate) and a depth of ~3.6 km. Second, we utilize the flux model to invert GPS time series from 8 summit stations. Results favor a shallow HMM-ERZ pathway an order of magnitude more hydraulically conductive than the deep SC-ERZ pathway. Further analysis shows that the HMM-ERZ pathway is required to explain the deformation time series. Given high-quality geodetic data, such an approach promises to quantify the connectivity of magmatic pathways between reservoirs in other similar volcanic systems.
Improved understanding of the impact of crystal mush rheology on the response of magma chambers to magmatic events is critical for better understanding crustal igneous systems with abundant crystals. In this study, we extend an earlier model by (Liao et al, 2018) which considers the mechanical response of a magma chamber with poroelastic crystal mush, by including poroviscoelastic rheology of crystal mush. We find that the coexistence of the two mechanisms of poroelastic diffusion and viscoelastic relaxation causes the magma chamber to react to a magma injection event with more complex time-dependent behaviors. Specifically, we find that the system’s short-term evolution is dominated by the poroelastic diffusion process, while its long-term evolution is dominated by the viscoelastic relaxation process. We identify two post-injection timescales that represent these two stages and examine their relation to the material properties of the system. We find that better constraints on the poroelastic diffusion time are more important for the potential interpretation of surface deformation using the model. We also find that the combination of the two mechanisms causes magma transport to reverse direction in the system, which would successively expose crystals to magma with different chemical compositions.
On January 15, 2022, Tonga’s Hunga Tonga-Hunga Ha’apai (HTHH) volcano violently erupted, generating a tsunami that killed three people. Acoustic-gravity waves propagated by the eruption and tsunami caused global complex ionospheric disturbances. In this paper, we study the nature of these perturbations from Global Navigation Satellite System observables over the southwestern Pacific. After processing data from 818 ground stations, we detect supersonic acoustic waves, Lamb waves, and tsunamis, with filtered magnitudes between 1 and 7 Total Electron Content units. Phase arrivals appear superpositioned up to ~1000 km from HTHH and are distinct by ~2200 km. Within ~2200 km, signals have an initial low-frequency pulse that transitions to higher frequencies. We note the presence of a faster perturbation generated one hour post-eruption which crosses the tsunami disturbance ~3000 km from HTHH, potentially contributing to premature land arrivals. Lastly, the arrival of tsunami-generated disturbances coincides with deep-ocean observations.
Volcanic eruptions are potentially hazardous natural events. Understanding how magma accumulates, migrates, and erupts is important to understanding and, eventually, predicting volcanic eruptions. However, the variation in the scale of volcanoes, co-occurrence of earthquakes, and the duration of the eruption makes understanding these events difficult. Ambient noise interferometry is becoming an increasingly more popular tool to study and monitor active volcanoes. We use this method to characterize the variations of subsurface seismic velocities associated with different stages of the eruption process at the Great Sitkin Volcano in the central Aleutian volcanic arc. This volcano initially erupted in May 2021 with elevated seismicity and gas release, followed by the formation of a new lava dome starting July 2021. The volcano had an increase in seismicity in February 2020 but without any eruption activity reported. Measuring the variation of seismic velocities from August 2019 to March 2022, we observe a local decrease in velocity leading up to the eruption and an increase in velocity following the emplacement of the lava dome. We do not observe any velocity variations preceding the non-eruptive increase of seismic activity in February 2020. Despite its remote location and relatively small scale, the findings of this study at the Great Sitkin volcano have significant implications for understanding volcanism and the development and prediction of volcanic eruptions in general.
Dynamic modeling of sequences of earthquakes and aseismic slip (SEAS) provides a self-consistent, physics-based framework to connect, interpret, and predict diverse geophysical observations across spatial and temporal scales. Amid growing applications of SEAS models, numerical code verification is essential to ensure reliable simulation results but is often infeasible due to the lack of analytical solutions. Here, we develop two benchmarks for three-dimensional (3D) SEAS problems to compare and verify numerical codes based on boundary-element, finite-element, and finite-difference methods, in a community initiative. Our benchmarks consider a planar vertical strike-slip fault obeying a rate- and state-dependent friction law, in a 3D homogeneous, linear elastic whole-space or half-space, where spontaneous earthquakes and slow slip arise due to tectonic-like loading. We use a suite of quasi-dynamic simulations from 10 modeling groups to assess the agreement during all phases of multiple seismic cycles. We find excellent quantitative agreement among simulated outputs for sufficiently large model domains and small grid spacings. However, discrepancies in rupture fronts of the initial event are influenced by the free surface and various computational factors. The recurrence intervals and nucleation phase of later earthquakes are particularly sensitive to numerical resolution and domain-size-dependent loading. Despite such variability, key properties of individual earthquakes, including rupture style, duration, total slip, peak slip rate, and stress drop, are comparable among even marginally resolved simulations. Our benchmark efforts offer a community-based example to improve numerical simulations and reveal sensitivities of model observables, which are important for advancing SEAS models to better understand earthquake system dynamics.
The 2020 Mw 6.8 Elazig earthquake was the largest along the Eastern Anatolian Fault (EAF) in over a century, providing valuable insights into its rupture behavior. We use satellite geodesy and seismology to detail the mainshock rupture, postseismic deformation and aftershocks. The mainshock propagated mostly westwards at 2 km/s from a nucleation point on an abrupt 10° fault bend. Only one end of the rupture corresponds to an established EAF segment boundary, and the earthquake may have propagated into the slip zone of the 1874 M 7.1 Golcuk Golu earthquake. It exhibits a pronounced (80%) shallow slip deficit, only a small proportion of which is recovered by early aseismic afterslip. The slow rupture velocity, shallow slip deficit and low afterslip are characteristic of earthquakes hosted by faults of low-to-intermediate structural maturity, indicating that faults continue to evolve in important ways even as they accrue cumulative off sets of tens of kilometers
Axial Seamount is a basaltic hot spot volcano with a summit caldera at a depth of ~1500 m below sea level, superimposed on the Juan de Fuca spreading ridge, giving it a robust and continuous magma supply. Axial erupted in 1998, 2011, and 2015, and is monitored by a cabled network of instruments including bottom pressure recorders and seismometers. Since its last eruption, Axial has re-inflated to 85-90% of its pre-eruption level. During that time, we have identified eight discrete, short-term deflation events of 1-4 cm over 1-3 weeks that occurred quasi-periodically, about every 4-6 months between August 2016 and May 2019. During each short-term deflation event, the rate of earthquakes dropped abruptly to low levels, and then did not return to higher levels until reinflation had resumed and returned near its previous high. The long-term geodetic monitoring record suggests that the rate of magma supply has varied by an order of magnitude over decadal time scales. There was a surge in magma supply between 2011-2015, causing those two eruptions to be closely spaced in time and the supply rate has been waning since then. This waning supply has implications for eruption forecasting and the next eruption at Axial still appears to be 4-9 years away. We also show that the number of earthquakes per unit of uplift has increased exponentially with total uplift since the 2015 eruption, a pattern consistent with a mechanical model of cumulative rock damage leading to bulk failure during magma accumulation between eruptions.
A new mission called GRACE Follow-On is now flying to continue the measurements started by the GRACE mission, and to test a laser interferometry system for making more accurate measurements of the satellite separation. In this paper we discuss the potential scientific benefit of strongly reducing the acceleration noise in a Next Generation Gravity Mission (NGGM), compared with that for GRACE and for GRACE Follow-On. A useful way of comparing the scientific benefits is from the view point of how well they can be used to test different procedures for estimating the changes in the geopotential based on sources of geophysical information other than satellite gravity results. In particular, changes in hydrology, the atmospheric density, and ocean conditions can make large and very non-uniform changes in the geopotential in short periods of time. To make the discussion as simple as possible, we consider mainly the variations in the geopotential at altitude along the satellite orbit for different ground tracks. For the NGGM, we initially assume laser interferometry between the two satellites but the same satellite acceleration noise level as for the GRACE-Follow-On mission. Then the total measurement noise level at long and medium wavelengths would be only moderately below the geopotential variation estimation uncertainty. However, if the acceleration noise level were sharply reduced by replacing the GRACE-type accelerometers by simplified gravitational reference sensors, it appears that considerably improved tests of different procedures for geophysical estimates of the geopotential variations could be made.
To date, approximately 20% of the ocean floor has been surveyed by ships at a spatial resolution of 400 m or better. The remaining 80% has depth predicted from satellite altimeter-derived gravity measurements at a relatively low resolution. There are many remote ocean areas in the southern hemisphere that will not be completely mapped at 400 m resolution during this decade. This study is focused on the development of synthetic bathymetry to fill the gaps. There are two types of seafloor features that are not typically well resolved by satellite gravity: abyssal hills and small seamounts (< 2.5 km tall). We generate synthetic realizations of abyssal hills by combining the measured statistical properties of mapped abyssal hills with regional geology including fossil spreading rate/orientation, rms height from satellite gravity, and sediment thickness. With recent improvements in accuracy and resolution, It is now possible to detect all seamounts taller than about 800 m in satellite-derived gravity and their location can be determined to an accuracy of better than 1 km. However, the width of the gravity anomaly is much greater than the actual width of the seamount so the seamount predicted from gravity will underestimate the true seamount height and overestimate its base dimension. In this study we use the amplitude of the vertical gravity gradient (VGG) to estimate the mass of the seamount and then use their characteristic shape, based on well surveyed seamounts, to replace the smooth predicted seamount with a seamount having a more realistic shape.
Bangladesh, a small and over populated country in Southeast Asia occupies most of the Bengal Basin that results from sediments derived from the collision of India with Asia. The basin is filled with a 19 km thick sequence of Cenozoic sediments deposited by the mighty rivers Ganges and Brahmaputra. Unconsolidated Holocene sediments susceptible to seismic amplification characterize the upper part of the Cenozoic sequence. Bangladesh sits a top on three tectonic plates; India, Tibet and Burma. The India plate is colliding with the Tibet subplate to the north, which gives rise to great Himalayas, while to the east it is subducting beneath Burma and Sunda slivers, which gave rise to Indo-Burma arc. The Surma basin of NE Bangladesh is being underthrust under the Shillong massif producing the 2-km high plateau. The Indo-Burma fold and thrust belt results from the oblique subduction of the thick sediments of the Bengal Basin on the India plate that has deformed into a series of north-south trending en-echelon folds and thrust faults. The faults rooting these folds and the underlying megathrust are capable of generating devastating earthquakes in and around Bangladesh. Past earthquakes have brought changes to the landscape, avulsion of rivers Brahmaputra and Meghna, migration of human settlements, and widespread sand liquefactions and sand and/or mud eruptions. Our GPS study demonstrated that the landward extension of Andaman-Sumatra subduction zone into Indo-Burma subduction in deltaic Bangladesh is active. The present day India-Burma oblique convergence rate is 17 mm/y and that the décollement beneath the fold-thrust belt is locked (Steckler et. al., 2016). The western part of the subduction zone over a shallow décollement shows little seismicity whereas the eastern part shows moderate seismicity of magnitude 4 to 6. Based on the GPS velocity across the fold belt and seismicity the Indo-Burma subduction zone can be potentially be divided into locked western segment and slipping eastern segment, analogous to Cascadia subduction zone. Fold belt parallel shortening across Dauki Fault in Shillong is 7 mm/yr, which is another potential source of a large earthquake. The huge population might be severely ravaged by devastating earthquakes from both these sources.
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.
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.
The triggering of large earthquakes on a fault hosting aseismic slip or, conversely, the triggering of slow slip events (SSE) by passing seismic waves involves seismological questions with important hazard implications. Just a few observations plausibly suggest that such interactions actually happen in nature. In this study we show that three recent devastating earthquakes in Mexico are likely related to SSEs, describing a cascade of events interacting with each other on a regional scale via quasi-static and/or dynamic perturbations. Such interaction seems to be conditioned by the transient memory of Earth materials subject to the “traumatic” stressing produced by the seismic waves of the great 2017 (Mw8.2) Tehuantepec earthquake, which strongly disturbed the aseismic slip beating over a 650 km long segment of the subduction plate interface. Our results imply that seismic hazard in large populated areas is a short-term evolving function of seismotectonic processes that are often observable.
On 12 January 2020, Taal volcano, Philippines, erupted after 43 years of repose, affecting more than 500,000 people. Using interferometric synthetic aperture radar (InSAR) data, we present the complete pre- to post-eruption analyses of the deformation of Taal. We find that: 1) prior to eruption, the volcano experienced long-term deflation followed by short-term inflation, reflecting the depressurization-pressurization of its ~5 km depth magma reservoir; 2) during the eruption, the magma reservoir lost a volume of 0.531 +/- 0.004 km^3 while a 0.643 +/- 0.001 km^3 lateral dike was emplaced; and 3) post-eruption analyses reveal that the magma reservoir started recovery approximately 3 weeks after the main eruptive phase. We propose a conceptual analysis explaining the eruption and address why, despite the large volume of magma emplaced, the dike remained at depth. We also report the unique and significant contribution of InSAR data during the peak of the crisis.
From 29 June to 1 July, 2015, a phreatic eruption occurred in Owakudani, the largest fumarole area in Hakone volcano, Japan. In this study, an interferometric synthetic aperture radar (InSAR) time series analysis of the Advanced Land Observing Satellite-2 (ALOS-2)/Phased Array type L-band Synthetic Aperture Radar-2 (PALSAR-2) data was performed to measure deformation after the eruption. The results show that the central cones of the volcano have subsided since the eruption and its deflation source is located beneath the previously estimated bell-shaped conductor, which is considered as a sealing layer confining a pressurized hydrothermal reservoir. Therefore, the InSAR results demonstrate the deflation of the hydrothermal system beneath the volcano. One possible cause of this deflation is compaction due to a decrease in pore pressure caused by rupture and fluid migration during and after the eruption.
Vertical land motion (VLM) from past and ongoing glacial changes can amplify or mitigate ongoing relative sea level change. We present a high resolution VLM-model for the wider Arctic, that includes both present-day ice loading (PDIL) and glacial isostatic adjustment (GIA). The study shows that the non-linear elastic uplift from PDIL is significant (0.5-1 mm/y ) in most of the wider Arctic and exceeds GIA at 15 of 54 Arctic GNSSsites, including sites in non-glaciated areas of the North Sea region and the east coast of North America. Thereby the sea level change from PDIL (1.85 mm/y) is significantly mitigated from VLM caused by PDIL. The combined VLM-model was consistent with measured VLM at 85% of the GNSS-sites (R=0.77) and outperformed a GIA-only model (R=0.64). Deviations from GNSS-measured VLM can be attributed to local circumstances causing VLM.
Different Earth orientation parameter (EOP) time series are publicly available that typically arise from the combination of individual space geodetic technique solutions. The applied processing strategies and choices lead to systematically differing signal and noise characteristics particularly at the shortest periods between 2 and 8 days. We investigate the consequences of typical choices by introducing new experimental EOP solutions obtained from combinations at either normal equation level processed by DGFI-TUM and BKG, or observation level processed by ESA. All those experiments contribute to an effort initiated by ESA to develop an independent capacity for routine EOP processing and prediction in Europe. Results are benchmarked against geophysical model-based effective angular momentum functions processed by ESMGFZ. We find, that a multi-technique combination at normal equation level that explicitly aligns a priori station coordinates to the ITRF2014 frequently outperforms the current IERS standard solution 14C04. A multi-GNSS-only solution already provides very competitive accuracies for the equatorial components. Quite similar results are also obtained from a short combination at observation level experiment using multi-GNSS solutions and SLR from Sentinel-3A and -3B to realize space links. For ΔUT1, however, VLBI information is known to be critically important so that experiments combining only GNSS and possibly SLR at observation level perform worse than combinations of all techniques at normal equation level. The low noise floor and smooth spectra obtained from the multi-GNSS solution nevertheless illustrates the potential of this most rigorous combination approach so that further efforts to include in particular VLBI are strongly recommended.