Forecasting mechanical failure and the 26 June 2018 eruption of Sierra Negra Volcano, Galápagos, Ecuador.

Using recent advancements in high-performance computing data assimilation to combine satellite InSAR data with numerical models, the prolonged unrest of the Sierra Negra volcano in the Galápagos was tracked to provide a fortuitous, but successful, forecast 5 months in advance of the 26 June 2018 eruption. Subsequent numerical simulations reveal that the evolution of the stress state in the host rock surrounding the Sierra Negra magma system likely controlled eruption timing. While changes in magma reservoir pressure remained modest (<15 MPa), modeled widespread Mohr-Coulomb failure is coincident with the timing of the 26 June 2018 moment magnitude 5.4 earthquake and subsequent eruption. Coulomb stress transfer models suggest that the faulting event triggered the 2018 eruption by encouraging tensile failure along the northern portion of the caldera. These findings provide a critical framework for understanding Sierra Negra's eruption cycles and evaluating the potential and timing of future eruptions.

INTRODUCTION

One of the great challenges in the field of volcanology is to develop quantitative models to investigate the processes that lead to volcanic eruptions and use these models to provide eruption forecasts (). Meeting this challenge requires the development of models capable of interpreting field observations and tracking the evolution of a magma system. One widely used observation of volcanic unrest is ground deformation. As magma accumulates in a subsurface reservoir, the overlying ground surface above it swells (). The surface signal, captured by satellite or ground measurements, can be inverted to estimate volume changes in a burried magma source. Deformation signals often begin months to years before an eruption, providing early warning of volcanic activity (). However, linking a surface deformation signal to eruption likelihood is challenging, and, frequently, the magnitude of the signal does not provide an accurate indication of eruption potential or timing (). In particular, volcano inflation does not always lead to an eruption and eruptions can occur when no preceding surface displacement is detected (). Hence, short-term observations, such as changes in seismicity on the time scales of minutes to hours, have typically been more successful predictors of an impending eruption (). However, clear signals are often lacking [e.g., the lack of seismic precursors at many Aleutian volcanoes () or the lack of precursory deformation at open volcanic systems ()], making eruption forecasting difficult.

Multiphysics-based numerical modeling approaches provide a means for investigating eruption catalysts by calculating the evolution of the host rock stress and reservoir pressure during periods of unrest (, ). A critical advancement is combining these models with geophysical observations to track a system’s evolution in real time. Model-data fusion frameworks, often used in climate modeling (), are key for investigating how a system evolves. The ensemble Kalman filter (EnKF), an ensemble-based Markov chain Monte Carlo (MCMC) sequential data assimilation approach (), has recently been adapted for tracking volcano system evolution (–). The EnKF has shown great promise for assimilating large geospatial data, such as satellite interferometric synthetic aperture radar (InSAR) and ground-based global navigation satellite system (GNSS) deformation data, into multiphysics volcano finite element method (FEM) models (). The high-performance computing (HPC), EnKF approach provides updates of the volcano system state through time, including an estimation of host rock stress, failure, and magma reservoir pressure. Estimations of stress and failure provide insight into the stability of a magma system and potential triggering mechanisms for magma migration and eruption, which is particularly beneficial in the absence of clear precursors ().

Sierra Negra, the most voluminous of the Galápagos volcanoes, is a 60 × 40 km basaltic shield volcano that occupies most of the southern portion of Isabela Island (Fig. 1) (). Sierra Negra has experienced at least seven historic eruptions since 1911, with an eruption occurring approximately every 15 years (, ). Sierra Negra’s prolonged intereruption cycle with extensive uplift and seismicity provides a unique natural laboratory to investigate stress evolution of a volcano and potential eruption precursors and catalysts while testing the EnKF approach. Sierra Negra experienced an extended period of unrest before its two most recent eruptions (in 2005 and 2018). Preceding its 2005 eruption, caldera-centered uplift >5 m was observed, culminating in a moment magnitude (Mw) 5.5 earthquake on the trapdoor fault system in the southern caldera floor followed within hours by an eruption on 22 October 2005 (, ). By the spring of 2018, the magnitude of the observed inflation since 2005 had reached 6.5 m (, ). The 26 June 2018 eruption commenced at 1340 Local Time (LT) and was preceded by a rapid increase in seismicity including a Mw 5.4 event that struck at 0315 LT, also along the southern side of the caldera trapdoor fault system (Fig. 1) (, ). The 2018 eruption lasted for 58 days and covered a 30.6-km2 area in fresh lava flows ().

Fig. 1.

Sierra Negra’s 26 June 2018 Mw 5.4 earthquake and eruption.

The eruption commenced on 26 June at 1340 LT from five fissures (white dashed lines, indicated by F1 to F5) with resultant lava flows indicated by red shaded regions (). Black dashed lines indicate the location of the trapdoor fault system on the southern and southwestern portion of the caldera. Focal mechanism shows the location of the Mw 5.4 earthquake.

Several questions surround the protracted unrest periods observed at Sierra Negra, including how its magma system endures such substantial and rapid inflation before eruption. Previous studies have pointed to the release of stress accumulated due to magma intrusion by slip along trapdoor faults in the overlying roof (, –) and the ductile nature of the warm host rock limiting stress accumulation near the reservoir (). These hypothesized mechanisms for stress relief and the rheological buffering of eruption-driving failure accumulation make Sierra Negra both a challenging and attractive target for testing new volcano forecasting techniques.

In this investigation, the EnKF approach is used to evaluate the large-magnitude uplift observed at Sierra Negra volcano leading up to its 26 June 2018 eruption. In January 2018, 5 months before the eruption, an initial EnKF forecast for Sierra Negra was completed using selected Sentinel-1 InSAR observations from 2014 to 2018. The forecast indicated that an eruption was likely to occur between 25 June 2018 and 5 July 2018 due to the significant and widespread accumulation of brittle failure in the host rock surrounding the magma reservoir. After the 2018 eruption, additional EnKF experiments were conducted using additional InSAR deformation data to evaluate the success of the pre-eruption forecast and provide strategies for conducting future forecasts. This paper explores how the January 2018 forecast for Sierra Negra was both a lucky accident and an encouraging sign for transformative advances in volcano forecasting. The primary goals of this study are (i) to evaluate the eruption precursor signals and eruption triggering mechanisms at Sierra Negra and (ii) to assess the ability of the EnKF approach to accurately track and forecast the system state through time.

RESULTS

In the EnKF analysis presented here, an ensemble containing 240 three-dimensional (3D) FEM models was updated sequentially through time as new InSAR observations became available. Each of the 240 models in the ensemble is unique, defined by their parameter values. The initial ensemble of models was generated using a Monte Carlo approach to choose parameter values. Hence, we only prescribed the magma reservoir as a spheroidal geometry but allowed the EnKF to determine the best-fit parameters for the ellipsoid shape (size, prolate versus oblate) and its location (see table S3 for parameter initial values). During each EnKF analysis step, the parameter values for all the models are updated to nudge the ensemble toward a better fit with the observations (see fig. S2 for details of the EnKF workflow) and, over time, they converge. Variation in the spread of the model parameters allows for an evaluation of the statistical probability of a particular model state and an indication of the EnKF performance. For example, a divergence in the ensemble results indicated by a sudden expansion of the parameter space might suggest that the EnKF is having trouble fitting the observations ().

Five months before the 26 June 2018 eruption of Sierra Negra, we completed an EnKF analysis to track the volcano’s stress evolution that turned out to successfully forecast the timing of the subsequent eruption. Because the initial forecast was conducted as a test using a simple elastic rheology and did not assimilate InSAR data before 2014 or after January 2018, we also conducted retroactive forecasts (which we refer to as “hindcasts”) to evaluate our findings. The following section details the results of four EnKF numerical experiments, which are evaluated on their ability to track the stability of the Sierra Negra magma system through time (Fig. 2 and table S3): (i) the pre-eruption forecast that assimilated ground deformation observations from descending InSAR tracks into an elastic FEM model with a constant Young’s modulus, which was conducted before the 2018 eruption; (ii) a post-eruption “hindcast” that assimilated all pre-eruption descending InSAR observations into an elastic FEM model with a constant Young’s modulus, “nTd” (non–temperature-dependent); (iii) a post-eruption hindcast that assimilated descending InSAR observations into an elastic FEM model with a temperature-dependent Young’s modulus, “Td”; and (iv) a post-eruption hindcast that assimilated both ascending and descending InSAR observations into an elastic FEM model with a temperature-dependent Young’s modulus, “Tot.” The hindcasts with only descending InSAR data are included to provide direct comparisons with the forecast. The full hindcast, Tot, is the most complete evaluation of the InSAR time series data.

Fig. 2.

EnKF parameter estimation.

(A) The 3D FEM model setup for a pressurized reservoir. Full model details are provided in the Supplementary Materials. Four EnKF model suites were run to track the evolving Sierra Negra magma system. (B) The ellipsoidal pressure source in the FEM is free to rotate and dip in up-down (U), east-west (E), and north-south (N) space. The pre-eruption forecast model uses selected Sentinel-1 InSAR data from December 2014 to 26 January 2018 (red dots indicating the ensemble mean, with orange error bars indicating 2-sigma SD). Subsequently, three hindcasts were conducted using additional Sentinel-1 InSAR data up to the eruption: a non–temperature-dependent elastic model that assimilates descending InSAR data only (nTd; black dots with gray error bars, depicting 2-sigma SD), a model with a temperature-dependent Young’s modulus that assimilates descending InSAR data only (Td; blue dots with blue error bars), and a temperature-dependent model which assimilates both ascending and descending InSAR observations (Tot; green dots with green error bars). (C to E) Predicted spatial parameters, X and Y, and depth, Z. (F and H) Geometrical constraints, R1 and R2. (G) The evolution of magma reservoir pressure, dP. A pressure evolution was produced for the forecast on the basis of the trajectory of the previous 2 years of deformation (red solid line). (I) The dip of the reservoir, θ. (J) The strike of the longer R1 axis of the ellipsoidal reservoir, ϕ.

Pre-eruption forecast

In the fall of 2017, Sierra Negra was chosen as a target for testing near real-time data assimilation using the EnKF due to its extensive, ongoing deformation signal observed by the Sentinel-1 satellite. Although GNSS data are also available for most of the pre-eruption deformation cycle (), our efforts focus on InSAR data assimilation because the ability to evaluate volcanic activity using satellite data in remote locations where ground-based observations are unavailable is critical for assessing hazards at many volcanoes (). In addition, InSAR data provide good spatial constraints reducing nonuniqueness in the model fit, and the lower temporal resolution is less computationally expensive. An important consideration when using data assimilation techniques, such as the EnKF, is that each time step involves significant processing time. Future efforts will incorporate GNSS into the deformation time series data analysis as the EnKF technique is further developed and computational efficiency is improved.

The pre-eruption forecast made in January 2018 used a time series of InSAR line-of-sight (LOS) displacement calculated from 69 descending acquisition observations (07 March 2015 to 26 January 2018) from the Sentinel-1 satellite. We did not include InSAR data available between 2005 and 2014 because this initial experiment with EnKF data assimilation was tailored to using Sentinel-1 InSAR data (only available after December 2014). Of particular interest was how the deformation source geometry varied through time because rapid changes in source geometry and magma input may affect its stability and potential for eruption. However, after an initial spin up period (~10 time steps) as the ensemble stabilized, the geometry and location of the model reservoir converged and remained relatively constant throughout the remainder of the data assimilation (Fig. 2, C to J). Pressure increase was estimated along the boundary of the reservoir (Fig. 2G). Note that, because the EnKF begins with the first InSAR observation in 2015, the full magnitude of pressure accumulated after the 2005 eruption was not tracked, but rather the change in pressure between 2015 and 2018 was evaluated.

After incorporating the 26 January 2018 InSAR observations, the conditions for failure around the magma reservoir were evaluated to determine the stability of the system (Fig. 3). Andersonian fault orientations were calculated from the modeled stress state in the regions of predicted Mohr-Coulomb failure [using a value of C = 10 MPa for cohesion, following previous analyses of the 2005 eruption ()] to investigate potential faulting (, ). In the mean model from the EnKF ensemble, some elements near the southern edge of the magma reservoir exhibited tensile stresses >1 MPa, areas of Mohr-Coulomb failure were observed in the overlying roof, and pressure change in the magma reservoir was ~8 MPa (Figs. 2G and 3A).

Fig. 3.

Sierra Negra pre-eruption failure forecast.

North-south cross sections are shown through the center of the deformation source calculated by the mean EnKF model (star in Fig. 5B). The mean EnKF model calculated on 26 January 2018 is propagated forward in time following the pressure trajectory determined by the previous 2 years of deformation (Fig. 2G) to provide calculations of the stress state on 1 June and 26 June 2018. Colors represent calculated Andersonian fault orientations () in regions of predicted Mohr-Coulomb failure, assuming a cohesion of C = 10 MPa. Black dashed outlines indicate regions of calculated tensile failure near the reservoir boundary, assuming a tensile strength of Tc = 1 MPa.

Fig. 5.

