Search

The submarine volcano Axial Seamount has exhibited an inflation predictable eruption cycle, which allowed for the successful forecast of its 2015 eruption. However, the exact triggering mechanism of its eruptions remains ambiguous. The inflation predictable eruption pattern suggests a magma reservoir pressure threshold at which eruptions occur, and as such, an overpressure eruption triggering mechanism. However, recent models of volcano unrest suggest that eruptions are triggered when conditions of critical stress are achieved in the host rock surrounding a magma reservoir. We test hypotheses of eruption triggering using 3-dimensional finite element models which track stress evolution and mechanical failure in the host rock surrounding the Axial magma reservoir. In addition, we provide an assessment of model sensitivity to various temperature and non-temperature-dependent rheologies and external tectonic stresses. In this way, we assess the contribution of these conditions to volcanic deformation, crustal stress evolution, and eruption forecasts. We conclude that model rheology significantly impacts the predicted timing of through-going failure and eruption. Models consistently predict eruption at a reservoir pressure threshold of 12–14 MPa regardless of assumed model rheology, lending support to the interpretation that eruptions at Axial Seamount are triggered by reservoir overpressurization.

Ensemble based data assimilation approaches, such as the Ensemble Kalman Filter (EnKF), have been widely and successfully implemented to combine observations with dynamic forecast models. In this study the EnKF is adapted to assimilate ground deformation observations from interferometric synthetic-aperture radar (InSAR) and GPS into thermomechanical finite element models (FEM) to evaluate volcanic unrest. Two eruption hindcasts are investigated: the 2008 eruption of Okmok volcano, Alaska and the 2018 eruption of Sierra Negra volcano, Galapagos, Ecuador. At Okmok, EnKF forecasts tensile failure and the lateral movement of the magma from a central pressure source in the lead up to its 2008 eruption indicating potential for diking. Alternatively, at Sierra Negra, the EnKF forecasts significant shear failure coincident with a Mw 5.4 earthquake that preceded the 2018 eruption. These successful hindcasts highlight the flexibility and potential of the volcano EnKF approach for near real time monitoring and hazard assessment at active volcanoes worldwide.

In recent years, the advent of ensemble-based methods in volcanology has greatly facilitated the use of numerical models within data assimilation frameworks that had previously been limited, either computationally or mathematically, to simpler analytical models. Because numerical models can simulate stress conditions throughout the model space, recent inversions based on assimilated volcanic deformation data are able to track not only the basic parameters of a magma reservoir, but also how those parameters affect the overall mechanical stability of the system. Although this approach has produced successful forecasts and hind-casts of volcanic eruptions, much work remains to be done in assessing its full capabilities and limitations. In particular, non-uniqueness in how source parameters are reflected in surface deformation can significantly impair the inversion’s ability to resolve the magma system’s true state and, by extension, the likelihood of eruption. While this problem is nearly intractable for deep reservoirs, for which changes in pressure and size are indistinguishable from deformation alone, preliminary synthetic tests at shallower systems have demonstrated a limited ability to resolve the main inflation mechanism. In this study, we investigate how the performance of an Ensemble Kalman Filter (EnKF) data assimilation framework varies under a wider range of experimental conditions than used in these initial investigations. In particular, we test how different mathematical implementations of the filter and how different levels of data availability affect the EnKF’s ability to distinguish inflation drivers and to accurately resolve reservoir parameters. To implement this experiment, two time series of synthetic GPS and InSAR data are generated, one in which deformation is driven by excess pressure and another in which it is driven by lateral expansion of the reservoir. For each filter implementation these datasets are down-sampled and given random noise prior to inversion, and after assimilation the resulting model is compared to the original synthetic conditions. We find that newer deterministic formulations of the EnKF are more accurate and consistent than the original stochastic implementation, although the improvement is relatively small. Moreover, some amount of parameter inflation is required to avoid model collapse, but more sophisticated adaptive inflation schemes do not produce better results than more basic formulations. Finally, we show that while increased data sampling does improve performance, this effect is subject to diminishing returns. In particular, data resolution near the center of inflation is more important than overall range of coverage. As new inversion techniques are developed or adapted from other fields, rigorous testing as demonstrated here will be a key step in being able to interpret future results and develop new forecasting frameworks for volcanic eruptions.

