Search

Understanding the stress evolution of extinct volcanoes can improve efforts to forecast flank eruptions on active systems. Field, petrographic, and seismic data are combined with numerical modeling to investigate the paleo-stress field of New Zealand's Akaroa Volcano, or Akaroa Volcanic Complex. Field mapping identifies 86 radially oriented dikes and seven lava domes found only within a narrow elevation range along Akaroa's erosional crater rim. These observations suggest that crater rim dike emplacement resulted from lateral deflection of vertically ascending intrusions from a centralized magma source, which in turn may have facilitated formation of the lava domes, as well as two scoria cones. We postulate that dike deflection occurred along a stress barrier, as neither a compositional change nor structural boundary are present. We use a finite element model (FEM) simulating Akaroa to test how different factors may have influenced the system's stress state and dike geometry. Elastic, non-flexural ("roller") model configurations containing a large, oblate, and shallow magma chamber produce stress barriers most conducive to radial dike emplacement along Akaroa's crater rim. These configurations also simulate rapid edifice construction above a preexisting lithospheric "bulge." Conversely, simulating flexural stresses exerted on the lithosphere by Akaroa's large mass hinder rather than promote radial dike emplacement. Temperature-dependent viscoelastic relaxation promotes gradual increases in stress barrier elevation, though this effect is strongly dependent on magma chamber parameters. These results suggest that Akaroa was constructed rapidly (within ∼100 kyr) prior to crater rim dike emplacement, which occurred throughout the volcano's remaining active lifespan.

Near-axis seamounts provide a unique setting to investigate three-dimensional mantle processes associated with the formation of new oceanic crust and lithosphere. Here, we investigate the characteristics and evolution of the 8˚20’N Seamount Chain, a lineament of seamounts that extends ~ 175 km west of the East Pacific Rise (EPR) axis, just north of the fracture zone of the Siqueiros Transform Fault. Shipboard gravity, magnetic, and bathymetric data acquired in 2016 are utilized to constrain models of seamount emplacement and evolution. Geophysical observations indicate that these seamounts formed during four distinct episodes of volcanism coinciding with changes in regional plate motion that are also reflected in the development of intra-transform spreading centers (ITSCs) along the Siqueiros transform fault (Fornari et al. 1989; Pockalny et al. 1997). Although volcanism is divided into distinct segments, the magnetic data indicate continuous volcanic construction over long portions of the chain. Crustal thickness variations along the chain up to 0.75 km increase eastward, inferred from gravity measurements, suggest that plate reorganization has considerably impacted melt distribution in the area surrounding the Siqueiros-EPR ridge transform intersection. This appears to have resulted in increased volcanism and the formation of the 8˚20’N Seamounts. These findings indicate that melting processes in the mantle and subsequently the formation of new oceanic crust and lithosphere are highly sensitive to tectonic stress changes in the vicinity of fast-spreading transform fault offsets.

In this study, synthetic numerical experiments are conducted to investigate how well the Ensemble Kalman Filter (EnKF) data assimilation approach distinguishes between two potential drivers of ground deformation at volcanic systems: pressurization and lateral reservoir expansion. Numerical models indicate that pressure-driven inflation creates larger radial displacements relative to inflation driven by lateral expansion. However, the introduction of noise can obscure these differences in simulated geodetic data. Although the EnKF does not fully reproduce the original synthetic models, it remains sensitive to changes in the magma reservoir's aspect ratio and is able to distinguish between the two inflation mechanisms. Ultimately, there remains significant non-uniqueness in how changes in reservoir pressure and size are reflected in surface deformation for any given aspect ratio, but future innovations may continue to improve filter performance.

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.

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

The accumulation of sizeable volumes of magma in the upper crust may produce plutons and/or result in supereruptions. Geophysical observations provide constraints on the rates, volumes, and melt distributions in magmatic systems, but they suffer from limited resolution and inherent nonuniqueness. Different, yet complementary, geophysical approaches must be combined with petrological, laboratory, and geochemical measurements. We summarize the results from such a combined approach from the central Andes. Taking a global perspective on large silicic systems reveals that several have >10% partial melt over large areas (10s of km2), and there may be localized zones with 50% or more.

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. Although this region has been investigated extensively, the ultimate controls on earthquake initiation and the duration of the seismicity remain unclear. In this study, we develop a finite element model for the Central United States to conduct a series of numerical experiments with the goal of determining the impact of heterogeneity in the upper crust, the lower crust, and the mantle on earthquake nucleation and rupture processes. Regional seismic tomography data (CITE) are utilized to infer the viscosity structure of the lithosphere which provide an important input to the numerical models. Results indicate that when differential stresses build in the Central United States, the stresses accumulating beneath the Reelfoot Rift in the NMSZ are highly concentrated, 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. By comparing the numerical results with three reference models, we argue that an extensive mantle low velocity zone beneath the NMSZ produces differential stress localization in the layers above. 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. These results show that, if boundary displacements are significant, the NMSZ is able to localize tectonic stresses, which may be released when faults close to failure are triggered by external processes such as melting of the Laurentide ice sheet or rapid river incision.

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.