Coulomb static stress transfer is calculated due to the 26 June 2018 Mw 5.4 earthquake, strike = 248, dip = 70 to the north, rake = 90, and Young’s modulus = 50 GPa, assuming receiver faults with strike = 105, dip = 70, and rake = 90.

(A) Cross section along x-x′ indicated on (B) through the assumed earthquake source fault (white line) and the location of fissure 1 (“F1”). The dashed black ellipse outlines the location of the pre-eruption forecast pressure source. The EnKF hindcast of Mohr-Coulomb failure (C = 1 MPa, gray outline) and tensile failure (green region) are shown for 26 June 2018. Dotted horizontal line indicates the 2-km depth of the Coulomb stress change plotted in (B). (B) Map view of Coulomb static stress transfer calculation at 2-km depth. Star symbol indicates the center of the forecasted pressure source that extends to the dashed black outline. The white circle indicates the center of the hindcast source, with its full extent outlined by the white line.

The calculated tensile stresses, pressure, and shear failure on 26 January 2018, 5 months before the eruption, were clearly not significant enough to drive eruption. A year-long forecast was produced by propagating the mean parameter values (and rates of change) from the EnKF ensemble forward through time to evaluate the stability of the system if it were to stay on this same trajectory. Because the pressure evolution was prescribed, the forecast model estimated magma system failure by investigating Mohr-Coulomb failure in the host rock and tensile failure along the magma reservoir boundary. As the model forecast progressed, calculated stress and failure accumulate in the roof above the magma reservoir and became more widespread (Fig. 3, B and C). The forecast model produced in January 2018 indicated tensile failure at the magma reservoir and through-going Mohr-Coulomb failure (i.e., continuous shear failure from the surface of the model to the boundary of the pressure source) was likely to occur between 25 June 2018 and 5 July 2018 (Fig. 3C). The 2018 forecast for Sierra Negra was presented in March 2018 at the UNAVCO (University NAVSTAR Consortium) Science Workshop as a rolling 10-day forecast. The forecast tracked the evolution of failure through 2018 and flagged the period of 26 June to 5 July 2018 as a potential time period for magma system failure (leading to eruption) due to through-going Mohr-Coulomb failure ().

Several caveats were discussed at the UNAVCO 2018 workshop including the lack of a temperature-dependent rheology in the forecast, the assumption of the failure criteria and host rock strength (e.g., cohesion and tensile strength), the lack of the consideration of the full stress evolution following the 2005 eruption (i.e., what was the initial stress state in early 2015?), and the assumption of the magma system maintaining its January 2018 trajectory going forward. Because our initial forecast for Sierra Negra was conducted as a test for the presentation at the March 2018 UNAVCO Science Workshop, we did not update the forecast in the months leading up to the 26 June eruption. In addition, we had not yet fully tested the temperature-dependent rheology FEM in the EnKF, so that is why it was not included in the initial forecast. Last, we had only used preliminary InSAR data in the forecast from descending observations for computational expediency. Given these many caveats, additional numerical experiments are necessary to evaluate the successful outcome of the forecast.

Post-eruption hindcasts

The Sierra Negra EnKF hindcasts (retroactive forecasts) use a time series of InSAR LOS displacement calculated from 98 descending acquisitions (Track 128, 13 December 2014 to 19 June 2018) and 42 ascending acquisitions (Track 106, 19 November 2016 to 18 June 2018) from the Sentinel-1 satellite (fig. S5). At each time step between InSAR scenes, the 240-member ensemble of 3D FEM models was calculated by COMSOL Multiphysics. As previous modeling efforts indicate that a temperature-dependent rheology may be necessary for Sierra Negra (), two of the hindcasts included an elastic rheology with a temperature-dependent Young’s modulus. One temperature-dependent hindcast used only descending InSAR observations for direct comparison with the pre-eruption forecast, while the other assimilated both ascending and descending observations into the temperature-dependent model (Tot). A non–temperature-dependent hindcast was provided for comparison. Since previous studies indicate that the viscous component of the rheology may be negligible for the time scale of the evaluated unrest period (), the EnKF has been developed to work with elastic rather than viscoelastic constitutive models in the finite element analysis.

