Yan Zhan integrates both observations and models to understand the mechanism driving volcanic unrests and eruptions. By using data assimilation technique, Yan aims to forecast a volcano’s move according to its previous behavior like a weather forecast.

Photo by Pixabay


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

Finite Element Model of Korovin Volcano

Numerical Models

Finite Elemenet Models quatify a volcanos’ 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 utilize models to optimize the stratergy of data collection

Numerical models for volcanic unrest

Before most volcanic eruptions, restless behaviors are usually observed, such as degassing, surface deformation, and earthquakes. It is important to use quantitative tools to understand the magma dynamics and its response of surrounding rocks. Therefore, we can forecast when a magma system is approaching the critical condition for eruption.

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 veritical diking event. The finite element model is calculated by COMSOL Multiphysics 5.5. (Zhan et al., in prep.)

As a postdoctoral fellow at Carnegie, Yan will work with Diana Roman and Hélène Le Mével on exploring the mechanism of precursory deformation, seismicity, and other behavior during volcanic unrest. Yan is developing numerical models based on the laws of poro-viscoelasticity, damage mechanics, and thermodynamics to quantify:

  1. deformation and stress buildup due to magma or volatile dynamics,
  2. failure initiation due to stress accumulation,
  3. stress evolution and failure propagation due to progressive damage,
  4. earthquake triggering due to failure processes,
  5. fracture healing after volcanic unrest, which affects the initial state for the next unrest cycle.

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

%d bloggers like this: