The forecast is carried out by the model integration under the control of the model operator that represents the physical process governing the system. The model state at the current instance is forecasted from the model state at previous instance according to the model operator as expressed in Equation (5).

f and a: denote the forecast and analysis, respectively, M: the model operator derived from Equations (1) and (2) that relates the state of the system at time, t_i to t_i+1 (see Appendices A,B), and q: the model error. Initial conditions of the model state variables are necessary to start the model.

At each time step, the error covariance P^f of the model forecast which is an N_x × N_x matrix can be obtained using:

The update stage consists of correction of the model forecasts by the observations in order to obtain a more precise estimation of the model state. The observations D, i.e., D=[d1d2…dNm] where N_m is the total number of observations, can be anything that are related to the true model state x^† by the observation operator H.

Note that the observation error ϵ allows the derivation of the error covariance of the observations R, i.e., R = 𝔼(ϵϵ^T).

Updating X^f is straightforward and can be performed by first computing the Kalman gain K using Equation (8). The Kalman gain simply represents the magnitude of which the incoming observation corrects the model estimates in Equation (5).

To formalize the update, we use Equation (9).

The result of the update is called the analysis X^a. It is the linear sum of the model forecast and the correction introduced by the observation given the difference between the observations and model forecasts. In general, the analysis must be at least statistically as accurate as any of the individual observation or the model forecast.

Lastly, the error covariance of the analysis P^a is calculated using Equation (10).

We perform a model simulation to produce synthetic data by setting each model parameter to a given value as defined in Table 1. In data assimilation, this step is referred to as the true run. These values are chosen such that they are consistent with the case of Grímsvötn volcano in Iceland. In particular, at Grímsvötn Volcano, the shallow storage zone has been well characterized by seismic and geodetic studies (Alfaro et al., 2007; Hreinsdóttir et al., 2014). Other geometrical parameters are chosen in consistent with the study of Reverso et al. (2014). The value taken for the shear modulus G however, represents an upper bound for the lower part of the crust (Auriac et al., 2014). This large value results in large magma overpressure values but it does not influence the interpretation of the assimilation technique.

To produce the synthetic overpressures, the analytical solution of Reverso et al. (2014) to the differential Equations (1) and (2) are used, i.e., x†=[ΔPs†ΔPd†]:

with A=ad 3γdas 3γs+ad 3γd[ΔPdt0†−ΔPst0†+(ρr−ρm)gHc−8γsQinμHcas 3πac 4(as 3γs+ad 3γd)]. The resulting overpressures are considered as the true values and are then used as input in Equations (3) and (4) in order to generate synthetic displacements dm†, i.e., dm†=[uRm† uzm†].

The observation is then generated by adding a white Gaussian noise, i.e., N(μu=0,σu2), to the synthetic displacements, where σ_u represents the instrument precision along the radial and vertical directions. We use the typical GNSS instrument errors that are 1 and 10 mm for u_R and u_z, respectively.

Figure 4 shows the observations (vertical and radial displacements) generated at various distances from the volcanic center (i.e., r = 1 to 4.9 km). Notice that for all the radial displacements and for those vertical displacements that are nearest to the volcano axis (i.e., r = 1 to 2.5 km), it is very difficult to discern the synthetic observation because the amplitude of the noise is very small relative to the signal. Note that in the results part (Section 5) we use a total of 80 observations (i.e., N_m = 80) to define the observation vector D, which comprises the vertical and radial displacements, uniformly (every 100 m) located at distance r = 1 to 4.9 km away from the volcano axis. The effect of the number of observations used as well as the frequency of incoming observations are discussed in Section 6.1.1.

The state variables are the overpressures ΔP_s and ΔP_d and we consider two parameters, the radius of the deep reservoir a_d and the basal magma inflow rate Q_in as uncertain. We focused on these particular parameters characterizing the deeper part of the magma plumbing system because they are the most difficult to recover using surface displacement data.

The time interval Δt is fixed to 2 days, which is much smaller than the true time constant τ (i.e., τ = 0.11 yr). The start of the assimilation begins just after an eruption when the reservoirs are refilling and terminates after 500 time-steps (i.e., t_f ~ 2.74 yr). The initial values of magma overpressures are set to zero (See Table 1), assuming that both reservoirs had been fully depressurized by the previous eruption (i.e., the magma pressures within the reservoirs equal the surrounding lithostatic one). The frequency of available observation f_obs is usually fixed to one; meaning that at each assimilation step (i.e., every 2 days) there is always an available observation. In the case where no observation (f_obs = 0) is available, the model forecast cannot be corrected by the observation. This particular case is hereafter called the free run and is presented for comparison.