Ensemble based data assimilation approaches, such as the Ensemble Kalman Filter (EnKF), have been widely and successfully implemented to combine observations with dynamic models to investigate the evolution of a system’s state. Such inversions are powerful tools for providing forecasts as well as “hindcasting” events such as volcanic eruptions to investigate source parameters and triggering mechanisms. In this study, a high performance computing (HPC) adaptation of the EnKF is used to assimilate ground deformation observations from interferometric synthetic-aperture radar (InSAR) into high-fidelity, multiphysics finite element models to evaluate the prolonged unrest and June 26, 2018 eruption of Sierra Negra volcano, Galápagos. The stability of the Sierra Negra magma system is evaluated at each time step by estimating variations in reservoir overpressure, Mohr-Coulomb failure in the host rock, and tensile stress and failure along the reservoir boundary. The deformation of Sierra Negra is tracked over a decade, during which almost 5 meters of surface uplift has been recorded. The EnKF reveals that the evolution of the stress state in the host rock surrounding the Sierra Negra magma reservoir likely controlled the timing of the eruption. While increases in magma reservoir overpressure remained modest (< 10 MPa) throughout the data assimilation time period, significant Mohr-Coulomb failure is indicated in the lead up to the eruption coincident with increased seismicity along both trapdoor faults within Sierra Negra’s caldera and along the caldera’s ring faults. During the final stages of pre-eruptive unrest, the EnKF models indicate limited tensile failure, with no tensile failure along the northern portion of the magma system where the eruption commenced. Most strikingly, model calculations of significant through-going Mohr-Coulomb failure correspond in space and time with a Mw 5.4 earthquake recorded in the hours preceding the 2018 eruption. Subsequent stress modeling implicates the Mw 5.4 earthquake along the southern intra-caldera trapdoor fault as the potential catalyst for tensile failure and dike initiation along the reservoir to the north. In conclusion, the volcano EnKF approach successfully tracked the evolving stability of Sierra Negra, indicating great potential for future forecasting efforts.

Ground deformation observed using near-real time geodetic methods, such as InSAR and GPS, can provide critical information about the evolution of a magma chamber prior to volcanic eruption. Rapid advancement in numerical modeling capabilities has resulted in a number of finite element models targeted at better understanding the connection between surface uplift associated with magma chamber pressurization and the potential for volcanic eruption. Robust model-data fusion techniques are necessary to take full advantage of the numerical models and the volcano monitoring observations currently available. In this study, we develop a 3D data assimilation framework using the Ensemble Kalman Filter (EnKF) approach in order to combine geodetic observations of surface deformation with geodynamic models to investigate volcanic unrest. The EnKF sequential assimilation method utilizes disparate data sets as they become available to update geodynamic models of magma reservoir evolution. While the EnKF has been widely applied in hydrologic and climate modeling, the adaptation for volcano monitoring is in its initial stages. As such, our investigation focuses on conducting a series of sensitivity tests to optimize the EnKF for volcano applications and on developing specific strategies for assimilation of geodetic data. Our numerical experiments illustrate that the EnKF is able to adapt well to the spatial limitations posed by GPS data and the temporal limitations of InSAR, and that specific strategies can be adopted to enhance EnKF performance to improve model forecasts. Specifically, our numerical experiments indicate that: (1) incorporating additional iterations of the EnKF analysis step is more efficient than increasing the number of ensemble members; (2) the accuracy of the EnKF results are not affected by initial parameter assumptions; (3) GPS observations near the center of uplift improve the quality of model forecasts; (4) occasionally shifting continuous GPS stations to provide variability in the locations of observations results in better model predictions than utilizing fixed locations when the number of available instruments is limited; (5) spotty InSAR data coverage on the flanks of a volcano due to topographic shadows and/or atmospheric interference does not adversely impact the effectiveness of EnKF if the available coverage is > 50%; and (6) snow or glacial obstruction in the center of uplift can adversely impact EnKF forecasts. By utilizing these strategies, we conclude that the EnKF is an effective sequential model-data fusion technique for assimilating multiple geodetic observations to forecast volcanic activity at restless volcanoes.

