Research

Integrate observations and numerical models by data assimilation technique, to understand the mechanism driving volcanic unrests, earthquakes,hydrothermal activities, and lithospherical deformation on and beyond the Earth.

Photo by Pixabay

Observations

Surface deformation, earthquakes, gravity changes, gas emissions, & rock records, etc. provide critical information on physical & chemical conditions in a volcanic system

Finite Element Model of Korovin Volcano

Numerical Models

Finite Element Models quantify a volcano’s behavior under a given initial & boundary condition when ruled by solid mechanics, fluid dynamics, & thermodynamics

Image by Altair Engineering

Data Assimilation

Monte Carlo-based Ensemble Kalman Filter not only upgrades numerical models by observations but also utilizes models to optimize the strategy of data collection

Anticipating Earth’s Fiery Breath

Volcanoes, often perceived as silent giants, whisper warnings before unleashing their fiery power. These precursory signals, ranging from subtle gas emissions and ground deformation to more dramatic earthquake swarms, offer a glimpse into the restless heart of our planet.

Understanding the language of volcanoes is crucial for predicting eruptions and mitigating their hazards. This is where the power of quantitative analysis comes into play. By applying sophisticated tools to decipher the subtle changes in magma dynamics and the surrounding rock’s response, we can begin to anticipate a volcano’s next fiery outburst.

Imagine peering beneath the Earth’s surface, tracking the movement of molten rock, and measuring the stress building within the crust. Through these quantitative insights, we can identify when a magma system is approaching a critical tipping point – a crucial step toward forecasting eruptions and safeguarding vulnerable communities.

Unlocking the Secrets of Magma Reservoirs: A Crystal Mush Perspective

Imagine a crystal mush – a slush of molten rock and crystal fragments – deep within the Earth. This unique perspective offers a window into the dynamic heart of a magma reservoir, where powerful forces drive volcanic eruptions. But what exactly are the conditions within this fiery cauldron?

Flow velocity (left) & Crystallinity (right) of a crystal mush initially with the temperature of 1000°C.

Through the power of numerical modeling, we can simulate the evolution of a magma reservoir over time, revealing its intricate dance with the surrounding crust. This allows us to visualize the build-up to an eruption and its impact on the Earth’s surface. By unlocking the secrets held within the crystal mush, we gain valuable insights into the awe-inspiring power of volcanoes.

Degassing and sealing

Temperature (Left), Steam velocity (Middle), and Vertical deformation (Right) from a finite element model of the Augustine Volcano, showing the host rock response to degassing before and after the upper conduit is sealed (Zhan et al., in prep.)

Dike propagation

Surface deformation, minimum principal stress, and pore pressure evolution during a vertical diking event. The finite element model is calculated by COMSOL Multiphysics 5.5. (Zhan et al., in prep.)

Data assimilation

Quantifying the eruption potential of a restless volcano requires the ability to model parameters such as overpressure and calculate the host rock stress state as the system evolves. A critical challenge is developing a model-data fusion framework to take advantage of observational data and provide updates of the volcanic system through time. The Ensemble Kalman Filter (EnKF) uses a Monte Carlo approach to assimilate volcanic monitoring data and update models of volcanic unrest, providing time-varying estimates of overpressure and stress. In Zhan and Gregg (2017), the EnKF was shown to be effective for forecasting volcanic deformation using synthetic InSAR and GPS data.

Illustration of the data assimilation strategy for magma reservoir modeling. The example model shows the predicted surface deformation (top surface) and von Mises stress due to a pressurized magma reservoir. The Ensemble Kalman Filter (EnKF) approach is modified from Gregg and Pettijohn (2016). (A) The initial ensemble of models is produced, including model uncertainties. The forecast ensemble (B) is sequentially derived from forward integration of the initial suite of models. (C) Data utilized in this implementation include InSAR and/or GPS, but could also include additional such as gravity, seismicity, and tomography. (D) When new data are available, an EnKF analysis is conducted to update the model parameters (E) and change the trajectory of the model. The updated model parameters are then used to create a new forecast ensemble. (Zhan and Gregg, 2017)

In Zhan et al. (2017), the EnKF is used to provide a “hindcast” of the 2009 explosive eruption of Kerinci volcano, Indonesia. A two-source analytical model is used to simulate the surface deformation of Kerinci volcano observed by InSAR time-series data and to predict system evolution. A deep, deflating dike-like source reproduces the subsiding signal on the flanks of the volcano, and a shallow spherical McTigue source reproduces the central uplift. EnKF predicted parameters are used in finite element models to calculate the host-rock stress state prior to the 2009 eruption. Models reveal that the host rock around the shallow magma reservoir was trending towards tensile failure prior to 2009, which may have been the catalyst for the 2009 eruption. Another application of the EnKF approach is utilizing high performance computing data assimilation to provide a successful forecast of the 26 June 2018 eruption of Sierra Negra volcano in the Galápagos, more than six months in advance (Gregg and Zhan et al., in prep.). Forecast models combined with satellite InSAR data reveal that the evolution of the stress state in the host rock surrounding the Sierra Negra magma system likely controlled the timing of the eruption.

Click here to learn more.

Volcanic unrest in Laguna del Maule volcano