There is generally very close agreement between the forecast and hindcasts as to the location and geometry of the pressure source (Fig. 2). Small differences are likely due to the improved InSAR time series data produced for the hindcast and improvements to the EnKF method. The largest differences between the four experiments appear in the estimation of pressure change (Fig. 2G), which varies from ~9 MPa at the time of the eruption for the total hindcast (Tot, green dots) to 30 MPa of pressure change in the temperature-dependent hindcast (blue dots). The increasing model covariance in the Td model, as observed by the expanding 2-sigma error bars, is due to a test of the EnKF for the Td model in which the parameters are normalized such that overpressure and radius have the same magnitude. Since parameter scaling did not improve the EnKF performance, and decreased its performance, it was not used in subsequent experiments. The spread in the ensemble from the Tot hindcast, which uses both ascending and descending InSAR observations, remains low, indicating a higher EnKF confidence.

An advantage of the ensemble modeling approach is that the percentage of models in failure can be tracked for a statistical evaluation of the potential eruption triggering mechanisms (Fig. 4). We use the total hindcast (Tot) to evaluate system evolution in the lead up to the eruption since it assimilated the most complete InSAR dataset and resulted in the best EnKF performance. The percentage of ensemble members experiencing Mohr-Coulomb failure in the host rock (Fig. 4A) and tensile failure along the reservoir boundary (Fig. 4B) is tracked at each data assimilation time step. As the cohesion (C) and tensile strength (Tc) of the rock are uncertain, several values are evaluated for each. In late 2017, a rise in the percent of models in the EnKF ensemble experiencing Mohr-Coulomb failure, using a cohesion value of C = 20 MPa, coincides with an increase in the recorded seismicity in the Sierra Negra caldera (). By the end of 2017, >60% of the models in the EnKF ensemble experienced Mohr-Coulomb failure (C = 20 MPa) in the roof above the reservoir. However, during this same period, <40% of the models in the EnKF ensemble experienced tensile failure along the reservoir boundary (Tc ≤ 1 MPa; Fig. 4B), and the mean ensemble model calculates no tensile failure.

Fig. 4.

Failure calculations from the Tot Hindcast.

(A) The percentage of models in the EnKF ensemble experiencing Mohr-Coulomb failure. A seismicity increase was observed in October 2017 by the el Instituto Geofísico de la Escuela Politécnica Nacional (IGEPN) seismic array () coinciding with a rapid increase in EnKF models in shear failure using a cohesion value of C = 20 MPa as indicated by the vertical gray line. Six months before eruption initiation (C) is indicated by the black dotted line, and the timing of the eruption (D) is indicated by the red dashed line. (B) The percentage of models in the EnKF ensemble experiencing tensile failure along the reservoir boundary indicates that 80 to 90% of the models are experiencing reservoir tensile failure in the days before the eruption for tensile strength, Tc ≥ 5 MPa. (C) The calculated failure for the mean EnKF model on 26 December 2017. The cross section runs south-north through the center of the pressure source. Colors represent calculated Andersonian fault orientations () in regions of predicted Mohr-Coulomb failure, cohesion, C = 10 MPa. No tensile failure is predicted along the reservoir boundary, Tc = 1 MPa. (D) The calculated failure for the mean EnKF model on 26 June 2018. The black dashed outline on the southern wall of the reservoir indicates the region of calculated tensile failure along the reservoir boundary.

In the lead up to the 26 June eruption, a greater percentage of models exhibit reservoir tensile failure (Fig. 4B) as the estimated change in pressure increases to 10 MPa (Fig. 2G). In the time steps before the eruption, >80% of models indicate tensile failure focused along the southern, shallower edge of the magma reservoir (Tc = 5 MPa; Fig. 4D), opposite of where most of the fissures erupted along the northern rim and flanks of the caldera (Fig. 1B). Curiously, the reservoir beneath the northern region of the caldera does not appear to be in tensile failure at the final time step in any of the ensemble models (Tc = 1 MPa). Rather, the northern side of the magma system is calculated to have remained in compression leading up to the eruption. Therefore, an additional catalyst or system wide stress change was apparently required to promote tensile failure and dike initiation to the north.

Andersonian fault orientations are calculated from the modeled stress state in the regions of predicted Mohr-Coulomb failure (here using a cohesion value of C = 10 MPa to be consistent with the pre-eruption forecast) to investigate potential faulting sources (Fig. 4, C and D) (, ). By late 2017, significant regions of the shallow roof experience shear failure (Fig. 4C). Normal faulting is predicted directly above the pressure source, with reverse faulting calculated above the outer edges. As the model progressed through 2018, the area of failure became more extensive and, by 26 June, shear failure was calculated to be through-going in the entire southern region of the system (C = 10 MPa).