Moving beyond the widely used kinematic models for the deformation sources, we present a new dynamic model to describe the process of injecting magma into an existing magma reservoir. To validate this model, we derive an analytical solution and compare its results to those calculated using the Finite Element Method. A Newtonian fluid characterized by its viscosity, density, and overpressure (relative to the lithostatic value) flows through a vertical conduit, intruding into a reservoir embedded in an elastic domain, leading to an increase in reservoir pressure and time‐dependent surface deformation. We apply our injection model to Interferometric Synthetic Aperture Radar (InSAR) data from the ongoing unrest episode at Laguna del Maule (Chile) volcanic field that started in 2007. Using a grid search optimization, we minimize the misfit to the InSAR displacement data and vary the three parameters governing the analytical solution: the characteristic timescale τ P for magma propagation, the maximum injection pressure, and the inflection time when the acceleration switches from positive to negative. For a spheroid with semimajor axis a = 6200 m, semiminor axis c = 100 m, located at a depth of 4.5 km in a purely elastic half‐space, the best fit to the InSAR displacement data occurs for τ P=9.5 years and an injection pressure rising up to 11.5 MPa for 2 years. The volume flow rate increased to 1.2 m3/s for 2 years and then decreased to 0.7 m3/s in 2014. In 7.3 years, at least 187 × 106 m3 of magma was injected., Key Points Our analytical and numerical model describes viscous magma propagation into a reservoirIncreasing conduit inlet pressure and volumetric flow rate accounts for the accelerating upliftAt least 187 million cubic meters of magma with viscosity 100 MPa s was injected between 2007 and 2014

We use numerical models to study the mechanical stability of magma reservoirs embedded in elastic host rock. We quantify the overpressure required to open tensile fractures (the failure overpressure), as a function of the depth and the size of the reservoir, the loading by the volcanic edifice, and the pore fluid pressure in the crust. We show that the pore fluid pressure is the most important parameter controlling the magnitude of the failure overpressure rather than the reservoir depth and the edifice load. Under lithostatic pore fluid pressure conditions, the failure overpressure is on the order of the rock tensile strength (a few tens of megapascals). Under zero pore fluid pressure conditions, the failure overpressure increases linearly with depth (a few hundreds of megapascals at 5 km depth). We use our models to forecast the failure displacement (the cumulative surface displacement just before an eruption) on volcanoes showing unrest: Sinabung and Agung (Indonesia) and Okmok and Westdahl (Aleutian). By comparison between our forecast and the observation, we provide valuable constraint on the pore fluid pressure conditions on the volcanic system. At Okmok, the occurrence of the 2008 eruption can be explained with a 1,000 m reservoir embedded in high pore fluid pressure, whereas the absence of eruption at Westdahl better suggests that the pore fluid pressure is much lower than lithostatic. Our finding suggests that the pore fluid pressure conditions around the reservoir may play an important role in the triggering of an eruption by encouraging or discouraging the failure of the reservoir.

Active volcanoes pose a constant risk to populations living in their vicinity. Significant effort has been spent to increase monitoring and data collection campaigns to mitigate potential volcano disasters. To utilize these datasets to their fullest extent, a new generation of model-data fusion techniques is required that combine multiple, disparate observations of volcanic activity with cutting-edge modeling techniques to provide efficient assessment of volcanic unrest. The purpose of this paper is to develop a data assimilation framework for volcano applications. Specifically, the Ensemble Kalman Filter (EnKF) is adapted to assimilate GPS and InSAR data into viscoelastic, time-forward, finite element models of an evolving magma system to provide model forecasts and error estimations. Since the goal of this investigation is to provide a methodological framework, our efforts are focused on theoretical development and synthetic tests to illustrate the effectiveness of the EnKF and its applicability in physical volcanology. The synthetic tests provide two critical results: (1) a proof of concept for using the EnKF for multi dataset assimilation in investigations of volcanic activity; and (2) the comparison of spatially limited, but temporally dense, GPS data with temporally limited InSAR observations for evaluating magma chamber dynamics during periods of volcanic unrest. Results indicate that the temporally dense information provided by GPS observations results in faster convergence and more accurate model predictions. However, most importantly, the synthetic tests illustrate that the EnKF is able to swiftly respond to data updates by changing the model forecast trajectory to match incoming observations. The synthetic results demonstrate a great potential for utilizing the EnKF model-data fusion method to assess volcanic unrest and provide model forecasts. The development of these new techniques provides: (1) a framework for future applications of rapid data assimilation and model development during volcanic crises; (2) a method for hind-casting to investigate previous volcanic eruptions, including potential eruption triggering mechanisms and precursors; and (3) an approach for optimizing survey designs for future data collection campaigns at active volcanic systems.

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. Although the EnKF has been proven effective to forecast volcanic deformation using synthetic InSAR and GPS data, until now, it has not been applied to assimilate data from an active volcanic system. In this investigation, the EnKF is used to provide a “hindcast” of the 2009 explosive eruption of Kerinci volcano, Indonesia. A two-sources analytical model is used to simulate the surface deformation of Kerinci volcano observed by InSAR time-series data and to predict the 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. Mohr-Coulomb failure models reveal that the host rock around the shallow magma reservoir is trending toward tensile failure prior to 2009, which may be the catalyst for the 2009 eruption. Our results illustrate that the EnKF shows significant promise for future applications to forecasting the eruption potential of restless volcanoes and hind-cast the triggering mechanisms of observed eruptions.