Two cases are considered here: case (A) where we only forecast the state variables and case (B) where the state variables and the two uncertain parameters (a_d and Q_in) are both estimated. We used the classic EnKF where the observations are perturbed prior to assimilation (Burgers et al., 1998). Below we describe in detail the practical implementation of EnKF as summarized in Figure 2:

The first test case tracks only the evolution of the overpressures in the shallow and deep reservoirs, i.e., Equation (A1), given an initial condition where the distributions of the uncertain parameters are unbiased (Figure 3C) or biased (Figure 5C).

In both cases, the EnKF performed well in forecasting the shallow and deep overpressures toward the end of the assimilation as evidenced by the red line (i.e., EnKF) that closely follows the black broken line (i.e., true value) in Figures 3, 5. In fact for the shallow overpressure, even with a poor knowledge about the uncertain parameters, the EnKF was able to catch almost perfectly the true overpressure. Unsurprisingly, the overpressure forecast for the deep reservoir is more likely affected by poor parameter initialization because the displacement induced by the deeper chamber is smaller, such that its overpressure is indirectly constrained by the dynamical model and consequently more influenced by the two uncertain parameters.

In Figure 5B, the result of EnKF is found closer to the free model forecast (i.e., results in green, labeled as free run) at the start of the experiment when the prior information about the uncertain parameters is far from their true values (see inset 2 for a magnified version). The model error seems smaller at the beginning as compared to the observation error, hence, the contribution of the model forecast is found to dominate the process. As the experiment goes on, the model error increases tremendously (i.e., the mean overpressure from the free run deviates away from the true overpressure) but thanks to the small measurement error, the analysis was able to converge closer to the true state of the system. Although overpressure estimation may have been unsuccessful toward the end of the assimilation, still, the resulting error differences between the true values and the EnKF-estimates are very small (i.e., 0.69% and 4.25% for the shallow and deep reservoirs, respectively). Note that this error difference may increase if the assimilation window is extended. The failure of estimation is due to the observation operator that relates the model and the observation, since it incorporates an uncertain parameter a_d which is fixed to a bad value in this case.

In Table S1, we present the summary of the synthetic results after performing the experiments for the state estimation case. As expected, there is a significantly higher percentage of error when the model is freely propagated forward in time without the correction of the observations (free run).

We estimated the displacements using the overpressures derived from the assimilation run, i.e., applying Equations (3) and (4). Figures 4, 6 show 10 out of the 80 combined radial and vertical displacements used during the assimilation. Notice that in the case where the uncertain parameters are well constrained (i.e., Figure 3C), the displacements are correctly forecasted even at distances where the signal-to-noise ratio starts to weaken (i.e., vertical displacements at r ≥ 3 km, see Figure 4). However, in the case of poorly initialized parameters (i.e., Figure 5C), forecasting the displacement is not favorable especially along the vertical component (Figure 6) where the result worsens as one goes farther away from the volcano axis.

There are two possible ways to solve the issue in forecasting the displacement: (1) extend the state vector X incorporating the model-forecasted displacements, i.e., X=[xfdf], where d^f = Hx^f or (2) extend the state vector incorporating the uncertain parameters (i.e., Equation A4). We opted for the second strategy because it will not only improve the field observation estimates but will also allow us to properly estimate the overpressures and to be able to constrain and gain more knowledge about the uncertain parameters. Furthermore, in performing the first option, the computational cost can increase significantly when we increase the number of observations used during the assimilation, i.e., using InSAR data.

Parameter estimation is very challenging especially when the parameters have no direct link to field observations and if there are no means to compare them to the actual values. The characteristics of the deep magmatic reservoir for example is poorly known in volcanology since it is buried very deep (i.e., >10 km) into the Earth.

We followed the same initial conditions performed in the first case given that: (1) the uncertain parameters are “unbiased” or have truncated-normal distributions, centered on their true values (Figure 3C) and (2) the uncertain parameters are “biased” or have normal distributions but are not centered on their true values (Figure 7C). We used the augmented state vector in Equation (A4) to include the uncertain parameters in the forecasting.

Results show that the filter tracks almost perfectly the true behavoir of the overpressures given the two different types of parameter initialization as shown for example in Figure 7, where the prior distributions of a_d and Q_in are biased. In fact, the final values of the overpressures have as little as 0.0−0.01% and 0.04−0.06% error difference with respect to their true values for the shallow and deep overpressures, respectively.