Ten hours before the 26 June eruption, a Mw 5.4 earthquake struck on the southern portion of the Sierra Negra caldera along a north-dipping reverse fault coincident with the region of extensive Mohr-Coulomb failure calculated by the EnKF hindcast (, ). The calculated moment tensor source of the earthquake is in agreement with model predicted fault orientations (Fig. 4D). On the basis of previous calculations indicating that the 2005 eruption of Sierra Negra may have been triggered by a Mw 5.5 earthquake in a similar location (), we investigated the stress change due to the 26 June earthquake using the U.S. Geological Survey (USGS) Coulomb 3.4 software (, ). Coulomb static stress change indicates that the region to the north of the earthquake may have experienced significant unclamping in response to the event (Fig. 5). Since the magma reservoir was near tensile failure but had not yet ruptured, the trapdoor faulting event almost certainly triggered eruption. In addition, compressional stress is estimated to increase directly above the reservoir, promoting dike deflection and fissure opening to the north of the caldera (Fig. 5A). The trapdoor earthquake appears to have buffered the model-predicted tensile failure along the south edge of the reservoir, allowing the coulomb static stress change to induced failure on the northern edge instead.

DISCUSSION

Seismic precursors and earthquake triggering

A key feature of the Sierra Negra magma system is the interplay between the caldera trapdoor fault system and the magma chamber (, , , ). In the total hindcast (Tot) of the Sierra Negra system, overpressure and tensile failure remain insignificant, while shear failure becomes widespread in the surrounding crust. Given the timing and spatial location of the eruption in context with the Mw 5.4 earthquake, and the similarity to the sequence in 2005, it is likely (almost certain) that the two are intrinsically linked. If the Mw 5.4 faulting event had not occurred, then the EnKF indicates that the magma system was tending toward increasing tensile stress along the northern and southern edges of the reservoir (Fig. 3, C and D) and an eruption triggered through dike initiation was increasingly likely. However, it appears that the timing was expedited by the trapdoor faulting event.

The initial stress state

A critical issue in volcano forecasting is determining the initial stress state of the system. Unfortunately, capturing a full eruption cycle from an initial ambient stress state through to eruption is a difficult prospect, given that few systems have long-term monitoring spanning multiple eruptions, and a system may not reach a fully relaxed stress state after an eruption before continued unrest. In the case of Sierra Negra, the model forecasts and hindcasts were initiated when Sentinel-1 InSAR data became available in 2014, thus neglecting the prior decade-long stress evolution, which followed its 2005 eruption. During this time period, upward of 4 m of uplift was recorded by GNSS (), which is absent in the Sentinel-1 InSAR only approach. Capturing a full eruption cycle with the complete magnitude of deformation is an important next step in the development and testing of the EnKF method.

Fortunately, in the absence of a full eruption cycle, the EnKF can be tuned to capture failure and stress change as the system evolves. The rheology of the host rock assumed in the model can be adjusted in the absence of information on the pre-existing stress state. In the case of a protracted unrest period, a weaker crust is necessary to inhibit failure and eruption of the system. In that sense, the weakened Young’s modulus provided by a temperature-dependent rheology was the key for the Sierra Negra hindcasts. In addition, the chosen failure criteria are equally important. Hence, the failure criteria used in this investigation likely reflect a minimum rather than the true failure envelope for the system. Were the entire deformation period tracked, a higher tensile strength and cohesion would have been necessary to match the EnKF forecasts with the observations of seismicity and the timing of eruption at Sierra Negra. Future rock deformation experiments are necessary to constrain parameters such as cohesion and tensile strength. Until then, tracking a variety of failure envelopes is required.

Pressure evolution in the lead up to the eruption

A critical or maximum overpressure is often cited as a means for triggering an eruption through the initiation and propagation of a dike (). Our approach does not preclude pressure as an eruption catalyst but rather postulates that pressure buildup in a magma system is the means for promoting tensile and/or shear failure, which, in turn, triggers magma migration that may lead to eruption. However, volume change without significant overpressurization will result in the same strain accumulation in the host rock, leading to tensile and/or shear failure. Unfortunately, it is difficult to differentiate between the two (volume versus pressure) because of inherent nonuniqueness of the modeling approach (). The key, however, is estimating the stress in the host rock to determine whether the areas surrounding a magma system are near to failure.

