Volcanic ash forecast using ensemble-based data assimilation: the Ensemble Transform Kalman Filter coupled with FALL3D-7.2 model (ETKF-FALL3D, version 1.0)

Abstract. Quantitative volcanic ash cloud forecasts are prone to uncertainties coming from the source term quantification (e.g., the eruption strength or vertical distribution of the emitted particles), with consequent implications for an operational ash impact assessment. We present an ensemble-based data assimilation and forecast system for volcanic ash dispersal and deposition aimed at reducing uncertainties related to eruption source parameters. The FALL3D atmospheric dispersal model is coupled with the ensemble transform Kalman filter (ETKF) data assimilation technique by combining ash mass loading observations with ash dispersal simulations in order to obtain a better joint estimation of the 3-D ash concentration and source parameters. The ETKF–FALL3D data assimilation system is evaluated by performing observing system simulation experiments (OSSEs) in which synthetic observations of fine ash mass loadings are assimilated. The evaluation of the ETKF–FALL3D system, considering reference states of steady and time-varying eruption source parameters, shows that the assimilation process gives both better estimations of ash concentration and time-dependent optimized values of eruption source parameters. The joint estimation of concentrations and source parameters leads to a better analysis and forecast of the 3-D ash concentrations. The results show the potential of the methodology to improve volcanic ash cloud forecasts in operational contexts.