In Figure 8, we plotted a_d and Q_in at each assimilation step given two different a priori assumptions for the uncertain parameters. Notice how the EnKF-estimated parameters converged well near their true values regardless of how they were initialized. In fact, at the beginning of the assimilation, the filter immediately recognizes its supposed trajectory hence decreases a_d and Q_in allowing convergence to their true values. Interestingly, we find Q_in more sensitive to the estimation than a_d as evidenced by the steep drop at the start especially when the prior parameter distribution is biased. To a greater extent, it fell beyond the true value but eventually recuperates and adjusts to its correct behavior. In Table S2 we summarized the results of synthetic case B, showing that there is a very good fit between the EnKF-estimated state variables and model parameters and their true values.

Figure 9 confirms our initial recommendation that augmenting the state vector to include the uncertain parameters in the estimation will correctly forecast the radial and vertical displacements.

When performing the comparison between EnKF and inversion, we want to make sure first that the model and the a priori information that we used are the same and are suitable for the two techniques. For example, if we start from a parameter distribution which does not cover the true values of the uncertain parameters (e.g., Figure 5C), the inversion will not work. On the other hand, if we start from a uniform distribution as we used to do in inversion, the EnKF estimation will not be the optimal solution since it requires a prior assumption that is Gaussian in nature. We then built the prior distributions for a_d and Q_in (Figure 14B) that agrees with the prerequisites of both the assimilation and inversion.

Gregg and Pettijohn (2016) have previously compared data assimilation with inversion using two different models–a viscoelastic assumption for EnKF and an elastic one for the inversion. In our case, to be fair and consistent, the same forward dynamical model (the two-magma-chamber model) is applied to the two techniques.

For the EnKF, we performed the state-parameter estimation strategy similar to that in synthetic case B. Whereas, we implemented a Bayesian-based estimation, i.e., Markov Chain Monte Carlo (MCMC), for the inversion. MCMC is most useful in cases where models are non-linear and expressing an analytical solution is not possible (Segall, 2013). Note that we used 80 synthetic observations (i.e., 40 radial and 40 vertical displacements) that are uniformly (every 100 m) located at a distance r = 1−4.9 km with frequency of available observation consistent with the time interval (i.e., Δt = 2 days; f_obs = 1) for the two techniques.

In performing MCMC, the posterior distribution of the uncertain parameters is sampled given their a priori distribution using the forward model and a proposed likelihood distribution. Models are selected based on the Metropolis-Hastings rules, which always accept models that fit better to observations than the previous iteration and randomly accept those that do not fit to avoid being trapped to a local minima (Segall, 2013). For example, at t_i, a set of a_d and Q_in values are drawn from their initial distributions, generating model forecasts using the two-chamber model. These model forecasts are then compared to the GNSS and/or InSAR data and are always accepted if the fit is better than the last sampling iteration, creating the so-called Markov-Chain. The sampling iteration can be performed thousand to million times in order to build a full posterior distribution. In our case, we performed 11,000 MCMC sampling iterations at each time-step and burned-in the first 10,000 so that in the end we have a similar ensemble size (e.g., 1,000) to that of the EnKF. Take note that unlike in EnKF which is sequential and only uses incoming observations to capture the temporal evolution of the overpressures, MCMC utilizes all the observations from t₀ up to the preferred time of observation t_i. We performed up to 400 time-steps with interval similar to that in EnKF (i.e., Δt = 2 days).

Figure 14A presents the resulting overpressures after performing EnKF and MCMC. The shallow overpressure is well-estimated by the two approaches, whereas the overpressures in the deep reservoir both started from a deviated value but eventually converged to their true values. The MCMC-estimates are also found smoother than the EnKF-forecasts especially in the deep reservoir. Since EnKF always assumes that the dynamical model is uncertain and needs to be corrected by incoming observations, it then tends to closely follow the behavior of the observations, hence we observe a noisier estimation with EnKF. It follows that MCMC cannot account for epistemic model errors or those uncertainties related to processes not included in the physical forward model (Segall, 2013).

If we observe how the uncertain parameters are estimated by the two techniques in Figure 14C, we will find that both are able to estimate accurately the true values, but MCMC is faster to converge. Take note however, that the uncertain parameters we used in this study are static parameters. In the case of an evolving model parameter, which could be the case for the basal magma inflow (i.e., Poland et al., 2012), the inversion may not be the optimal method to use for estimation.

The model that we used here is a highly simplified view on how a volcano plumbing system works based on idealized assumptions (e.g., elastic medium in a homogeneous half-space, incompressible magma). It can be implemented easily for real-time forecasting of effusive eruptions as contrary to the finite element model of Gregg and Pettijohn (2016). Although we only considered two uncertain model parameters (a_d and Q_in) in this study, since they are the most difficult to constrain using geodetic data through conventional approach, future work can be extended to explore other parameters such as the depth and the shape of the reservoirs and the strength of the hydraulic connection between them.

