Research ArticleGEOPHYSICS

Stress inversions to forecast magma pathways and eruptive vent location

See allHide authors and affiliations

Science Advances  31 Jul 2019:
Vol. 5, no. 7, eaau9784
DOI: 10.1126/sciadv.aau9784

Abstract

When a batch of magma reaches Earth’s surface, it forms a vent from which volcanic products are erupted. At many volcanoes, successive batches may open vents far away from previous ones, resulting in scattered, sometimes seemingly random spatial distributions. This exposes vast areas to volcanic hazards and makes forecasting difficult. Here, we show that magma pathways and thus future vent locations may be forecast by combining the physics of magma transport with a Monte Carlo inversion scheme for the volcano stress history. We validate our approach on a densely populated active volcanic field, Campi Flegrei (Italy), where we forecast future vents on an onshore semiannular belt located between 2.3 and 4.2 km from the caldera center. Our approach offers a mechanical explanation for the vent migration over time at Campi Flegrei and at many calderas worldwide and may be applicable to volcanoes of any type.

INTRODUCTION

A basic, poorly investigated problem in volcano hazard studies is that we do not know where the next eruptive vent might form. This problem affects, to some extent, all active volcanoes, as there is always at least some chance of eruptive fissures opening in unexpected locations distant from the volcano summit. At central volcanoes with a characteristic cone shape (e.g., Hawaii or Etna), many eruptions are expected to occur on, or close to, the volcano summit. Eruptions also punctuate rift zones branching from the volcano summit, showing that new fissures may open at low altitude along the rifts, endangering populated areas. Some types of volcanoes, however, do not show any cone-shaped edifice and lack a summit focusing eruptive activity. This is the case of calderas, kilometer-sized subcircular depressions resulting from the drainage of magma chambers and collapse of their roof. Calderas have fed some of the most catastrophic eruptions on Earth and are extremely hazardous (1). However, their eruptions are generally few and far apart; thus, hazard is often underestimated by the local population, which, at some calderas, approaches 1 million. As calderas are associated with large volcanic fields, with past eruptive vents scattered throughout, the problem of forecasting the location of future eruption is extremely challenging.