Recent analytical investigations have suggested that magma buoyancy is critical for triggering catastrophic caldera forming eruptions. Through detailed assessment of these approaches, we illustrate how analytical models have been misapplied for investigating buoyancy and are, therefore, incorrect and inconclusive. Nevertheless, the hypothesis that buoyancy is the critical trigger for larger eruptions warrants further investigation. As such, we utilize viscoelastic finite element models that incorporate buoyancy to test overpressure evolution and mechanical failure in the roof due to the coalescence of large buoyant magma bodies for two model cases. In the first case, we mimic empirical approaches and include buoyancy as an explicit boundary condition. In the second set of models, buoyancy is calculated implicitly due to the density contrast between the magma in the reservoir and the host rock. Results from these numerical experiments indicate that buoyancy promotes only minimal overpressurization of large silicic magma reservoirs ( 100 km3) they surpass a rheological threshold where their subsequent evolution is controlled by host rock mechanics. Consequently, this results in a thermomechanical division between small systems that are triggered “internally” by magmatic processes and large systems that are triggered “externally” by faulting related to roof uplift or tectonism. Finally, critical assessment of recent analytical approaches illustrates that care must be used when applying previously derived analytical solutions to ensure that assumptions used in the original formulation are not violated during application to new geologic problems.

Large magma bodies that feed super-eruptions and build batholiths are not instantaneously emplaced. Many accumulate over time scales of 10 5 to 10 6 years as part of magmatic episodes that last 10 7 years and propagate a thermal and magmatic front through the crust to stabilize the reservoirs in the upper crust. This history imposes a rheological and thermodynamic conditioning on the host rocks that sets in motion three feedbacks that promote growth and longevity of large silicic magma reservoirs. Herein we review the development of ideas about the thermomechanical evolution of large silicic magma systems and explore the feedbacks and their implications for the growth, longevity, and evolution of large silicic magma reservoirs. Feedback 1 promotes increasing temperatures and consequent lower viscosities in the host rocks and the development of a ductile halo. Feedbacks 2 and 3 are feedbacks that result from the thermal dependence of the rheological properties of this ductile halo. In feedback 2 low wall rock viscosities lead to dissipation of strain in the host rocks reducing the likelihood of eruption. Feedback 3 is a negative loop between volume change and pressurization also reducing the likelihood of wall rock failure and eruption. We show that these feedbacks are most pronounced in larger reservoirs (> 500 km 3 ) and conspire to promote reservoir growth. Predicted imprints of these feedbacks are extended melt present lifetimes, complex heterogeneous age records and crystal-rich magma in some large silicic magma reservoirs. In this framework, interruption of the slow steady progress towards viscous death and solidification manifests as a supereruption. Second boiling and recharge (including buoyancy effects) acting in concert or independently lead to roof uplift and extension and eruptions are finally triggered by downward propagating faults from the extended and weakened roof. This connotes a thermomechanical division of calderas into those where eruptions are triggered “internally” by magmatic processes and those that are triggered “externally” by faulting related to roof uplift and attenuation. The division is controlled by size of magma reservoir, although, true to nature, exceptions exist, demonstrating interruption of the feedbacks by other processes like tectonism.

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.

Approximately every 100,000 years the Earth experiences catastrophic caldera-forming “supereruptions” that are considered to be one of the most hazardous natural events on Earth. Utilizing new temperature-dependent, viscoelastic numerical models that incorporate a Mohr-Coulomb failure criterion, we find that eruptive failure of the largest magma chambers is a function of the geometry of the overlying roof and the location of the brittle-ductile transition. In particular, the ductile halo created around the hot magma chamber buffers increasing overpressures and prevents pressure relief via magmatic injection from the magma chamber. The numerical results indicate that as chamber volume increases, the higher temperatures in the host rock and the decrease in the roof aspect ratio cause a shift from reservoir-triggered eruption to an external roof-triggered mechanism. Specifically, as overpressure increases within the largest magma chambers, extensive uplift in the overlying roof promotes the development of through-going faults that may trigger eruption and caldera collapse from above. We find that for magma chamber volumes > 103 km3, and roof aspect ratios (depth/width)

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.