In addition to the simplicity and specificity of the model that we tested, we base the potentiality assessment of data assimilation on a synthetic case that is consistent with observations recorded at a specific volcano –Grímsvötn in Iceland. In fact, similar deformation behavior has been observed at other basaltic volcanoes: Kilauea and Mauna Loa in Hawaii (Lengliné et al., 2008, Westdahl Volcano (Lu et al., 2003), and Axial Seamount Volcano (Nooner and Chadwick, 2009) such that it can be accepted as generic for this type of frequently erupting volcanoes. However, except for Axial Seamount, the time constant derived at other places are larger than the one observed at Grímsvötn (Reverso et al., 2014). This parameter is expected to only influence the impact of the temporal frequency of available observations. In the case of Grímsvötn, it is indeed expected to be more restrictive regarding the importance of a high temporal resolution dataset. It was beyond the scope of this first paper to modify the parameters chosen for the synthetic case, but a systematic exploration of the dimensionless parameters will be required in further studies.

We emphasize that the focus of this work is to give a preliminary assessment of how EnKF can be utilized in eruption forecasting rather than validating the dynamical model that we considered. Unlike in the field of ocean-atmosphere science where models are more advanced and established, realistic and generic physics-based models of volcanoes are still in progress. However, one of the main interests in using data assimilation is that it takes into account the fact that models are not perfect and are mostly based on the simplification of the complex reality. This is represented by the model error, q, as depicted in Equation (5). Evensen (2003) have shown that the model error can be accounted in the EnKF scheme using the following expression:

where Qti=qti2¯ is the model covariance. Pti+1f represents the accumulated model errors from the beginning of the assimilation until the instant t_i+1. In order to overcome the difficulty in quantifying directly the model error at each time-step, the EnKF represents the model error by an ensemble of model state generated from a large number of perturbations of uncertain model parameters. The model error covariance P in the EnKF practice is an approximation of the real model error. In case of infinite ensemble members, P is considered to equal the real model error. For this reason, a large ensemble size is always required. Any dynamical model can actually be used in data assimilation as long as there is a link between the model and the observations and the model can be restarted at any instant. However, take note that models that are too far from the reality would result in large model errors that would be difficult or impossible to correct by the observations, especially when the condition of the observations is not good enough (i.e., in terms of quantity, distribution and accuracy). The use of more realistic physics-based models that could better interpret field observations such as those that could account for magma rheology and compressibility (e.g., Rivalta and Segall, 2008; Anderson and Segall, 2011, 2013; Segall, 2016; Got et al., 2017) are then highly encouraged. Data assimilation can also be extended to models representing other plumbing mechanism such as magma reservoirs recharged by dikes at depth (e.g., Karlstrom et al., 2009) or even those eruptions that are related to dike intrusions.

While parameter-estimation allows us to gain more knowledge about the plumbing system and the behavior of the volcano, in real-time crisis, one of the key variables to infer an impending effusive eruption is the overpressure. The EnKF strategy presented here uses a simple dynamical model that can be easily integrated with large amount of real-time geodetic data, allowing to quickly and accurately track the value of the overpressures both in short-term and long-term periods. Assuming a statistic distribution for the threshold magma overpressure leading to reservoir wall rupture, based on previous eruption for instance, the updated overpressure provided by EnKF can be used to estimate the timing of an impending eruption. Although the critical overpressure value is not always known, especially for volcanoes that do not erupt frequently, it depends on the rock strength that can be estimated and on the local stress field strongly influenced by the edifice geometry (Pinel and Jaupart, 2003).

Another important challenge of the assimilation approach is the availability of frequent data. Although for GNSS we can obtain daily observations, InSAR data are less frequent and are still dependent on the quality of interferograms produced. Also, in reality, the 3D-displacement field vector from InSAR are not always retrievable due to the satellite's acquisition geometry. Furthermore, while we only used deformation data alone, the observation vector can include gas emission and seismic data for a more deterministic approach in forecasting (i.e., seismic data in particular can be used to estimate the timing of an eruption), as long as they can be related to the dynamical model used.

Interestingly, when it comes to near-real time monitoring, it may be possible to use inversion and data assimilation jointly in order to accommodate vast amount of incoming data. MCMC allows faster forecasting of non-evolving model parameters whereas EnKF via state-estimation is easy to implement and only requires incoming observations. Meaning, the uncertain model parameters can be first constrained by MCMC and the overpressures will then be forecasted using EnKF via state-estimation strategy (e.g., synthetic case A).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer GMC and handling Editor declared their shared affiliation, and the handling Editor states that the process nevertheless met the standards of a fair and objective review.