The pressure state of the magma system may be linked to the intensity of the subsequent eruption. Hence, a key observation for constraining the pre-eruption magma reservoirs pressure state might be the magnitude of the eruption. In the case of the 2018 eruption of Sierra Negra, pressure within the reservoir was significant enough to drive diking and produce multiple fissure openings to the north and west of the caldera (). The apparently modest change in pressure estimated by the EnKF (~10 MPa) may be an appropriate order of magnitude and adequate to drive the 2018 eruption. However, significantly more research must be done to quantify pressure variations within a magma system to better understand the role of magma pressure in triggering dike initiation and eruptions. Until pressure variations can be better constrained, models of stress change and failure in the surrounding host rock, which can be linked more directly to observations of seismicity and deformation, provide an important approach for investigating eruption triggering mechanisms.

Forecasting the 26 June 2018 eruption

In many ways, the forecast provided by the EnKF 5 months ahead of the 26 June earthquake and eruption () can be chalked up to “accidental good fortune.” The forecast was based on rough estimates of physical properties, required the system to remain on the inflation trajectory determined on 26 January 2018, and relied on an assumption for what constitutes system failure. Specifically, system failure was flagged when Mohr-Coulomb failure calculated in the host rock was through-going, from the surface of the model to the magma chamber boundary (), using a cohesion of 10 MPa that appeared to work well for the 2005 eruption of Sierra Negra () and a constant Young’s modulus. It is unclear whether the through-going failure flagged by the EnKF forecast was forecasting the potential of the 26 June earthquake or the eruption. We posit that the more important outcome is the success of the EnKF to quantify deformation, stress, and failure as indicators to track the evolution of the system. The apparently successful eruption forecast of Sierra Negra illustrates the potential for evaluating magma system stress evolution in real time using the EnKF approach. This framework has transformative implications for forecasting volcanic unrest with higher fidelity in the future, which is particularly important in densely populated areas near active volcanoes.

MATERIALS AND METHODS

FEM approach

Our numerical approach uses previously developed and benchmarked thermomechanical FEM models (). COMSOL Multiphysics is used to calculate the stress, strain, and temperature variations because of a pressurized magma chamber in a 3D linear elastic space (Fig. 2A). A free surface is assumed at the top of the model space, and roller boundary conditions are applied on the side and base of the model. The magma chamber is represented by a pressurized ellipsoid that is free to move in space and rotate in strike and dip (Fig. 2B). Model parameters and variables are provided in tables S1 and S2, respectively. The geometrical and spatial parameters describing the magma reservoir are varied by the EnKF analysis described below.

The mechanical behavior of the model is governed by the quasi-static conservation of momentumwhere σ is the Cauchy stress tensor and b is the body force density vector.

The presented COMSOL models use the COMSOL Floating Network License for cluster computing, the Heat Transfer Module, and Structure Mechanics Module. EnKF results are plotted using Python, and the COMSOL MATLAB LiveLink is used for model visualization.

Thermal model for temperature-dependent hindcasts (Td and Tot)

A new steady-state thermal structure is calculated for each model in the EnKF ensemble for each time step (fig. S1A). The steady-state thermal structure is solved numerically by COMSOL from the steady-state heat conduction equationwhere k is the thermal conductivity, T is temperature, and Q is the crustal volumetric heat production, assumed to be zero. A magma temperature of 1100°C is assumed along the reservoir boundary. A background geotherm of 30°C/km is assumed. Although volcanic systems are transient and unlikely to reach a steady-state thermal equilibrium, this provides an end-member starting point. A nTd, elastic hindcast is provided for comparison.

We are particularly interested in the impact of the thermal structure on the elastic properties of the host rock and the resultant model predictions. Hence, we have incorporated a temperature- and depth-dependent Young’s modulus ()where a far-field, depth-dependent Young’s modulus (E = − 6.9z2 − 1.3 × 106z + 5 × 1010 Pa) is assumed in the brittle region of the model space (), Tm is the magma reservoir temperature, and Tgeo is the geothermal gradient (fig. S1B). Equation 3 provides a smooth transition between the brittle and ductile Young’s modulus to minimize computational issues and mimic nature, which likely has a transition in material properties rather than a sharp boundary.

Failure estimation

Failure in the host rock surrounding a reservoir is critically important for determining the stability of the system and potential for eruption. We use a combination of two approaches to predict magma chamber stability. First, we investigate faulting and failure in the brittle portions of the model space using a Mohr-Coulomb failure criterionwhere τ is the shear stress at failure, C is cohesion, f is the internal friction coefficient, and σn is the mean stress normal to the failure plane (). Second, we investigate the evolution of tensile stresses, σts, which are defined as the least compressive stress along the magma chamber boundary. In application, as a magma system grows and inflates, the expansion results in flexure and uplift of the overlying roof, promoting faulting and brittle failure. The Andersonian fault orientations of the model elements in the predicted region of failure are tracked to evaluate the fault types predicted during magma system evolution (, , , ). Simultaneously, tensile stresses along the chamber boundary can result in mode I failure and dike initiation ().