In this study, we utilize thermomechanical models to investigate how magma chambers overpressurize as the result of either magmatic recharge or volatile exsolution. By implementing an adaptive reservoir boundary condition we are able to track how overpressure dissipates as the magma chamber expands to accommodate internal volume changes. We find that the size of the reservoir greatly impacts the resultant magma chamber overpressure. In particular, overpressure estimates for small to moderate-sized reservoirs (1–10 km 3 ) are up to 70% lower than previous analytical predictions. We apply our models to Santorini volcano in Greece where recent seismic activity and ground deformation observations suggested the potential for eruption. The incorporation of an adaptive boundary condition reproduces Mogi flux estimates and suggests that the magma reservoir present at Santorini may be quite large. Furthermore, model results suggest that if the magma chamber is >100 km 3 , overpressures generated due to the high magma flux may not exceed the strength of the host rock, thus requiring an additional triggering mechanism for eruption. Although the adaptive boundary condition approach does not calculate the internal evolution of the magma reservoir, it represents a fundamental step forward from elastic Mogi models and fixed boundary solutions on which future investigations of the evolution of the magma can be built.

We present new analyses of volatile, major, and trace elements for a suite of glasses and melt inclusions from the 85°E segment of the ultra-slow spreading Gakkel Ridge. Samples from this segment include limu o pele and glass shards, proposed to result from CO2-driven explosive activity. The major element and volatile compositions of the melt inclusions are more variable and consistently more primitive than the glass data. CO2 contents in the melt inclusions extend to higher values (167–1596 ppm) than in the co-existing glasses (187–227 ppm), indicating that the melt inclusions were trapped at greater depths. These melt inclusions record the highest CO2 melt concentrations observed for a ridge environment. Based on a vapor saturation model, we estimate that the melt inclusions were trapped between seafloor depths (∼ 4 km) and ∼ 9 km below the seafloor. However, the glasses are all in equilibrium with their eruption depths, which is inconsistent with the rapid magma ascent rates expected for explosive activity. Melting conditions inferred from thermobarometry suggest relatively deep (25–40 km) and cold (1240°–1325 °C) melting conditions, consistent with a thermal structure calculated for the Gakkel Ridge. The water contents and trace element compositions of the melt inclusions and glasses are remarkably homogeneous; this is an unexpected result for ultra-slow spreading ridges, where magma mixing is generally thought to be less efficient based on the assumption that steady-state crustal magma chambers are absent in these environments. All melts can be described by a single liquid line of descent originating from a pooled melt composition that is consistent with the aggregate melt calculated from a geodynamic model for the Gakkel Ridge. These data suggest a model in which deep, low degree melts are efficiently pooled in the upper mantle (9–20 km depth), after which crystallization commences and continues during ascent and eruption. Based on our melting model and the assumption that CO2 is perfectly incompatible, we show that the highest CO2 concentrations of the melt inclusions (∼ 1600 ppm) are consistent with the calculated CO2 concentrations of primary undegassed melts. The highest measured CO2/Nb ratio (443) of Gakkel Ridge melt inclusions predicts a mantle CO2 content of 134 ppm and would result in a global ridge flux of 2.0 × 1012 mol CO2/yr.