The Laguna del Maule (LdM) volcanic field is a large rhyolitic magmatic system in the Chilean Andes, which has exhibited frequent eruptions during the past 20 ka. Rapid surface uplift (>20 cm/yr) has been observed since 2007 accompanied by localized earthquake swarms and micro-gravity changes, indicating the inflating magma reservoir may interact with a pre-existing weak zone (i.e., Troncoso fault). In this investigation, we model the magma reservoir by data assimilation with InSAR data. The reservoir geometry is comparable to the magma body inferred by seismic tomography, magnetotelluric, and gravity studies. The models also suggest that a weak zone, which has little effect on surface displacement, can significantly affect the stress field and host rock failure development, especially with high pore-pressures (Fig. 1).  In particular, concentrated dilatancy within the weak zone facilitates the microfracture formation during reservoir inflation. High-pressure fluid can inject into the weak zone from the magma reservoir to trigger earthquakes, and further migrate upwards to create positive gravity changes by occupying unsaturated storages. The pore pressure will then decrease, halting the seismicity swarm until the next cycle. This “hydrofracturing” process may release some accumulated stress along the magma reservoir delaying an eventual eruption in turn. Besides, the resultant models are propagated forward in time to evaluate potential stress trajectories for future unrest.

A 3D cartoon illustrating the interaction between the magma reservoir and hydrothermal system. (a) Microfractures open or get connected due to the inflation of the magma chamber. (b) A high-pressure fluid from the magma reservoir then injects into the dilatant fault zone triggering seismicity. (c) The fluid moves to unsaturated shallow storages to create a positive gravity anomaly, and the pore pressure will decrease after that to stop the seismic swarm until the next cycle. (Zhan et al., 2019)

Click here to learn more.

Effects of material proporties on volcano modeling

The effects of Young’s modulus and Poisson’s ratio on brittle failure with variations in the source depth to center. The CMSU is the maximum value of the surface uplift that can be observed before any failure is initiated around the magma chamber.
Mohr circle diagrams illustrating that the initial failure type is controlled by host rock stiffness and strength. The Mohr circles represent of the stress state at the tip of the chamber as shown in the inserted plot. (a) For the host rock with relatively low Young’s Modulus (E = 20 GPa) and cohesion (20 MPa), the Mohr circle increases and touches the Coulomb failure envelop first. (b) The high cohesion (80 MPa) of the stronger rock (E = 80 GPa) allows the Mohr circle to reach tensile failure first. (Zhan and Gregg, 2019)

Forecasting the onset of a volcanic eruption from a closed system requires understanding its stress state and failure potential, which can be investigated through numerical modeling. However, the lack of constraints on model parameters, especially rheology, may substantially impair the accuracy of failure forecasts. Therefore, it is essential to know whether large variations and uncertainties in rock properties will preclude the ability of models to predict reservoir failure. In Zhan and Gregg (2019), a series of 2-dimensional, axisymmetric models are used to investigate sensitivities of brittle failure initiation to assumed rock properties. The numerical experiments indicate that Poisson’s ratio and internal friction angle have much less effect on failure forecasts than Young’s modulus. Variations in Young’s modulus significantly affect the prediction of surface deformation before failure onset when Young’s modulus is <40 GPa. Earlier precursory volcano-tectonic events may occur in weak host-rock (E <40 GPa) due to an early stage of Coulomb failure prior to tensile failure. Thus, combining surface deformation with seismicity may enhance the accuracy of eruption forecast in these situations.

Click here to learn more

Intraplate earthquake

Geological background of the Central United States. The tectonic provinces and three ancient rifts are modified after Hoffman (1989) and Whitmeyer and Karlstrom (2007). Black arrows show region horizontal compressional stress (Heidbach et al., 2008). Red circles are epicenters of M ≥ 3.0 events from 1811 to 2015
A geodynamic model to explain the process of earthquake triggering due to stress accumulation. A wedge-like mantle low-velocity zone (MLVZ) makes the stress localized above it and a relative strong Lower Crust exhibiting high seismic velocity (LCHV) enables the stress building beneath the ancient Reelfoot Rift. The reactivation of the pre-existing fault within the rift triggers earthquakes in the NMSZ. (Zhan et al., 2016)

The mechanism of intraplate earthquakes is not well understood. For example, the New Madrid Seismic Zone (NMSZ) in the Midwestern United States was the site of several major M 6.8–8 earthquakes in 1811–1812 and remains seismically active, but the ultimate controls on earthquake initiation and the duration of the seismicity remain unclear. In Zhan et al. (2016), we develop a series of numerical experiments to determine the impact of lithospheric heterogeneity on earthquake nucleation and rupture processes (Fig. 2). The rheological structure of the lithosphere is inferred by seismic tomography. Results indicate that due to a relatively weak upper mantle, the differential stress prefer to accumulate beneath the crust of NMSZ, whereas the stresses below the geologically similar Midcontinent Rift System are comparatively low. The numerical observations coincide with the observed distribution of seismicity throughout the region. Furthermore, the relatively strong crust in this region, exhibited by high seismic velocities, enables the elevated stress to extend to the base of the ancient rift system, reactivating fossil rifting faults and therefore triggering earthquakes.

Click here to learn more