emission rates and/or on model initialization. A large amount of scientific research has been conducted in recent years to: i) better quantify the amount of ash emitted, its vertical distribution across the column and the related uncertainties; ii) obtain data on the 3D structure of ash clouds, particularly using ground, aircraft, and space-based instrumentation; iii) improve model representation of physical processes occurring within ash plumes and clouds and; iv) transfer of scientific outcomes into operations. However, despite the substantial advances in model formulation and initialization, it is estimated that, in operational 5 contexts, forecasted ash concentrations still can have an uncertainty as large as one order of magnitude (e.g., IVATF, 2011).
Epistemic uncertainties in ash dispersal forecasts may have different origins, including: i) uncertainties in the source term (i.e. eruption column height, mass eruption rate, particle grain size distribution); ii) uncertainties in the atmospheric model driving dispersal simulations (e.g. wind velocity and direction, small scale turbulence intensity, atmospheric temperature and humidity) and; iii) uncertainties in model parameterizations of physical processes occurring both in the eruptive column and 10 during subsequent passive transport (e.g. ash settling and removal processes, particle aggregation, interaction with meteorological clouds, etc.). In addition to these, aleatoric uncertainties exist always regarding the future evolution of the Eruption Source Parameters (ESPs) when an eruption is on-going at the time of running a forecast. Several studies (e.g., Zehner, 2010;Kristiansen et al., 2012) have concluded that the main source of epistemic uncertainty in ash dispersal forecasts comes from the ESPs that, very often, are not well constrained in real time. 15 Inverse modeling and, in particular, data assimilation methods, are techniques that can be used to estimate the state of dynamical systems based on partial and noisy observations. In a broad sense, these techniques build on continuous or quasicontinuous observations to produce model initial conditions (analyses) that can be used to better predict the future state taking into account uncertainties in observations and model formulation. Data assimilation methods have been successfully applied to the estimation of the state of the ocean or the atmosphere (e.g., Kalnay, 2003;Carrassi et al., 2018) as well as for the optimiza-20 tion of uncertain model parameters (e.g., Ruiz et al., 2013). More recently, applications have been extended to atmospheric constituents (e.g., Bocquet et al., 2010;Hutchinson et al., 2017), including ash dispersion models with the purpose of estimating the 3D distribution of ash concentrations to be used as initial conditions for forecasts. Surprisingly, examples of application of data assimilation techniques to volcanic ash dispersion are scarce and still mainly limited to a research level. For example, Wilkins et al. (2015) implemented a data insertion methodology to improve the initial conditions of ash concentrations based 25 on satellite estimations of ash mass loadings in a Lagrangian dispersion model. Fu et al. (2015Fu et al. ( , 2017a

applied an Ensemble
Kalman Filter technique to the estimation of ash concentrations in an Eulerian dispersion model based on flight concentration measurements and satellite estimations using idealized experiments and real observations. Their results showed that both observational sets (flight measurements and satellite mass loads) reduced forecasts errors, in their particular case attributed to a wrong model representation of ash sedimentation processes. One important issue when using satellite estimates of ash mass 30 loadings is that observations only provide a 2-D distribution of ash mass while models usually require the vertical profile of ash concentrations. Fu et al. (2017b) presented a modified approach for the comparison between model and observations in the context of the ensemble Kalman filter that try to deal with this limitation.
Uncertainties in the source parameters can be circumvented in part by using inverse modeling techniques for the optimization of these parameters. Eckhardt et al. (2008) implemented a source parameter optimization approach based on the definition of a cost function which measures the departure of ash concentrations from observed values and the departure of the estimated parameters from their a-priori values. This allowed to reconstruct the full emission profile using data from different sensors. Stohl et al. (2011);Kristiansen et al. (2012);Denlinger et al. (2012); Pelley et al. (2015); Steensen et al. (2017) discussed further developments and evaluations of the proposed approach. In particular, Pelley et al. (2015) describe the operational implementation of this algorithm at the London VAAC. In Chai et al. (2017), the optimal parameters are found using a quasi-5 Newton minimization approach of a similar cost function and Lu et al. (2016) use a similar approach in the context of an Eulerian model. Finally, Zidikheri et al. (2017a, b) presented an optimization algorithm based on a systematic search of the optimal parameter values for both qualitative and quantitative ash forecasts and evaluate the performance of the technique for different cases showing a positive impact on forecast quality. Wang et al. (2017) performs idealized experiments in which a particle filter and an expectation maximization algorithm are used for the estimation of ash source parameters obtaining 10 promising results.
The goal of this paper is to contribute on the development of data assimilation methods to improve quantitative ash dispersion forecasts. To this end, we propose an ensemble-based data assimilation system for volcanic ash combining an Ensemble Transform Kalman Filter (ETKF) (Ott et al., 2004;Hunt et al., 2007) and the FALL3D ash dispersal model (Costa et al., 2006;Folch et al., 2009), named ETKF-FALL3D. This system produces a joint estimation of 3-D ash concentration and critical ESPs, 15 that can improve the performance of the classical ash dispersion forecast strategies. This manuscript presents a first analysis of the ETKF-FALL3D system using different Observing System Simulation Experiments (OSSEs) in which synthetic observations of ash column mass loadings are simulated and assimilated. The system is evaluated under constant and time-dependent ESPs and the sensitivity of the system performance to parameter uncertainty, ensemble size, and observations uncertainty is explored and discussed. Additionally, some impacts of the Gaussian assumptions underlaying the ensemble Kalman filter in the 20 present case are discussed. A description of the methodology is presented in section 2, the experimental setup of the sensitivity experiments is described in section 3, the results are discussed in section 4, and the final conclusions are outlined in section 5.

The FALL3D model
FALL3D is an Eulerian atmospheric dispersal model that solves the advection-diffusion-sedimentation equation for a set of 25 particle classes (bins), each characterized by a particle size, density and shape factor. Given an eruption source term and meteorological variables, FALL3D solves the 4D ash concentration for each particle class, from which the total and the fine ash column mass loadings are diagnosed performing a vertical integration. The meteorological fields must be furnished off-line by a Numerical Weather Prediction (NWP) model forecast or from a re-analysis dataset. The source term determines the amount of tephra injected to the atmosphere, its vertical distribution along the eruption column, and the fraction of mass associated to 30 each particle bin. This term can be parameterized using different schemes available in the model for the Mass Eruption Rate (MER) (e.g., Mastin et al., 2009;Degruyter and Bonadonna, 2012;Woodhouse et al., 2013) and the vertical mass distribution (e.g., Pfeiffer et al., 2005;Folch et al., 2016). For simplicity and without loss of generality, we will assume here a MER given by the Mastin et al. (2009) scheme, which depends on the fourth power of the eruption column height and does not account for wind effects, and a Suzuki vertical mass distribution (Pfeiffer et al., 2005) depending on two shape parameters A and λ. The Suzuki parameter A controls the relative height at which emission is maximum whereas the parameter λ controls the width of the distribution of mass around that level. A previous sensitivity test (Osores, 2018) has shown that the two FALL3D model parameters that affect most the model results are the eruption column height h and the parameter A in the Suzuki distribution.

5
For this reason, these two parameters will be used in the following sections to define the ETKF-FALL3D system experiments.

The ETKF-FALL3D System
In operational applications, data assimilation is implemented sequentially to provide an estimation of the state of the system at a series of times in the so-called "data assimilation cycle". Each data assimilation cycle consists on two steps: a first step in which the numerical model is used to provide an a priori estimation or forecast of the state of the system and its uncertainty, 10 followed by a second step in which the prior estimation is combined with observations (which are also considered uncertain), to obtain a posterior estimation or analysis. These two steps are repeated sequentially in order to propagate forward in time information from past observations. Let us assume that the state of a system at time t is represented by a state vector x t that, in our particular case, consists on the values of ash concentration at each model grid point and for each particle class. In other words, x t is a column vector with 15 n elements, being n the total number of state variables in the FALL3D model (i.e. the total number of grid points times the number of particle bins). For parameter estimation, model parameters θ, e.g. those defining the characteristics of the source term, are also considered part of the state of the system and thus assumed uncertain. For the sake of simplicity, we limit the FALL3D source term parameters to the eruption column height h and the Suzuki distribution A − Suzuki parameter, but the methodology that follows can easily be extended to any other set of model input parameters. The augmented state vector s t 20 at time t is defined as the concatenation of the state vector x t and the (time-dependent) estimated model parameters θ, that is, is a column vector with n s = n + 2 elements.
In the Ensemble Kalman filter the time-dependent uncertainty in the state variables and parameters is estimated using a Monte-Carlo approach through an ensemble of augmented states. Let us assume that we start at time t-1 with an ensemble of initial conditions and model parameters. Then, our forecast of the state of the system at time t, is obtained by integrating in 25 time the FALL3D model for each ensemble member: where M t represents the FALL3D model operator which integrates the model in time for the i-th ensemble member starting from the i-th initial conditions (analysis) x a t−1 and fixing the model parameters to θ a t−1 during the time integration interval. Note that a persistence model is assumed for the model parameters (i.e. θ f t = θ a t−1 ) since no information about its variations 30 is available yet during the forecast. Following the assumptions of the ensemble Kalman filter, the joint a priori probability distribution of the augmented state at time t is assumed Gaussian, with a mean and a covariance matrix estimated from the ensemble of forecasts: (2) where s f t is the ensemble forecast mean, P f t is the ensemble forecast covariance matrix (a square matrix of dimension n s × n s ), and S f t is the ensemble forecast perturbation matrix whose i-th column is computed as S t . Note that the integration of the ensemble in time propagates the uncertainty on the initial conditions and parameters at time t-1 into the future in order to provide a time dependent estimation of the forecast uncertainty. This is a key feature that makes these methods particularly appealing for the estimation of uncertain model parameters (e.g. Aksoy et al., 2006;Ruiz et al., 2013) and for an accurate quantification of concentration.

10
At the analysis step a set of observations is available which is related to the true state of the system by the following expression: where y t is a m-size column vector containing the value of the m observations at time t and x true is the true model state (assumed to be unknown). H is a (usually non-linear) transformation that maps the state variables (i.e. ash concentrations 15 for different particle sizes) into the observed quantities (e.g. cloud column mass load) and the vector represents the error in the observations. This error is typically assumed to be a zero-mean Gaussian random variable with covariance matrix R (dimensions of m x m). The errors in the observations are assumed to be uncorrelated in time and independent of the state of the system. Under these assumptions, the information provided by the forecast and the observations can be combined to obtain an estimation of the augmented state that minimizes the root mean square error with respect to the unknown truth (e.g., Kalnay, where s a t is the a posteriori estimation of the augmented state (i.e. the analysis) and H t is the tangent linear of the observation operator. The factor P f t H T t (H t P f t H T t + R) −1 is usually referred as the Kalman gain. The Kalman Filter equations also provide a way to estimate the uncertainty of the analysis. After the assimilation of the observations, the augmented state 25 covariance matrix is updated to: where P a t is the posterior or analysis augmented state covariance matrix. Note that (5) and (6) can be used to obtain an ensemble of analyses for the state variables and the parameters whose ensemble mean is equal to s f t and the perturbations are sampled from a Gaussian distribution with zero mean and covariance matrix equal to P a t . These equations can be difficult to solve explicitly for high dimensional systems due to the large size of P t and R t but several methods have been proposed to address this issue and to implement the ensemble Kalman Filter in high dimensional systems. In the present work, we use the 5 ETKF approach which solves the ensemble Kalman Filter equations in a sub-space defined by the ensemble members. Details about the equation that arises from this particular implementation can be found in Hunt et al. (2007), but a summary is given in Apendix A. One of the main advantages of this approach is that finding the analysis ensemble mean requires to invert a k x k matrix, which is significantly cheaper than inverting the n x n matrix for the case in which k « n (which is usually the case for high dimensional applications of the filter).

10
The process is schematically shown in Figure 1. The cycle starts with an estimation of the mean parameters, assuming they have a Gaussian distribution, k random samples are taken. Each parameter sample is used in one of the ensemble members integrated with the dispersion model. When an observation is available, it is combined with the ensemble forecasts using the ETKF equations. From this combination an ensemble of analysis is obtained with a set of optimized parameters that also has a Gaussian distribution. Then the next cycle starts from the ensemble of analysis and the set of optimized parameters, to produce 15 a new ensemble forecast. When a new observation is available, the assimilation method is applied, and the cycle continues so on.

ETKF-FALL3D experimental set-up
To explore the capability of the ETKF-FALL3D system we use an OSSE approach, in which a long model integration is performed and regarded as the true evolution of the ash cloud. This model integration will be referred as the nature run. 20 Observations are simulated from the nature run and then assimilated with the ETKF-FALL3D system. The June 2011 Puyehue-Cordón Caulle eruption (Osores et al., 2012;Collini et al., 2013) has been selected for the generation of the nature run. The Particle Total Grain size distribution (GSD) is represented by 12 classes with diameters between 2 mm (-1φ ) and 1µm (10φ) and densities ranging between 400 for the larger particles to 2100 kg m −3 for the smaller ones (Bonadonna et al., 2015).

Ash mass loading observation simulations
The vertical distribution of the source is parameterized using the Suzuki scheme considering λ=5, the settling velocity model is that of Ganser (Ganser, 1993), and the vertical and horizontal turbulent diffusion are parametrized by the similarity (Ulke, Forecasting System (GFS) analysis with an horizontal resolution of 0.5º, a temporal resolution of 6 hours and 27 constant pressure vertical levels.
The simulated observations represent ash mass column loads (vertically integrated ash mass per unit area) estimates retrieved from satellite radiances (e.g., Prata and Prata, 2012;Francis et al., 2012;Pavolonis et al., 2013). Simulation of satellite retrievals are chosen since these observations are available almost globally and have a high spatial and temporal resolution making them 5 particularly appealing for the implementation of operational data assimilation systems. To represent some of the limitations of current satellite-based ash mass load retrievals, the simulated observations are available only where the true load values are between 0.2 g m −2 and 10 g m −2 (e.g., Wen and Rose, 1994;Prata and Prata, 2012;Pavolonis et al., 2013). The observational error is assumed to have a random Gaussian distribution, with a standard deviation of 0.15 of the ash mass load.
For the sake of simplicity observations are assumed to be co-located with the model grid points, we also assume that 10 observation errors are uncorrelated (i.e. R is diagonal) and that observations are unbiased. All observations are generated assuming a clear sky condition both above and below the ash cloud. Two nature runs were generated to evaluate the ETKF-FALL3D system: one with constant emission profiles and another with time-varying emission profiles.

Constant emission profile
This nature run simulation considers a source term that remains constant during all the simulated period, with an eruption 15 column height of 8.5 km above the vent and a A − Suzuki parameter of 5.5 (Figure 3a). Figure 3b shows the ash mass loading from the nature run on 7 th June at 12:00 UTC for illustrative purposes. The number of available observations (which depends on the thresholds described in the previous section) is time dependent (ranging from 27 to 52 grid point observations) and, in this particular case, is primarily affected by the atmospheric circulation which produces variations in the 3D ash concentration within the model domain.

Variable emission profile
In this experiment, h and A − Suzuki are time dependent (Fig. 3a). In order to represent a realistic variability of the source parameters, the h evolution corresponds to the estimated heights during the 2011 Puyehue-Cordón Caulle eruption (Osores et al., 2014). Since the A − Suzuki parameter can not be directly estimated, the evolution of this parameter is simulated using an autoregressive model (Fig. 3). 25 In Fig. 3c, the 7 th June at 12:00 UTC ash mass loading field is shown. In this experiment the number of observations that are assimilated depends on the emission profile as well as the wind field, and it can range from 15 (on 11 th June 06:00 UTC) to 86 (on 11 th June at 18:00 UTC).

Data assimilation experiments setup
In the data assimilation experiments performed in this work, the simulated observations are assimilated every 6 hours. The 30 number of ensemble members in the experiments is set to 32 (unless stated otherwise). In most experiments source parameters are assumed to be unknown and estimated within the data assimilation cycle. The model grid, boundary conditions and all other model parameters and configuration options are as in the nature run. The ensemble at the first assimilation cycle is initialized using zero ash concentrations for all members and a set of parameters which are sampled randomly from a Gaussian distribution. The relaxation to prior spread inflation approach (RTPS, Whitaker and Hamill (2012)) with a parameter of α = 0.5 is applied to the state variables to reduce the impact of sampling error. For the parameters the conditional inflation approach of 5 Aksoy et al. (2006) is implemented which prevents the parameter ensemble spread to fall below a prescribed threshold. Since the domain used in the data assimilation experiments is small, no covariance localization is used in the estimation of the state variables or parameters.
Given that in the ensemble Kalman filter the distribution of ash concentration and parameters is assumed to be Gaussian, negative ash concentration or unphysical parameter values can result from the assimilation of observations. These unphysical 10 solutions must be corrected before using the analysis ensemble as initial conditions for the next ensemble forecast cycle. For ash concentration, negative values are turned into zero concentrations. In the case of eruption source parameters, unphysical values are replaced by a random sample from the analysis ensemble mean and covariance. If the randomly sampled value is still outside the physically meaningful range for the parameter, the process is repeated. The physically meaningful range for model parameters is set to 0-20 km and 0-15 for h and A − Suzuki respectively. 15 One of the main hypothesis of the Kalman filter is that state variables and parameters are approximately linearly correlated with the observations. This is not true for the h parameter since in the Mastin et al. (2009) emission scheme the source strength is proportional to the fourth power of h. For this reason, instead of estimating h, we estimate h 4 so that the estimated parameter is more linearly correlated with the observations. In this work, several experiments are performed to study the convergence of the filter and its sensitivity to some key parame-20 ters. Two experiments are performed using the constant parameter nature run to assess filter convergence. The first experiment starts with source parameters that are higher than the true value and will be referred as CONSTANT-UPPER, the second starts with an under estimation of the source parameters and will be referred as CONSTANT-LOWER. The initial parameter spread   Table 2. Finally, a set of simulation experiments are carried out using a larger domain to evaluate the impact of the optimized parameters upon the 5 simulation of the ash cloud farther from the vent.

Performance metrics
The evaluation of the FALL3D-ETKF system is achieved by comparing the 3D ash concentration forecast (and analysis) against the nature run, and also by measuring the consistency between the estimated and the actual forecast uncertainties. The comparison is based on the RMSE, error bias and the ensemble spread of either the forecast or the analysis which are given by 10 the following expressions: where x f,i is either the forecast or analysis ensemble mean ash concentration at time and location i and x

Constant emission profile experiments
In these experiments we explore the impact of data assimilation and parameter estimation in the steady parameter scenario.  Fig. 4 shows also the parameter ensemble spread. In these experiments, the ensemble almost always contains the true parameter value, meaning that the parameter uncertainty is well captured by the ensemble. However, it should be noted that, in these experiments, the ensemble spread 5 of the model parameters is prescribed a priori to a value which may not be the optimal one under different conditions (e.g. if the optimal parameters are time dependent or if other sources of uncertainty like errors in the atmospheric circulation, are present ). Sensitivity experiments to the parameter ensemble spread will be discussed in the following sections.  The fact that a biased parameter ensemble (i.e. the underestimation of h observed in Fig. 4a) produces a less biased estimation of ash concentrations (Fig. 5c) may be related to the non-linear relationship between h and the total mass emission at the source.
Since the emitted mass depends on h 4 , positive perturbations in h are associated with a much larger emission rate and are thus farther from the observations than ensemble members with negative perturbations in h. This creates a bias in the estimation 30 of the concentrations because, even if the ensemble is centered at the true h value, positive perturbations are farther from observations than the negative ones and therefore the ensemble mean tends to over-estimate concentrations. ETKF tries to compensate for this effect converging to a slightly biased parameter set which reduces the error bias and the RMSE.
As observed in Figure 5d, the analysis error in ash concentration is below the forecast error. This indicates that the ETKF method is efficient in reducing the error in the 3-D concentration field based in the information provided by a 2-D observation.
This is a remarkable result in a context where most of the observations are 2-D whereas operational requirements are 3-D. This finding will be the base to use the analysis as a better diagnose of the state of the plume to improve the forecasts. The reason behind this lies in the structure of the forecast error covariance matrix, which is estimated from the ensemble of forecasts.
This matrix contains the information about the covariances between mass loading (which is the observable quantity) and the concentration at different heights from which the mass loading is obtained and which are not directly observed. In this work, 5 reliable covariances between 3D ash concentrations and mass loadings are obtained by taking into account the uncertainties associated with the source parameters.

Time dependent emission experiments
These experiments use the observations simulated from the nature run with time-varying parameters (Fig. 3). The parameter ensemble is initialized with a mean h of 11 km and a mean A − Suzuki of 7 and standard deviations of 0.5 km and 2.0 10 respectively. Figure   are larger during the periods of higher h, leading to a larger spread in the concentration field. Fig. 7d shows the spatially averaged reduction in the RMSE for the total ash concentration between the forecast and the analysis. The RMSE is reduced between the forecast and the analysis at almost all vertical levels, indicating that the vertical covariance structure between mass loadings and ash concentrations at different levels is well estimated leading to accurate 3D ash concentration estimations. 30 In order to assess the impact of treating the parameters as a time dependent variable, this experiment is compared with an experiment in which data assimilation is performed but only the ash concentration field is updated. In this case, source parameters are keep constant in time at a value which is equal to the time average of the true parameters (CONTROL-NOEST, Fig. 6). This value is chosen in order to obtain a solution that is as close as possible to the one obtained with the time dependent parameters. Fig. 7 shows that the forecast RMSE and bias in the 3D ash concentration is almost always larger in the CONTROL-NOEST experiment with respect to the CONTROL experiment. The error in the CONTROL and CONTROL-NOEST experiments becomes similar around day 3 and after day 8 because at those time instants the source parameters are close to each other ( Fig. 6). Moreover, the ensemble spread for the CONTROL-NOEST experiment is almost constant in time and, because of that, changes in the forecast uncertainty are not captured (Fig. 7b). This is because time variations in the ensemble spread are 5 mainly associated with changes in the mean values of parameters. These experiments suggest that performing data assimilation for the estimation of 3D ash concentrations is not sufficient to properly constrain 3D ash concentration values and that source parameters have also to be taken into account, particularly close to the source where these parameters rapidly impact on concentrations.
As an example, Fig. 8 shows the ensemble forecast mean for the CONTROL and CONTROL-NOEST experiments and for 10 the nature run at FL200 at the 12 th assimilation cycle. The ash concentration pattern at this particular level is well represented by the simulation that estimates the source parameters whereas, in the CONTROL-NOEST experiment, there is a significant underestimation of the concentrations due to the underestimation of the column height at this particular time. Note that data assimilation is being performed to correct the 3-D ash concentrations in both experiments.

Sensitivity experiments 15
This section discusses the sensitivity of the analysis and the forecast to the parameter ensemble spread, the ensemble size, and the observation uncertainty. The purpose is to identify which are the potentially more important tuning parameters for the optimization of the system and how robust the system is with respect to errors in observations, which are known to exist in satellite-based ash mass loading estimations.
To explore the sensitivity to the parameter ensemble spread, the experiments CONTROL, HI-SPREAD and LOW-SPREAD 20 with different parameter spreads (Table 2) are compared. Figure 9 shows the estimated h obtained in these experiments as well as the total ash concentration RMSE and bias. As observed, the CONTROL experiment gives the more accurate estimation of h and the minimum RMSE and bias. When the parameter ensemble spread is larger than in the CONTROL experiment, parameter values are systematically under-estimated. As previously discussed, this can be explained by the non-linear dependence between h and the total emitted mass. However, what is relevant from this experiment is that increasing the ensemble spread 25 degrades the quality of the estimation and increase the impact of non-linearities. Higher dispersion in h increase the magnitude of positive h perturbations leading to a larger error bias, particularly above and below the maximum concentration (Fig. 9c). In the case of the LOW-SPREAD experiment, results are closer to the CONTROL experiment. However, this experiment shows a slower convergence with larger h estimation errors during the first days of the experiment. Slower or lack of convergence is expected when the parameter uncertainty is under estimated. In this case the ETKF does not allow for large corrections in 30 the parameter values based on the observations, basically because the error in the parameters is assumed to be small. These experiments show that the system is particularly sensitive to the parameter ensemble spread that has to be specified a priori. Moreover, in these idealized experiments, the optimal parameter ensemble spread is determined by the uncertainty in the observations and with no information regarding the changes of the true parameters in time.
As discussed in Section 4.1, parameters are estimated based on their covariance with the observed quantities. In the ensemble based data assimilation methods these covariances are estimated directly from the ensemble, so they can be affected by sampling errors. To asses the impact of these sampling error on the analysis, quality assimilation experiments with different ensemble sizes have been performed. Three experiments with 8, 16 and 32 ensemble members are presented (ENS-8, ENS-16 and CONTROL respectively). Fig. 10 shows the results in terms of h estimation and total ash concentration RMSE and bias. The 5 CONTROL experiment shows more accurate h estimation and, consistently, lower RMSE and bias values. However, results are not very sensitive to the size of the ensemble. The lack of sensitivity to the ensemble size might be surprising, particularly considering that no spatial localization is being used in order to reduce the impact of sampling errors. However, note that in this case the only source of uncertainty in the system comes from the uncertain parameters. Based on this, the effective dimension of the space in which uncertainties has to be constrained is two. This is confirmed by the strong covariances that exist between 10 the parameters and ash concentration within the domain (not shown). This effective low dimensionality is reinforced by the fact that the domain is small and close to the source and, because of that, ash concentration at most grid points is strongly correlated with the value of the uncertain source parameters.
The last sensitivity experiment looks into the issue of observation errors in satellite retrievals of mass loadings. In the experiments presented so far, the standard deviation of the observation errors has been assumed to be 15% of the mass loading increased. This behavior is probably a consequence of the non-linear effects arising from the non-linear relationship between 25 h and ash emission rate that has been previously discussed.

Ash concentration simulations in an extended domain simulation
The experiments discussed so far have been performed in a relatively small domain surrounding the vent. In most applications however, it is expected that forecasts over larger domains are required. In this section we explore the adequacy of the parameter estimation approach to generate a good estimation of ash dispersion over larger domains in an idealized context in which the 30 atmospheric flow is perfectly known. To this purpose a nature simulation over a larger domain (Fig. 12) is performed. This nature run is forced with the same evolution of parameters of the time dependent parameter nature run and spanning the same time period.
To see if the estimated parameters can be used to reconstruct the ash cloud far from the source, the ensemble mean estimated parameters from the CONTROL experiment are used to produce an estimation of the ash cloud. Figure 12 shows a snapshot of the nature run and the experiment forced with the estimated parameters. Note how both ash clouds are very similar even far from the source, indicating that the estimated parameters are sufficient for the reconstruction of the ash plume. Figure 12 also shows ash distribution at 12 and 24 hours forecast times. These forecasts are initialized from the model 5 run forced with the estimated parameters and use the estimated parameters at the initialization time from the CONTROL run.
Parameters remains constant during the forecast since there is no predictive model available for these parameters. Fig. 12 shows that there is a good agreement between forecasts and nature run. This suggests that initializing a forecast from a long run forced with the optimized parameters can be a cost-effective strategy to generate short lead ash concentration forecast over a relatively large domain.

10
Although these results are encouraging, it should be taken into account that in more realistic situations, other sources of uncertainty (e.g. uncertainty in the flow or model errors) can significantly affect the evolution of the ash plume far from the source. In this case, the forecast quality can suffer from the estimation of the 3D ash concentration over the entire domain based on the assimilation of mass loading observations.

15
In this work we have presented an ensemble-based data assimilation system coupled with the FALL3D ash dispersion model that is able to use ash mass loading observations to simultaneously estimate the 3D ash concentration field and to constrain source term parameters. An OSSE preliminary evaluation of the system has been presented and some sensitivity tests performed. The experiments focused on two FALL3D model parameters, one that defines the vertical emission profile and the eruptive column height (and related emitted mass). The ETKF-FALL3D system shows a robust convergence to the optimal 20 parameters and an accurate reconstruction of the 3-D ash concentrations based on noisy and partial 2-D ash mass loading observations.
Estimation of time-dependent source parameters has been successful within the OSSE context. The ensemble not only produces an estimation of the covariances between the observed variables and the parameters, but also provides a time-dependent estimation of the forecast uncertainty which resembles the time evolution of the forecast errors. The strong time variability of 25 the ensemble spread is mainly associated with the relationship between column height and emitted mass.
Sensitivity experiments has been conducted to investigate how the parameter ensemble spread, the ensemble size and the observation errors affect the results. The parameter ensemble spread produces a significant impact on the quality of the estimated concentrations and parameters. Larger ensemble spreads lead to stronger biases, both in concentrations and parameters, whereas lower ensemble spreads produce an over-confident ensemble and slower converge rates that degrade the estimation 30 results. It is important to note that the optimal parameter ensemble spread can depend on the time variability of the estimated parameters and on other sources of uncertainty like errors present in the observations and in the model. The sensitivity to the ensemble size revealed that, even for this low dimensional estimation problem, ensemble sizes up to 32 members show some improvement with respect to ensembles of 16 and 8 members although the impact of increasing the ensemble size is smaller than the impact associated with changes in the parameter ensemble spread.
The sensitivity to the observation errors shows a particular behavior with an increase in systematic errors both in the parameters and in the concentrations with increasing observational errors. When observation errors reach 40% of the true ash loadings, the estimated parameters fail to converge during the first days of the experiment leading to significantly larger errors 5 in the ash concentration forecasts.
Experiments presented in this work are limited to a small domain surrounding the vent. Experiments on a larger domain, show that the optimized parameters can be used to force an ash dispersion simulation that can reproduce the ash cloud properties far from the vent as long as the atmospheric circulation is accurately known. This simulations can be used to initialize ash dispersion forecasts over larger domain as a computationally cheaper alternative to running a data assimilation system with