High-resolution magnetic data collected along the axis and the south flank of the Puna Ridge, the submarine extension of the East Rift Zone of Kilauea Volcano, Hawaii, aid in the interpretation of magmatic processes. Crustal magnetization varies along the crest of the Puna Ridge. Three prominent magnetization highs are located between the depths of 1500 and 2500 m. Each of the magnetization highs is associated with distinct volcanic features: a flow field on top of fissured terrain, a lava plateau, and a staircase of lava terraces. The crustal magnetization highs can be explained by the presence of these volcanic features, taking into account water depth and the size of the volcanic features and making simple assumptions about the rock magnetization. No measurements of rock magnetization have been made on Puna Ridge samples. We assume two values: 50 and 10 A/m, based on values from mid-ocean ridge basalts and the estimated ages of Puna Ridge rocks [Clague et al., 1995]. The two shallowest magnetization highs are adjacent to a landslide scarp. Lava flows overtop the scarp suggesting that a dike has been emplaced and eruptions have occurred since the landslide formed. A possible scenario is that the collapse of the flank led to dike intrusion and eruption at this location because of a reduction in the horizontal stress within the edifice. Our observations suggest that flows are starting to fill in the gap caused by the landslide to rebuild this section of the north flank.

[1] We examine mantle melting, fractional crystallization, and melt extraction beneath fast slipping, segmented oceanic transform fault systems. Three-dimensional mantle flow and thermal structures are calculated using a temperature-dependent rheology that incorporates a viscoplastic approximation for brittle deformation in the lithosphere. Thermal solutions are combined with the near-fractional, polybaric melting model of Kinzler and Grove (1992a, 1992b, 1993) to determine extents of melting, the shape of the melting regime, and major element melt composition. We investigate the mantle source region of intratransform spreading centers (ITSCs) using the melt migration approach of Sparks and Parmentier (1991) for two end-member pooling models: (1) a wide pooling region that incorporates all of the melt focused to the ITSC and (2) a narrow pooling region that assumes melt will not migrate across a transform fault or fracture zone. Assuming wide melt pooling, our model predictions can explain both the systematic crustal thickness excesses observed at intermediate and fast slipping transform faults as well as the deeper and lower extents of melting observed in the vicinity of several transform systems. Applying these techniques to the Siqueiros transform on the East Pacific Rise we find that both the viscoplastic rheology and wide melt pooling are required to explain the observed variations in gravity inferred crustal thickness. Finally, we show that mantle potential temperature Tp = 1350°C and fractional crystallization at depths of 9–15.5 km fit the majority of the major element geochemical data from the Siqueiros transform fault system.

Gravity and bathymetry data from transform fault systems associated with mid-ocean ridges that are spreading at different rates indicates that crustal thickening may occur under transform faults associated with intermediate- and fast-spreading ridges. If confirmed, this observation would require a change in our understanding of how magma accretion occurs along rapidly spreading mid-ocean ridges. Mid-ocean ridge morphology and crustal accretion are known to depend on the spreading rate of the ridge. Slow-spreading mid-ocean-ridge segments exhibit significant crustal thinning towards transform and non-transform offsets1,2,3,4,5,6,7,8,9,10,11,12, which is thought to arise from a three-dimensional process of buoyant mantle upwelling and melt migration focused beneath the centres of ridge segments1,2,4,5,6,7,9,10,12. In contrast, fast-spreading mid-ocean ridges are characterized by smaller, segment-scale variations in crustal thickness, which reflect more uniform mantle upwelling beneath the ridge axis13,14,15. Here we present a systematic study of the residual mantle Bouguer gravity anomaly of 19 oceanic transform faults that reveals a strong correlation between gravity signature and spreading rate. Previous studies have shown that slow-slipping transform faults are marked by more positive gravity anomalies than their adjacent ridge segments1,2,4,6, but our analysis reveals that intermediate and fast-slipping transform faults exhibit more negative gravity anomalies than their adjacent ridge segments. This finding indicates that there is a mass deficit at intermediate- and fast-slipping transform faults, which could reflect increased rock porosity, serpentinization of mantle peridotite, and/or crustal thickening. The most negative anomalies correspond to topographic highs flanking the transform faults, rather than to transform troughs (where deformation is probably focused and porosity and alteration are expected to be greatest), indicating that crustal thickening could be an important contributor to the negative gravity anomalies observed. This finding in turn suggests that three-dimensional magma accretion may occur near intermediate- and fast-slipping transform faults.