Ensemble Kalman filter

Model-data fusion strategies are critical for producing model forecasts of complex system behavior. We have adapted the EnKF (), an ensemble-based MCMC sequential data assimilation approach to assimilate large geospatial data into multiphysics volcano FEMs (, ). The ensemble-based EnKF can be applied with FEMs and circumvents the linearity and computational issues inherent to other Kalman filtering approaches (). In addition, the EnKF is highly parallelizable. The workflow (fig. S2) has been adapted for HPC using a handshake between Python and COMSOL Multiphysics. The HPC EnKF approach is highly scalable, and individual FEM models are distributed across compute nodes for swift, simultaneous calculation at each data assimilation time step. In practice, every time a new InSAR observation becomes available from Sierra Negra, the new data are assimilated to provide parameter updates for the models in the EnKF ensemble. The updated models are then propagated forward in time to provide updated forecasts of the volcanic system state. As more data, D, become available, the model errors are reduced and the forecasts are refined. Measurements and models are combined in the EnKF analysis step to provide the analysis ensemble, Aawhere X is the ensemble covariance matrix, C is the measurement covariance matrix, and H is the model operator matrix (–).

For the Sierra Negra implementation, 240 COMSOL FEM models are calculated at each time step on a cutting-edge HPC system at the National Center for Supercomputing Applications. A COMSOL Cluster Sweep is used to distribute the parameter values estimated by the EnKF analysis step to produce 240 models, which are distributed across compute nodes and CPUs (central processing units). Because of the inherent overhead of the COMSOL software, 240 models provide an optimal speedup. In addition, previous testing has indicated that as few as 100 ensembles provide a sufficient convergence for the EnKF approach (). Future EnKF implantations with open source modeling approaches will allow for more flexibility.

Coulomb static stress transfer

To investigate the static stress change resulting from the 26 June 2018 Mw 5.4 earthquake, we use the USGS Coulomb 3.4 Coulomb static stress software (, ). The Coulomb stress change is defined aswhere Δτ is the change in shear stress (positive in the slip direction), μf is the apparent friction coefficient, and Δσn is the change in normal stress (positive indicates unclamping).

The fault plane solution for the 26 June 2018 was determined by seismic wave form and InSAR analysis (). The scalar moment is estimated as 1.73 × 1024 dyne/cm, with a strike of 248°, dip of 65°, and rake of 90°. The location is thought to be along the southern trapdoor fault evidenced by a sinuous ridge in the Sierra Negra caldera. The location is further corroborated by ground observations collected after the event (). We calculate the Coulomb stress change for receiver faults oriented with a strike of 105° and dip of 65°, the complementary orientation of faults on the opposite side of the Sierra Negra caldera.

InSAR data processing

To measure the surface deformation over the Sierra Negra caldera, using MintPy software (https://github.com/insarlab/MintPy), we apply the small baseline InSAR time series analysis approach () to the Sentinel-1 descending track 128 subswath 1 dataset from 13 December 2014 to 26 January 2018 (69 acquisitions) and from 13 December 2014 to 16 June 2018 (98 descending and 42 ascending acquisitions) for the hindcasts. We generated a network of interferograms with five sequential connections for each acquisition using the stack Sentinel processor () within Jet Propulsion Laboratory/CalTech’s InSAR Scientific Computing Environment (ISCE) software (). We multilook each interferogram by 15 and 5 looks in range and azimuth direction, respectively, filter using a Goldstein filter with a strength of 0.2. We remove the topographic phase component using Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) [SRTMGL1, ~30 m, 1 arc sec with void filled ()]. The interferograms are phase-unwrapped using the minimum cost flow method (). We correct the displacement time series for the stratified tropospheric delay using the ERA-Interim weather reanalysis dataset () and for the topographic residual (). Reliable pixels are selected using a temporal coherence threshold of 0.7 ().

It is computationally prohibitive to assimilate data from the entire InSAR database. Hence, a QuadTree algorithm based on root mean square error of the displacement values is applied to reduce the number of samples for each epoch of InSAR data from ~150,000 to ~500.