Volcano hazard models, including lava and pyroclastic flow or plume expansion and fallout models, have reached high levels of sophistication but remain poorly constrained because of the large uncertainties on where magma will breach Earth’s surface to create eruptive vents, especially at calderas (2). The need for probabilistic maps of future eruptive vents has typically been addressed empirically, on the basis of the surface distribution of past vents (35). This is generally more or less implicitly justified on the basis of two different underlying assumptions. On one hand, it is sometimes assumed that previous vents mark weaknesses or paths in the host rock that will guide ascending magma. This assumption is not supported by observations, as many volcanoes are punctuated by monogenetic (e.g., used by only one eruption) vents, surrounded by tens of other monogenetic vents. According to a second rationale, the observed vents’ patterns reflect an unknown physical controlling mechanism (5), thereby justifying a data-driven approach to create probabilistic maps. An animation of the eruptive history in the last 15 thousand years (ka) before present at Campi Flegrei caldera (http://hazard-mapping.org/Campi-Flegrei.html) illustrates that after the initial phase, where vents are scattered, most vents open relatively close to clusters of previous vents. This qualitatively shows how data-driven approaches work better and better as the density distribution gets populated. However, the number of observed vents at volcanoes rarely exceeds a few tens; hence, the underlying physical mechanism remains sampled by very few eruptions to provide a sharp representation of the vent distribution probability, resulting in very coarse maps. The animation shows that occasional eruptions (including the last 1538 Monte Nuovo eruption) hit locations distant from previous vents: Forecasts for such low-probability locations will be affected by large uncertainties, as interpolation or extrapolation is needed where data are scarce or lacking. Such eruptions would have been hard to anticipate. Some studies have attempted to complement the limited statistics by considering fractures (6, 7), seismicity, tomographic images, geochemistry, and gravity (8, 9). This is also problematic, as there is no clear evidence that including such additional information improves vent forecasts. Fumaroles and fractures may have been caused by, rather than having been a driver of, dike propagation, thus forcing interpretation to a wrong conclusion. These approaches may be difficult to validate or falsify retrospectively.

Forward validation is an outstanding issue for vent forecasts, as these maps are not routinely subject to testing. Performance testing generally involves partitioning the data into calibration and validation datasets. In the case of probabilistic vent opening maps, this translates into recalculating the maps on the basis of the earliest eruptions and verifying the performance on the most recent eruptions. Calculating maps for an arbitrary point in time may not be straightforward if the calculations rely on poorly dated datasets; some of these predictive models may be thus inherently hard to test. Moreover, partitioning the data would lower the number of points used to set up the map, exacerbating the issues linked to the need of extrapolating or interpolating vent density. Despite these limitations, no better method has been proposed to address vent forecasting, hampering progress in volcano hazard modeling.

Here, we propose a new concept to vent forecasting based on the combination of physical models of magma transport with Monte Carlo statistics. Numerous theoretical and field studies have established that host rock stresses dictate magma pathways: Magma-filled fractures (dikes) feed eruptions, and during ascent to the surface, the dikes align themselves with the most energy-efficient orientation, which is roughly perpendicular to the least compressive principal stress axis, σ3 (1016). Contrary to intuition, preexisting faults appear of subordinate importance in guiding magma (1013), as their orientation in respect to the stress field is optimized for shearing movements; thus, opening along such planes is inefficient. In this framework, dike trajectories can be predicted, provided that we know, with sufficient resolution, the volcano stress field and the dike starting location (the magma chamber). Stress magnitudes and directions in crustal rocks, however, are generally very poorly constrained. Despite this, we show that magma trajectories, and thus eruptive vent locations, are so sensitive to stress variations that the previous vent locations at a volcano can be used to constrain the stress field to a sufficient degree of accuracy to render reliable physics-based vent forecasts possible.

RESULTS AND DISCUSSION

Concept of the vent forecast scheme

We propose a scheme based on three independent blocks: (i) a deterministic model for magma trajectories and thus vent opening where the dike trajectories intersect Earth’s surface; (ii) a deterministic model for the volcano stress field, needed to calculate the dike trajectories; and (iii) a probabilistic scheme to constrain, on the basis of the observed vents, the posterior distributions of all models’ parameters. In other words, under the assumption that crustal stresses govern dike trajectories, we will retrieve the model parameters’ distributions that lead magma to propagate from the known magma reservoir’s location to the location of the observed vents. Last, probabilistic forecasts are obtained by combining the deterministic models with the posterior parameter distributions. In the following, we first introduce the general scheme and outline some of the options available to deal with the three blocks described above. Next, we apply the scheme to vent forecasting at Campi Flegrei caldera. For this first application, we opt for the simplest deterministic and statistic models to enhance ease of explanation and model transparency.

State-of-the-art models for magma propagation trajectories are based on maximizing tensile stresses at the dike tip line (15) or energy release rate on prospective dike incremental elongations (14, 16). These models account for magma buoyancy pressure and allow for mixed-mode propagation and are thus capable of accurately treating the effects of layering or the free surface. These methods, however, involve several parameters. In our application to Campi Flegrei, we use a simpler option: Magma propagates strictly perpendicular to the least compressive principal stress axis σ3. According to this approximation, magma-driven fractures will be pure opening fractures (10, 11). In other words, development of mixed-mode fractures resulting, e.g., from the interaction of the dike with layering or the free-surface will be neglected; the resulting bias on the trajectories is discussed below. In three dimensions, the directions perpendicular to σ3 in every point of space identify a family of surfaces on which magma is assumed to propagate; these will intersect Earth’s surface as a curve on which the eruptive fissure will lie. In two dimensions, the surfaces become σ3-perpendicular streamlines, i.e., curves aligned in every point to σ1. These surfaces or curves are obtained from the eigenvalues and eigenvectors of the stress tensor in every point of the crustal volume under investigation.

Deterministic simulations of magma trajectories rely on a well-balanced stress model for the volcano. Layering affects edifice stresses and magma trajectories (16), but the benefit of considering it should be weighed against the need of including detailed information on numerous, often poorly constrained, structural parameters. For our application to Campi Flegrei, we take a simpler approach and assume a homogeneous half-space.

For the parameter inversions, efficient Bayesian resampling methods, such as Markov chain Monte Carlo (17), are probably the best option for the general case. In our application to Campi Flegrei, we resort to the computationally inefficient but simple Sampling/Importance Resampling (SIR) algorithm (18, 19). The SIR is a noniterative Bayesian Monte Carlo “perfect sampler” (17): A large number of magma propagation trajectories are calculated at the beginning of the procedure on the basis of prior distributions for the model parameters; from this set, a subset is selected (resampling) so that the modeled arrivals at Earth’s surface exactly reproduce the empirical distribution of observed vents. The resampling results in posterior distributions of the model parameters, which can then be used in forward mode to compute a large number of magma trajectories and thus well-populated probabilistic vent forecasts.

In the following, we describe the details of the general conceptual scheme and present a simple application to Campi Flegrei. Last, we discuss the potential and limits of the approach, considering the physical plausibility of the inversion results and testing the forecasts’ performance.

Parametric volcano stress model

We assume that the background state of stress (i.e., before the establishment of any volcanic edifice) is Andersonian (i.e., Earth’s surface is stress free and one of the principal stress axes is vertical). Under this assumption, given that the equations of elasticity are linear, trajectory calculations may rely just on the perturbations from a lithostatic stress state (15). We assume that the local state of stress will evolve at any time, t, and point in space, (x,y,z > 0), as the linear superposition of both slow and sudden processes affecting the volcano historyσTOT(x,y,z,t)=σT(x,y,z,t)+σC(x,y,z,t)+σL(x,y,z,t)+σU(x,y,z,t)+σI(x,y,z,t)+σE(x,y,z,t)(1)where σTOT(x,y,z,t) is the supralithostatic full stress tensor; σT(x,y,z,t) is the regional tectonic stress tensor; σC(x,y,z,t) is the stress perturbation due to magma chamber pressurization (20, 21); σL(x,y,z,t) is the stress perturbation due to edifice load (13, 14, 22); σU(x,y,z,t) is the stress perturbation due to unloading caused by mass redistributions owing to, e.g., icecap melting (23), excavation of a caldera (24), or a flank collapse (25, 26); σI(x,y,z,t) is the stress tensor due to previous magmatic intrusions; and σE(x,y,z,t) is the stress tensor due to previous large earthquakes or slow slip events. Multiple or additional terms may be appropriate in special cases. Coupling between these terms, where, for example, edifice loading influences magma chamber shape and thus its stress perturbation, are neglected in this first-order linear approach. Fully interacting approaches, relying on, e.g., boundary element calculations, are also compatible with our general scheme and can be considered in future studies.

The terms in Eq. 1 are not all equally important: They are roughly ordered according to a decreasing extent of the stressed rock volume, since these stresses decay away from the source over a distance that scales with the size of the pressurized area (13). This translates into the list being ordered roughly according to expected relevance so that σT, σL, and σU are often dominant. Stresses around a pressurized magma chamber control the location of magma chamber rupture and have been commonly considered the dominant component of edifice stresses (2022). However, there is accumulating evidence suggests that their importance in determining magma trajectories may have been overestimated. Recent volcano stress models considering a combination of tectonic and edifice loading and unloading stresses have successfully explained magma trajectories, dike orientation, and vent location in diverse settings (13, 2326), as opposed to corresponding ones involving only magma chamber stresses (20). Loading/unloading stresses scale with the height of the edifice or the caldera effective depth (i.e., including the contribution of low-density infill). Depending on the case, this may amount to large loading or unloading; e.g., a 4000-m-tall edifice or a 1000-m-deep caldera corresponds to ~100 and ~25 MPa loading or unloading, respectively. Moreover, loading/unloading decays approximately as R−1 over a vertical distance R that scales with the edifice or caldera radius. The distance over which magma chamber pressurization stresses decay depends on magma chamber shape, but in general, they are only intense in the proximity of the chamber and decay very fast away from it, e.g., as R−2 for a spherical chamber (13). We further stress the relative importance of the magma chamber and loading/unloading by calculating a simple finite element (FE) model comparing stress trajectories due to a shallow topographic depression of ~200 m with those around a pressurized sill and those with both depression and sill (Fig. 1, A to C). Consistent with previous studies comparing the effect of edifice loading and magma chamber stresses on magma trajectories (13), results show that, even for such a shallow caldera, the model that only contains unloading stresses has trajectories like those of the composite model, while the magma chamber stress model does not. Thus, the purely unloaded model is preferable, as it describes magma trajectories satisfactorily, with a lower number of parameters. Moreover, magma chambers depressurize upon magma transfer into a dike, so that magma chambers stresses further lose importance over the course of the pre-eruptive phase. In conclusion, tectonic stresses and stresses due to mass redistribution at the surface often dominate in volcano edifices except when very close to the magma chamber.

Fig. 1 Numerical Finite Elements (FE) models of principal stress orientation and resulting magma trajectories.

σ3 is represented by the black segments, the thick red line represents the location of a sill-shaped pressurized melt lens, and the thin red lines represent σ3-orthogonal streamlines. (A) “Unloading scenario” where we applied 9 MPa of vertically oriented tensional stress due to overburden removal (blue arrows). (B) “Inflating sill scenario” where we applied 5 MPa pressure to a thin flat cavity. (C) “Unloading + inflating sill” scenario where we combined the two previous cases. We added a weak (1 MPa) horizontal stretching to all models.

Similar to magma chamber pressurization stresses, σI and σE are very intense in the near field but decay rapidly with distance and, thus, are generally of smaller magnitude than σL and σU except when very close to the intrusion or fault (13). In other words, it takes many intrusions or large earthquakes to compensate for stresses induced by edifice growth or partial collapse. However, the effect of previous intrusions may still be important to consider: Subsequent intrusions may preferably arrange in a complementary fashion to cumulatively compensate for strains caused by, e.g., loading or flank dynamics. Unfortunately, old intrusions are often poorly constrained, and we neglect this effect in our application to Campi Flegrei.

Even if the stress model is well balanced, magma trajectories calculated deterministically on the basis of Eq. 1 may poorly match observations, as stresses are generally ill-constrained. During the lifetime of volcanoes, loading/unloading stresses may evolve in complex ways, with stress-releasing and homogenizing processes such as earthquakes or magma intrusions periodically alternating with collapse episodes and stress buildup processes, e.g., superposition of layers of erupted products. Thus, loading stresses calculated on the basis of the three-dimensional (3D) topography of the volcano, as if it were built instantaneously, are generally overestimated (20, 24, 26). We propose a probabilistic scheme to constrain poorly known stresses based on empirical data. To this aim, we rewrite Eq. 1 in parametric form asσTOT(x,y,z,t)=σT(t)+PC(t)GC(x,y,z)+PL(t)GL(x,y,z)+PU(t)GU(x,y,z)+PI(t)GI(x,y,z)+τE(t)GE(x,y,z)(2)where PC, PL, and PU are scalar chamber supralithostatic pressurization and gravitational loading and unloading pressures, respectively; PI is magma pressure minus the stress normal to the intrusions; τE is the static shear stress released by any substantial earthquake; and GC,LU,I,E(x,y,z,t) incorporate the variability in space and time of the stress perturbations induced by the individual processes, respectively. σT, PL, PU, PC, PI, and τE are the parameters to be inverted for. The assumption behind Eq. 2 is that the uncertainties on the stress field are represented by uncertainties on the scaling factor of each contribution, while their spatial variability is known to a much higher degree of certainty.

Location of dike nucleation

Eruptive vents are often interpreted as originating from vertical magma propagation from a “mirror” location directly below the vent. Thus, when the vent distribution shows migration patterns, this is often attributed to a migration of the magma source. Here, we argue that the curvature of stress-controlled magma pathways may offset significantly vent location from their source at depth (14, 24, 26), adding uncertainty to the starting location of the dikes. In our approach, the dike starting location can be either assumed on the basis of independent observations or treated as an additional unknown to be inverted for. We will show below that stresses and magma starting locations cannot be both determined independently with a high degree of accuracy, as they suffer from a strong trade-off. Thus, incorporating independent information on the magma starting location may be critical to an accurate resolution of edifice stresses. Crustal deformation may provide the most reliable constraint, as magma chamber shape strongly influences the location of magma chamber rupture (21, 22). At many calderas, crustal deformation is often consistent with thin horizontal penny-shaped cracks as sinks of melt accumulation. Such penny-shaped cracks are expected to rupture somewhere at their tip line, at a depth d and radial distance from the center r.

Probabilistic scheme and inversion procedure

First, the stress field due to the individual contributions in Eq. 2 is estimated from the 3D topography, edifice history, and structural information, e.g., by applying a distribution of loads on Earth’s surface to mimic topography (13, 14, 2426). Once all G terms in Eq. 2 have been estimated, prior distributions for the P terms, p(θ), are defined according to plausible ranges or any other prior information available. On the basis of those prior distributions, a set of random parameters are drawn, and the stress tensor (Eq. 2) is calculated and diagonalized to obtain the principal stresses. Trajectories are then calculated starting from the known location of the magma chamber, or, alternatively, depth d and radius r of dike nucleation can be included among the parameters inverted for. The dike arrivals at Earth’s surface for all combination of random parameters result in a modeled vent distribution p(x). A Bayesian resampling scheme will return the likelihood p(x|θ) and the posterior distributions of the model parameters p(θ|x).

In the SIR scheme, the distribution p(x) is resampled by only retaining in each bin xi a fraction of modeled ventsw(xi)=[Yo(xi)/Ym(xi)]/Σj[Yo(xj)/Ym(xj)](3)where w(xi) are the importance weights (18, 19) and Yo(xi) and Ym(xi) are the distribution of the observed and modeled vents, respectively. The resulting likelihood exactly replicates the distribution of the observations and restitutes posterior probabilities for the model parameters. Combining such posteriors with the deterministic stress and trajectory models results in a probabilistic stress field model that can be used to compute forecasts.

Explanatory models and forecasts

We now illustrate how our stress inversion procedure can help both improve the understanding of vent migration patterns during a volcano’s history and produce probabilistic forecasts that are testable at least in retrospective. We illustrate below this explanatory and predictive potential with application to Campi Flegrei. For our explanatory model, we invert for stress parameters and location of dike nucleation. For our forecasts, we will reduce the set of parameters by fixing the magma chamber depth and the radius to take advantage of information from inversion of crustal deformation data.

A precondition for a stress inversion is that tectonic and edifice stresses were roughly stable over the period of formation of a set of vents. Such a stationary stress state is an idealization of reality, as each batch of magma will change the stress field with its eruption deposits or from permanent host rock deformation. Nevertheless, as discussed above, a volcano eruptive history is often punctuated by a few dominant events. These events ideally partition the volcano’s eruptive history in a series of approximately stationary stress epochs that we can use to partition the vent data for the inversions. Inversion results from different epochs can then be compared with the edifice history to learn about how major edifice-modifying events are reflected in the estimated stresses. In particular, estimates of PL and PU obtained from the inversion, PL,m and PU,m, are expected to amount to a fraction of PL,o = ρLghL and PU,o = ρUghU, i.e., those expected from the observed 3D volcano structure, where ρL and ρU are the average density of the added crustal material and the excavated volume, respectively; g is the acceleration due to gravity; and hL and hU are the height of the edifice and the thickness of the overburden removal, respectively. The fractions RL = PL,m/PL,o and RU = PU,m/PU,o, where superscripts m and o indicate pressures resulting from the inversion and modeled on the basis of the stratigraphy, respectively, are a proxy for the longevity of elastic stresses in the volcanic edifice.

Application to Campi Flegrei caldera

We test our approach against the high-risk Campi Flegrei caldera (Italy). Campi Flegrei (Fig. 2A) caldera formed during the eruptions of the Campanian Ignimbrite ~39 ka ago and the Neapolitan Yellow Tuff (NYT) ~15 ka ago (27, 28). Postcaldera volcanism developed >70 monogenic vents focused in the northeast (NE), presently onshore, sector of the caldera. Eruptive activity migrated progressively inward over epochs 1 (15 to 9.5 ka ago), 2 (8.6 to 8.2 ka ago), and 3 (4.8 to 3.7 ka ago) (Fig. 2, A and B) (2932). Deposition of 25- to 300-m-thick, 3- to 25-m-thick, and 5- to 80-m-thick eruptive products over the three epochs, respectively, has partially refilled the inner caldera (29). Coeval resurgence uplifted the caldera central sector of ~180 m, of which >60 m ensued in the last ~5 ka ago (33). The last eruption occurred at Monte Nuovo, in 1538 (29, 34), whereas the most recent activity consists of four unrest episodes in 1950 to 1952, 1969 to 1972, 1982 to 1984, and 2005 to present, with an uplift of ~0.7, ~1.7, ~1.8, and ~0.4 m, respectively (35, 36). Inflation of a caldera-centered oblate spheroidal magma chamber at a depth of ~3.5 km is consistent with the deformation in the last ~600 years at least (3437). The inward migration of post-collapse volcanism (Fig. 2B) and the onshore focusing of vents (Fig. 2A) are both currently unexplained. A progressive shrinking of the magma chamber (38), invoked to explain the observed inward vent migration, is inconsistent with the approximately constant eruptive rates in the last 15 ka ago (Fig. 2C) (35) and the inferred size of the stationary shallow magmatic source in the last ~5 ka ago (34). Also, the onshore focusing of the vents, previously explained by the activity of a tilted resurgent block (39), is in contrast with later studies highlighting a nontilted resurgent dome (40).

Fig. 2 Eruptive history of the Campi Flegrei caldera during the last 15 ka.

(A) Shaded relief map of the caldera with location of eruptive vents (31). (B) Boxplot showing the radial distribution of eruptive vents. The boxes have horizontal lines at the first quartile q1, median, and third quartile q3. The whiskers correspond to 99.3% data coverage. The red circle at epoch 1 is an outlier [i.e., data point that is smaller than q1 − w × (q3 − q1), where w is the whisker length]. The horizontal size of the boxes represents the time interval of the epoch. (C) Volume of eruptions and cumulative (gray line) erupted volume (31).

Forward explanatory model for Campi Flegrei

Our working hypothesis is that the inward vent migration observed over the last 15 ka has been caused by stress variations. We drastically simplify Eq. 2 with the aim of keeping the number of parameters low. Loading from the volcanic edifice may be neglected because of the lack of an important topography [we explore the effect of the existing approximately SW (southwest)–NE topography gradient below]. We also neglect previous intrusions, earthquakes, and magma chamber stresses for the reasons explained above. Moreover, we assume axisymmetric geometry so that Eq. 2 will writeσTOT(r,z,t)=σT(t)+PU(t)GU(r,z)(4)where r is the distance from the caldera axis. The only stresses left in the equation are a homogeneous tectonic stress and a uniform unloading stress applied to the caldera floor because of the caldera excavation. We assume that σT(t) and PU(t) are piecewise constant functions with discontinuities at the time of transition between different epochs. We analytically calculate the function GU by means of equations describing the stress perturbation due to a rectangular negative (unloading) pressure source located at the surface. We fix the caldera radius (a = 7.5 km) and select as free parameters PU and the ratio σt/PU (the latter ratio for reasons explained later). In this section, we consider d and r as additional free parameters.

We run 107 Monte Carlo trajectory simulations based on random parameter sets drawn from uniform distributions: PU = −15 to 0 MPa, σT = 0 to 3 MPa, d = 0 to 5.5 km, and r = 0 to 6 km. The resampled distributions of PU, σT/PU, d, and r for the three epochs are bimodal (Fig. 3), indicating two distinct solution subsets according to whether σT or PU is dominant. If σT > 2PU/π (Fig. 3, B, G, and L), dikes ascend vertically (Fig. 3, E, J, and O, gray) from the edge of a magma reservoir that shrinks laterally over the three epochs (Fig. 3, D, I, and N, gray). We deem such reservoir shrinking scenario unlikely, as discussed above. Conversely, σT < 2PU/π corresponds to an “unloading-dominated scenario.” Here, dikes nucleate at shallow depth at r < 4 km (Fig. 3, C, D, H, I, M, and N, colored) and propagate to the surface in trajectories with upward concavity (Fig. 3, E, J, and O, colored). The likeliest value (mode) of PU decreases from −5.75 to −5.25 MPa from epochs 1 to 3 (the average PU varies from ~−7.9 to −7.7 MPa for epochs 1 to 3), equivalent to depositing 10 to 25 m of volcanic products over 200 to 300 m of overburden removal following the NYT eruption. The ratio between the inverted unloading pressure PU,m and that estimated from stratigraphy PU,o is RL ~ 0.8, without significant variations between the three epochs, suggesting that inelastic processes have partially relaxed and homogenized stresses (20, 26). The mode and average of σT/PU are 0.275 and ~0.26, so that σT = 1.5 to 2 MPa, consistent with evidence of weak extension from earthquake focal mechanisms (35) and fracture patterns (41). The wide spread around the likeliest values results largely from a trade-off between stress parameters and d (explored below), from simplifying the geometry to axisymmetric, from neglecting local stresses due to second-order topographic features, and from merging vents in discrete epochs instead of considering them individually. Notwithstanding such spread, the mean values of the stress probability density function resolve a small unloading stress change of ~0.2 to 0.5 MPa. This is consistent with the reloading expected from a thickness of ~10 to 25 m of eruptive deposits revealed by field and borehole data (29, 39) assuming a density of 2000 to 2500 kg m−3. Thus, the “unloading-dominated” solutions describe a shallow source of substantially stable size and depth, broadly consistent with geophysical and geological evidence (33, 35), overlaid by a progressively refilling caldera depression. The parameters appear affected by some trade-off (figs. S1 to S3). In particular, d is correlated with both PU and σT/PU, demonstrating that the two problems of magma source location and stresses cannot be easily separated on the basis of eruptive vent location. A trade-off between PU and σT is removed by considering PU and σT/PU instead (fig. S4, A to F). Moreover, d appears stable with respect to r, consistent with a melt accumulation zone over a relatively narrow depth range (figs. S1 to S3). Parameter trade-off appears limited for a narrow range (e.g., 0.5 km) of d (fig. S5), suggesting that a robust inversion is obtained if d is well constrained. We conclude that our model is capable of capturing the stress variations following the NYT eruption and the subsequent caldera infilling, only based on vent distribution.

Fig. 3 Forward model.

Summary of the resampled Monte Carlo simulations shown as histograms for epoch 1 (red), epoch 2 (green), and epoch 3 (blue). Gray and colored bars refer to simulations where σt is larger and smaller than 2PU/π, respectively, corresponding to “shrinking magma chamber” and “unloading controlled” scenarios. (A to D, F to I, and K to N) The resampled distributions for the input parameters as indicated. (E, J, and O) Representative streamlines for the two scenarios and the frequency of observed eruptive vents in the upper part of the panel.

Discussion on model assumptions

Sharp contrasts between layers of different rigidity and the interplay between propagation and the free surface affect dike trajectories (15, 16). Layering is common in volcanic areas, including Campi Flegrei (38). If rigidity decreases with decreasing depth (as commonly is the case), then trajectories become increasingly vertical, and vice versa. Ignoring such transitions may introduce a bias in the inversion results: Propagation into less rigid layers may be mapped into a larger extensional stress, and vice versa. The inversions will average out similar effects and return “effective” stress parameters. This is similar to the bias resulting from assuming the melt source deeper than it actually is: It will be mapped into biased stress parameters. If the inversion results are then used to make physical inferences on the stress field, then such bias should be considered and discussed. However, if the stress parameters are only used for predictive models, then the bias is not critical; in contrast, it is an element of strength of the proposed approach: Instead of complicating the model by accounting for the internal structure of the volcano, which is generally very poorly constrained, the simplest model possible will have the advantage of incorporating multiple effects into few effective parameters that will reflect a number of known and unknown processes into the forecast.

If the deep structure of the volcano was known in detail, then this could be incorporated into the method, but then simple principal stress trajectories would need to be abandoned for a dike propagation approach that can deal with crustal inhomogeneities by allowing for mixed-mode propagation, with the price of increased computing time and number of parameters. The question of whether accounting for layering and mixed-mode propagation results in an increased predictive power for the proposed approach is an interesting point for future studies.

As for the interplay between dike propagation and the free surface, numerical and analog studies have shown that the free surface induces mixed-mode propagation and increased curvature toward the free surface into a saucer-shaped intrusion (15, 42). Strong deviations occur when horizontal intrusions (sills) are very close to the free surface; i.e., their diameter is four to five times the thickness of the overburden (42), but trajectories start to bend slightly when their diameter is double their depth (15). σ3-perpendicular propagation pathways are pure opening pathways and thus cannot account for this complexity. However, this issue, while important, is not very critical in most cases, including Campi Flegrei, as the diameter of the sill-shaped magma chamber is similar to its depth. Cases of sill-shaped magma chambers with diameters larger than three times their depth can be addressed by using propagation models that can account for mixed-mode propagation (1416).

Forecast for Campi Flegrei

We now construct and test a predictive version of our approach. We reduce the number of free parameters by fixing the probability distribution for magma chamber depth to a generalized beta function (α = β = 2) ranging between 3 and 4 km, based on the inversion results in (38), and the nucleation radius to r = 0, thus inverting for only PU and σT/PU. Generalized beta functions are chosen here as they seem to provide good fits. Note that no loss of generality occurs owing to the r = 0 condition, because dike trajectories are horizontal for small radii; i.e., they lie on the assumed sill-shaped melt lens before bending toward the surface.

We use the data from epoch 1 to illustrate a stable-stress forecast, assuming that the condition of a stable state of stress is met. We use a random subset (two-thirds of the sample size, thus 20 of 31 vents) of the observed vent radii to resample PU and σT/PU. On the basis of the obtained posterior distributions, we produce a forecast by running 104 Monte Carlo simulations. Last, we test the forecast with the remaining one-third of the vents (Fig. 4, top). We find that this 2D forecast matches the distribution of radii well for those vents not used to set up the forecast. Note that having the magma chamber geometry fixed leads to stress parameters that are numerically different from those obtained with the explanatory model above; this is expected because of the trade-off highlighted above and does not invalidate either model.

Fig. 4 Vent forecasts.

Top: Stationary stress forecast. (A) Thin red lines are forecasted magma trajectories, and red stairs indicate the distribution of arrival radii based on 66% of the observed vents of epoch1. The blue and red bars show the distribution of observed vents and a random 66% of them, respectively. (B and C) Marginal distributions of inverted PU and σT/PU (red bars) fitted by beta functions (magenta lines), respectively. Bottom: Evolving vent forecasts. (D) Observed vent radii and resampled trajectories for epoch 1 [we excluded the outlier vent of Rione Terra (caldera center)]. (E and F) Resampled distributions of PU and σt/PU for epoch 1 fitted with a generalized beta function [magenta lines for (E), (F), (H), and (I)]. (G to I) Same as (D) to (F) but for epoch 2. (J and K) Distributions for PU and σt/PU projected for epoch 3. (L) Vent forecasting for epoch 3 (gray stairs) compared with observed vent radii (blue bars). (M to O) Same as (G) to (I) but for the Monte Nuovo cone.

We have demonstrated that, on the basis of only a few stress parameters (here two parameters), it is possible to forecast the location of future vents, provided that the stress field has remained approximately stable over the time interval used for the stress inversion. However, what if we want to forecast vents following a major stress-modifying event for the volcano, with no observed vents available to invert those new stresses and constrain such forecast? With this purpose in mind, we have developed a forecast approach to extrapolate the modeled stresses into the future. The approach takes advantage of the link between inverted stresses and stresses modeled according to the observed stratigraphy. We illustrate and test such evolving-stress forecast by using vent maps for epochs 1 and 2, along with stratigraphy information, to produce forecasts for epoch 3 and Monte Nuovo. We first resample the distributions of PU and σT/PU for epochs 1 and 2 using all available vents (Fig. 4, bottom, D to I). We then assume a linear relationship between inverted unloading and excavation stresses estimated from stratigraphy (and any other geological or geophysical information). This is equivalent to assuming that RU is approximately constant over the volcano history. We then assume that the average of the unloading pressure distribution PU,m and its support (minimum and maximum stress value) evolve as RUPU,O, which is proportional to the unloading stress due to the deposition observed during the successive epochs (Fig. 4, E, F, H, and I, and fig. S6).

PU0 is calculated as followsPU0=g(dbρb+Σ hiρi)(5)where g is the acceleration due to gravity; db = 2 km is the 40 ka Campanian Ignimbrite collapse depth; hi and ρi are the thickness and density of the post-Campanian Ignimbrite deposits (ρi ranges from 2400 to 2000 kg m−3), respectively; and ρb = 2600 kg m−3 is the density of the basement (table S1). Besides updating the average of the unloading distribution (fig. S6A), we also update the lower limit of the support of PU,m, Min(PU,m), by applying the same shift 〈ΔPU,m〉, on the assumption that a reloading of the caldera will affect both best guess and its lower limit. Since two possible updates of Min(PU,m) are available (starting from epochs 1 and 2), we take their average. The upper limit remains 0. The estimated 〈PU,m〉 for epoch 3 and Monte Nuovo are used to update σT/PU for the respective periods, assuming that the product 〈PU,m〉*〈σT/PU〉 remains constant (fig. S6B).

The forecast vent radii for epoch 3 have 56% of the observed vents falling within the first and third quartile (where 50% are expected) of the forecast distribution (Fig. 4L). For Monte Nuovo, the mode of the resampled distribution, corresponding to a 3.25-km radius, coincides with the observed distance of Monte Nuovo from the caldera center (Fig. 4O). While only one vent is insufficient to validate a forecast, we highlight that Monte Nuovo is located in a low-probability area from previous studies and would have been hardly anticipated if those maps had been available before the eruption.

Explaining nonaxisymmetric volcanism

Our axisymmetric approach returns a forecast for the radial distance, but not for the sector of the caldera (radial location) where the vent may open. This information may be approximately obtained from our 2D model by combining the stresses that we inverted for the end of epoch 3 (PU = −5.0 MPa and σT/PU = −0.4) with a first-order 2D stress perturbation created by the asymmetric topography of the caldera. We consider two end-member cross sections: a S-N offshore section with maximum unloading (negative relief) and a SW-NE onshore section including the highest relief outside the caldera (Camaldoli). We consider an average bathymetry of −100 m (offshore area) and an average topography of 300 m for Camaldoli (Fig. 5A). We use plane strain approximation and assume a homogeneous elastic half-space. The computed minimum compressive stress under the caldera floor is subvertical below 1.5 km and rotates progressively to subhorizontal toward the caldera rim (Fig. 5B), leading to concave trajectories. The larger surface loading toward the NE rim results in an asymmetric stress pattern, which drives most streamlines northeastward away from the magma lens. Eruptions occur mainly within this NE, onland part of the caldera (Fig. 5B). Only a few streamlines propagate southward and erupt outside the offshore caldera rim, consistent with the low number of observed vents (Fig. 2A). This suggests that the nonaxisymmetric distribution of the topography around the caldera may control the inhomogeneous distribution of its volcanism.

Fig. 5 Explaining asymmetric volcanism.

(A and B) Topographic/loading profiles (A, magenta and blue lines, respectively) and modeled stress field below Campi Flegrei caldera for the two half cross sections highlighted in the inset of (B). The black segments represent the direction of σ3. The thick red line represents the pressurized melt zone according to (38) and likely active in the last 5 ka at least (33). Thin red lines are σ3-perpendicular streamlines. Elev. a.s.l., elevation at sea level.

Our 2D results are broadly consistent with previous maps that indicate high probability of vents opening in the NE part of the caldera (6, 7). Yet, our results imply substantial vent opening density in locations where previous studies give a low probability, i.e., the first and third quartile of the forecast describe a predominantly onshore, approximately annular belt located between 2.3 and 4.2 km from the caldera center, on which the Monte Nuovo vent lies.

Forecast potential for other calderas

Our method may be applied to other calderas and volcanoes of different shapes. Below, we propose a general diagram forecasting the location of future vents for calderas based on their geometry (radius and depth), inferred depth to the magma reservoir, and stress conditions. This is intended to be purely demonstrative as important stress sources are ignored, including the presence of any stratocone, the density profile below the caldera floor, not to mention the high uncertainty related to tectonic stresses. We calculate how the chamber depth/caldera radius ratio, d/a, and the tectonic/unloading stress ratio, σT/PU, determine the arrival distance of the dikes feeding the vents; this arrival distance is classified into intracaldera, along rim (i.e., 0.85a to 1.15a), or off-caldera; these three domains correspond to the colored background in Fig. 6. We then locate on the diagram 17 approximately circular calderas (table S2) with a sufficient knowledge of the abovementioned parameters (caldera geometry, depth to the magma reservoir, vent location, and regional stress). PU was calculated considering only the maximum caldera depression and taking into account any water level. The regional stress was set more arbitrarily, assigning a stronger extension of 1 to 3 MPa to calderas in extensional settings and a weaker extension of 0.2 to 0.5 MPa induced by inflation or resurgence to calderas in neutral tectonics (1, 2). Despite its simplicity, our forecast broadly matches the observed dominant vent location (intracaldera, along rim, and off-caldera, shown as pie diagrams for each caldera). In particular, depending mostly on σT/PU, vents distribute from predominantly intracaldera to predominantly along the caldera rim and finally to off-caldera, confirming the importance of the stress ratio in determining their location. The considered calderas appear shifted to the left with regard to the theoretical domains, possibly depending on our limited and indirect knowledge of σT. This shift may also explain, at least partly, the two outliers from the overall pattern (Crater Lake and Santorini). These outliers may result from underestimating σT or magma buoyancy, or stress conditions different from those at the time of the past eruptions.

Fig. 6 Performance of the unloading model for notable worldwide calderas.

The background color represents the phase diagram for intracaldera, rim, or off-caldera vent location according to model prediction, as specified. The caldera rim field is bounded by 0.85a and 1.15a (inset). Estimates for calderas in nature are represented by pie plots according to the proportion of observed intracaldera, rim, or off-caldera vents, with the same color coding as the background. Uncertainties are shown as black lines. The star highlights the position of Campi Flegrei (CF) according to the σt/PU calculated in this study. The inset shows streamlines of selected models with d/a = 0.4 and d/a = 1.4 and variable σt/PU as a reference. d/a is depth d normalized to the caldera radius a, while σt/PU is the normalized stress controlling trajectory concavity. Caldera acronyms are as follows: FE, Fernandina; WO, Wolf; DA, Darwin; SiN, Sierra Negra; AL, Alcedo; CeA, Cerro Azul; AS, Aso; CrL, Crater Lake; RO, Rotorua; SA, Santorini; BO, Bolsena; RA, Rabaul; AI, Aira; VA, Valles; DO, Dolomieu.

CONCLUSIONS

In contrast to traditional approaches for vent forecasting, our method directly links observations, i.e., the distribution of observed vents, to properties of the volcanic edifice and the magmatic system, providing a physics-based forecast of future vent locations. The ability of our approach to resolve small stress variations developing during the eruptive history of a caldera represents a first validation of its potential. In the ideal case of detailed knowledge on the (re)distribution of surface loads during each individual eruption, and even stress change caused by any intruded magma, it would be possible to perform a stress inversion for each eruption. This could provide a quantitative approach to understand the evolution of stress over time and recalculate vent radius probability based on any future stress perturbation.

Unlike previous models, our approach performs well on retrospective tests: Our high-probability forecast for the location of Monte Nuovo was obtained after removing knowledge of its location from the model. Other geological and geophysical information on the structure and dynamics of a caldera may be included to better constrain the stress field in the volcano edifice and, in particular, the terms GL(x,y,z) and GU(x,y,z) from Eq. 2, which should improve the quality of the forecasts.

By improving the individual components of the present forecasting model (3D stress computation, including an accurate distribution of deposit coverage in time), our method may also be used to produce maps of probability of vent opening at calderas and in volcanic edifices of different shapes. Combining our method with hazard simulation tools may result in more reliable hazard assessment, especially for systems characterized by large vent location variability.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/7/eaau9784/DC1

Fig. S1. Probability distributions for model parameters of epoch 1.

Fig. S2. Probability distributions for model parameters of epoch 2.

Fig. S3. Probability distributions for model parameters of epoch 3.

Fig. S4. Covariance distributions of σt/PU versus PU and σt versus PU for a set of simulations with starting depth homogeneously distributed between 3 and 4 km and radius equal to 0 km.

Fig. S5. Impact of variability of starting depth.

Fig. S6. Parameters projection for the time-varying stress forecast.

Table S1. Thicknesses and densities of deposits filling Campi Flegrei caldera with associated loads.

Table S2. Properties of notable worldwide calderas.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank M. Bagnardi and D. Bindi for discussion. Funding: E.R. and F.C. were funded by the European Union Supersite MED-SUV project, grant agreement no. 308665. L.P. was funded by the German Federal Foreign Office through the German Humanitarian Assistance program, grant S05-41-321.50 IDN 03/16. V.A. was funded by the DPC-INGV project V2, “Eruptive Precursors.” T.D. was funded by the DFG-ICDP, grant agreement no. RI 2782/3-1. M.A.D.V. was funded by the DPC-INGV projects V1 and V2. Author contributions: E.R. conceived the study and the conceptual model. E.R. and V.A. coordinated the work and wrote the manuscript. F.C. and E.R. carried out the numerical simulations. L.P., F.C., E.R., and T.D. developed the statistical concept. M.A.D.V. provided geological and stratigraphic data. All authors have read and revised the manuscript and contributed ideas to the research. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Navigate This Article