Particle filter based volcanic ash emission inversion applied to a hypothetical sub-Plinian Eyjafjallajökull eruption using the chemical component of the Ensemble for Stochastic Integration of Atmospheric Simulations (ESIAS-chem) version 1.0

A particle filter based inversion system to derive timeand altitude-resolved volcanic ash emission fluxes along with its uncertainty is presented. For the underlying observation information only vertically integrated ash load data as provided by retrievals from nadir looking imagers mounted on geostationary satellites is assimilated. We aim to estimate the temporally varying emission profile with error margins, along with evidence of its dependencies on wind driven transport patterns within variable observation intervals. Thus, a variety of observation types, although not directly related to volcanic ash, can 5 be utilized to constrain the probabilistic volcanic ash estimate. The system validation addresses the special challenge of ash cloud height analyses in case of observations restricted to bulk column mass loading information, mimicking the typical case of geostationary satellite data. The underlying method rests on a linear-combination of height-time emission finite elements of arbitrary resolution, each of which is assigned to a model run subject to ensemble-based space-time data assimilation. Employing a modular concept, this setup builds the Ensemble for Stochastic Integration of Atmospheric Simulations (ESIAS-chem) 10 that comprises a particle smoother in combination with a discrete-grid ensemble extension of the Nelder-Mead minimization method. The ensemble version of the EURopean Air pollution Dispersion Inverse Model (EURAD-IM) is integrated into ESIAS-chem but can be replaced by other models. The performance of ESIAS-chem is tested by identical twin experiments. The application of the inversion system to two notional sub-Plinian eruptions of the Eyjafjallajökull with strong ash emission changes with time and injection heights demonstrate the ability of ESIAS-chem to retrieve the volcanic ash emission fluxes 15 from the assimilation of column mass loading data only. However, the analysed emission profiles strongly differ in their levels of accuracy depending of the strength of wind shear conditions. Under strong wind shear conditions at the volcano the temporal and vertical varying volcanic emissions are analyzed up to an error of only 10 % for the estimated emission fluxes. For weak wind shear conditions, however, analysis errors are larger and ESIAS-chem is less able to determine the ash emission flux variations. This situation, however, can be remedied by extending the assimilation window. In the performed test cases, the 20 ensemble predicts the location of high volcanic ash column mass loading in the atmosphere with a very high probability of >95 %. Additionally, the ensemble is able to provide a vertically resolved probability map of high volcanic ash concentrations to a high accuracy for both, high and weak wind shear conditions. 1 https://doi.org/10.5194/gmd-2021-30 Preprint. Discussion started: 10 May 2021 c © Author(s) 2021. CC BY 4.0 License.


Introduction
Apart from water and carbon dioxide, volcanic eruptions release enormous amounts of harmful trace gases and particulate 25 matter. Chemistry transport models have limits in estimating the emission strength of those events in an accurate temporal and spatial resolution because the emission characteristics, i. e. the emission strength and its temporal and vertical distribution within the eruption plume, are typically not well known (e. g., Bursik et al., 2012). Therefore, special assessment methods for unexpected emission events are necessary. Because of the potential threat to human health, economy, and climate of volcanic eruptions, a detailed knowledge not only about the spatial and temporal variations of the emissions and its strength is required 30 but also accurate information about the analysis error.
Typically, volcanic eruptions occur as sequences of emissions with highly varying ejection mass and height. Only limited observations of volcanic ash emission parameters are available (e.g. eruption plume heights retrieved form radar measurements, Arason et al. (2011)). Thus, eruption models are applied to simulate volcanic emissions instead. These may be inferred from statistical data (e. g., Sparks et al., 1997;Mastin et al., 2009) or physical processes (e. g., Woodhouse et al., 2013;Folch et al., 35 2016). Statistical models are based on observational data from only a few, highly heterogeneous volcanic eruptions, while physical plume-scale models require vent and magma details, which are poorly known, and thus making these models highly uncertain.
Another potential way to constrain volcanic ash emissions is the use of observations of volcanic ash in the atmosphere. Advanced numerical analysis techniques for quantitative and stochastic estimation of volcanic ash concentrations and emissions 40 use mostly satellite observations of column mass loading via data assimilation methods (e. g., Wilkins et al., 2016a). In contrast to lidar observations, e. g. from the Cloud Aerosol LIdar with Orthogonal Polarization (CALIOP) instrument on board CALIPSO satellite (Winker et al., 2009) or the ground (e. g. via the European Aerosol Research Lidar Network, EARLINET 1 ), column mass loading observations do not provide any information about the vertical distribution of volcanic ash. Therefore, multiple data assimilation methods make assumptions about the vertical extent of the volcanic ash cloud (e. g., Schmehl et al.,45 2012; Wilkins et al., 2016b;Zidikheri et al., 2016). In addition to column mass loading observations of volcanic ash clouds, Zidikheri et al. (2017a) suggested to use brightness temperature measurements to distinguish regions with high mass load from those with low mass load. The great advantage of using column mass loading observations as available from, for example, the Spinning Enhanced Visual and InfraRed Imager (SEVIRI) on board Meteosat Second Generation (Schmetz et al., 2002) is the horizontally more complete picture of the volcanic ash extent. Observations of column mass loading and lidar measurements 50 are restricted to areas without water cloud cover. Finally, the challenge of estimating the three-dimensional volcanic ash field from two-dimensional volcanic ash column mass loading observations remains to be solved.
First estimations of volcanic ash emissions from the 2010 Eyjafjallajokull eruption in a high temporal and vertical resolution were made by Stohl et al. (2011) and later by Kristiansen et al. (2012) and Kristiansen et al. (2015). Their algorithm bases on the inversion technique of Eckhardt et al. (2008), in which an optimal combination of distinct emission packages is estimated 55 using a least squares method. Further developments of assimilation systems for volcanic ash emissions include the "genetic 1 https://www.earlinet.org/index.php?id=earlinet_homepage 2 https://doi.org/10.5194/gmd-2021-30 Preprint. Discussion started: 10 May 2021 c Author(s) 2021. CC BY 4.0 License. algorithm variational approach" (Schmehl et al., 2012) and an adjoint free ensemble version of the 4d-var data assimilation method (Lu et al., 2016). Zidikheri et al. (2016) and later Zidikheri et al. (2017b) developed an assimilation system that aims to analyze the horizontal distribution of volcanic ash column mass loading rather than the emission strength. This was extended by Zidikheri et al. 60 (2017a) to additionally estimate the height and the particle size distribution of volcanic ash emissions using a parameter refinement method. Wilkins et al. (2014) used the "data insertion" method in which observed volcanic ash column mass loading act as virtual sources for volcanic ash with a predefined vertical distribution. The algorithm was successfully applied to the eruptions of Eyjafjallajökull, Iceland, 2010 (Wilkins et al., 2016b) and Grímsvötn, Iceland, 2011(Wilkins et al., 2016c. A general framework for calculating uncertainties of volcanic ash concentrations for constant volcanic ash emissions given 65 "any model and any observational data" was proposed by Denlinger et al. (2012). Bursik et al. (2012) applied the "polynomial chaos quadrature weighted estimate" (PCQWE) method to volcanic ash emissions. This approach was further extended by Stefanescu et al. (2014) to take uncertainties in the wind fields into account as well. An extension of the polynomial chaos quadrature method was proposed by Madankan et al. (2014) to generate hazard maps of volcanic ash in the atmosphere. The PCQWE method aims to map uncertainties in the input parameters of a volcanic eruption onto the volcanic ash concentrations 70 without accounting for constraining volcanic ash observations. Additionally, Dare et al. (2016) investigated the influence of meteorological ensemble forecasts on the dispersion of volcanic ash. They found that not only the ensemble statistics should be evaluated but also may the single ensemble members contribute significant information to the distribution of volcanic ash.
Although these studies applied highly advanced data assimilation methods for analyzing the emission strength and its uncertainty, a joint assessment of both, emission strength and its uncertainty, in a high temporal and vertical resolution has not yet 75 been evaluated. Thus, this contribution aims to fill this gap by providing such estimates of emission strength and its uncertainty in a high temporal and spatial resolution resulting in height-resolved probability maps of volcanic ash concentrations.
Section 2 describes the full stochastic assimilation system ESIAS-chem and the methods applied. The potential and limitations of ESIAS-chem are shown by identical twin experiments in section 3 for predictions of the ensemble mean volcanic ash concentration and the posterior probability of high volcanic ash content. Conclusions will be given in section 4. The Ensemble for Stochastic Integration of Atmospheric Simulations (ESIAS) is designed to control simultaneous and interactive runs of ultra-large ensembles of complex atmospheric models, with special dedication to numerical weather prediction models and chemistry transport models, along with extensions to include data assimilation modules. There exist a meteorological (ESIAS-met) and an atmospheric chemical part (ESIAS-chem). The emphasis of this paper is placed on the latter for 2.1 Ensemble of emission packages 90 The main assignment of ESIAS-chem is to simulate normalized emissions as a set of N emis pairwise distinct emission packages (emission scenarios in the terminology of Stohl et al. (2011)) each of which with a unit mass of ash. More specific, the eruption plume is discretized into N T time steps and K max vertical layers. Emission strength and plume height of a volcanic eruption may vary quickly, which may result in multiple maxima in the vertical distribution of the volcanic ash emissions within the discretized emission profile. This is accounted for in our approach. Each ensemble member simulates the volcanic ash 95 concentration released by a single emission package for an individual time step and height interval. A similar approach for estimating the volcanic ash column mass loading was used by Stohl et al. (2011) and Kristiansen et al. (2015) aiming to estimate the optimal emission profile but not its uncertainty.
The optimal combination of the emission packages is calculated by minimizing the cost function J(a) a = argmin (J(a)) = argmin x k ∆z k , where x k is the modeled concentration of volcanic ash, given in [µgm −3 ] and ∆z k is the thickness of model layer z in [m]. It is noticed that the vertical and temporal resolution of the emissions can be varied by changing the parameters K max and N t , thus, making them independent from the vertical and temporal resolution of the model. Although the algorithm is designed for column mass loading observations, it is as well applicable to observations of any volcanic ash related quantity, e. g. brightness 110 temperature as proposed by Zidikheri et al. (2017a).

Nelder-Mead algorithm
The minimization problem posed by (1) is solved using the Nelder-Mead algorithm (Nelder and Mead, 1965). The Nelder-Mead minimization algorithm is proved to be robust, especially in cases where the function to be minimized has discontinuities or the function values are noisy (see McKinnon, 1998), which is expected to be likely in highly variable volcanic eruptions.
The idea of the algorithm is to move a simplex on the surface of the cost function to find its minimum in a N -dimensional space. The version of the Nelder-Mead method used in this study follows Gao and Han (2012) and utilizes adaptive parameters controlling the step size for each iteration of the minimization. The version has been implemented for parallel operation (Klein and Neira, 2014;Lee and Wiswall, 2007).
The minimization was performed in N N 0 , which has been found to be more effective than the minimization in R N . Thus, only integer values are accepted for the scaling factors a i of the emission packages. By applying this constraint it is assumed that the introduced errors are of lower order than the error introduced by the temporal resolution of the emission packages. Further, in order to generate an ensemble of analyses and because the minimization depends on the initial simplex, the minimization tasks 125 are performed for each ensemble member with different initial simplices. Thus, the minimization algorithm is called hereafter discrete-grid ensemble Nelder-Mead method (DENM).

Particle filter
The particle filter methodology, also known as sequential importance resampling, is used as a non-Gaussian data assimilation method for large ensemble simulations. The particle filter method was proposed by Gordon et al. (1993) and further popularized 130 in oceanography and meteorology by van Leeuwen (2009). It develops from Bayes' Theorem where p(·) denotes the probability density function (PDF), y the observations, and x the model state. The a priori PDF is approximated by an ensemble of N ens model runs

135
where δ(·) denotes the Kronecker delta function and x i is the model state of particle (ensemble member) i.

Applying (4) results in
where an individual weight w i is applied to each ensemble member. Thus, each ensemble member is weighted by the normalized likelihood of its current model state. It is noted that in the particle filter method no assumptions of the error statistics of 140 the model state and the observations were made. Further, the particle filter method is directly applicable to nonlinear models.
In particle filters, filter degeneracy often occurs (cf. Snyder et al., 2008;Bickel et al., 2008), especially in high dimensional problems. Several methods exist to reduce filter degeneracy (see e.g. van Leeuwen (2009) for a review) and the reader is referred to the original papers for more information. In this model setup for improving the ensemble skill, residual resampling (Liu and Chen, 1998) is used, in which the ensemble members with high weights are duplicated and perturbed, 145 replacing ensemble members with vanishing weights. Thus, the statistical meaning of the ensemble is preserved.

ESIAS-chem
ESIAS-chem is designed as a flexible analysis system for quantitative volcanic ash assessments, along with a quantification of the uncertainty of the optimal emission flux profile. Further, it is applicable to other accidentally released matter and con- stituents, suitably constraining observations are available. In order to account for meteorological uncertainties, ESIAS-chem 150 is capable to be coupled with ensembles of meteorological fields. The workflow of the system is illustrated in Fig. 1. Once a volcanic eruption is detected (t = t 0 ), the system is started by generating the ensemble of emission packages (cf. section 2.1) of default mass of volcanic ash.
As long as no observations are available, this ensemble serves as an estimation predictor of the maximum possible volcanic ash extension without providing quantitative volcanic ash estimates. Once new observations become available, the system is 155 restarted at time t = t 0 . The previously calculated ensemble of emission packages is reused, now integrated forward in time until the observation time (t = t i+1 ). Further ensemble members are included to account for the latest emissions in the interval The ensemble is compared with the observations to calculate (1), which is to be minimized using the DENM method. As this algorithm is optimized for the estimation of column integrated ash loading in a considerably underdetermined control system, 160 a regularization term is added to the cost function (1), leading to This choice restricts the scaling factors a to vary too strongly. As each volcanic eruption has its own individual emission characteristics, no fixed assumptions have been made for matrix B, which accounts for correlations of the individual emission packages. Those emission characteristics may be (partial) eruption column collapse, abrupt change in emission strength, or 165 multiple vents within one model grid box erupting in different strength and height. Here, B is chosen as diagonal matrix B = diag(10) but other choices may also be possible. This value for B was found to be optimal for this study by several tests.
Thus, the regularization term was chosen in order to maintain a suitable spread of the analysis ensemble. To account for 170 common eruption column profiles, it has been further tested to start the minimization with an ensemble of predefined umbrellashaped emission profiles of different emission strength and temporal and vertical distributions (not shown). At least for the current study it has been found that the analysis using these predefined emission profiles performed worse than using arbitrary emission profiles, mainly because in the former case the minimization was too strongly restricted by the predefined emission profile shape and suppressed variations in the emission profile. Thus, the current analysis is performed using an arbitrary set 175 of initial emission profiles. In addition, the observation errors are represented by perturbed observations (cf. Houtekamer and Mitchell, 1998), which are assimilated by ESIAS-chem leading to larger ensemble spreads.
Finally, a particle filter step is applied. The weights, which result from the filtering step, are applied to the optimized emission profiles. If new observations become available the assimilation window may be elongated to use the new observations for updating the emissions within the whole assimilation window [t 0 , t i+1 ].

Identical twin experiments
The utility of the system to provide quantitative estimates of the volcanic ash emission uncertainty is explored by identical twin experiments. Identical twin experiments are necessary, yet not sufficient standard test procedures for validating spatio-temporal data assimilation and inverse modelling set-ups. They rest on the "perfect model assumption" and its analog for the data side: exactly known accuracy and representativity. This provides a total knowledge of the "synthetic truths" as given by simulations 185 with the same model and extraction of artificial "measurements/soundings" thereof. Given the identical twin assumption the experiment is then to be made realistic in all other respects, as the two different weather conditions on our case. Daley (1991) concludes however, that identical twin experiments "err on the optimistic side" (loc. cit.). Yet, the applicability of ESIAS-chem to real volcanic eruptions will be shown in a companion paper. In this study, ESIAS-chem is coupled to the regional air quality  (1) is scaled by the observation error covariance matrix, a relative error of low volcanic ash values overemphasizes the system to these observations of low volcanic ash content. Further, no relative error for observations of no-volcanic ash observations can be obtained. Therefore, it was decided to use the following expression for the observation error σ y i with a minimum error of 0.1 gm −2  The performance of ESIAS-chem for different assimilation window lengths is tested for a hypothetical eruption of the Icelandic volcano Eyjafjallajökull under different, yet real meteorological conditions. The Hovmoeller-like plot of the nature run emission profile is shown in Fig. 2. It shows the variable emission rate by a height-time graphic above the volcano. The se-215 lected sub-Plinian eruption type (Bardintzeff and McBirney, 2000) is characterized by two short explosive phases between 02-04 UTC and 06-08 UTC reaching a height of approx. 8 km above the volcano.
The length of the assimilation window influences the performance of the data assimilation algorithm due to the influence of vertical and horizontal mixing and vertical wind shear. Hence, the performance of ESIAS-chem is tested for different assimilation window lengths. All assimilation windows start at 00 UTC for the specific day and last for 06-36 hours. Fig. 3 illustrates 220 the different assimilation window lengths. The assimilation windows differ in length by multiples of 06 hours. Certainly, by increasing the assimilation window length the observations include more information, as the residence time of volcanic ash in the atmosphere is increased. Contrary, vertical and horizontal mixing of volcanic ash emitted may limit the benefit that is gained by increasing the assimilation window length. Once the volcanic ash emissions are optimized, a forecast is appended until 36 hours after the simulation start (hatched areas in Fig. 3). Thus, beneficial impacts of the assimilation results for the 225 analyses with differing assimilation window lengths can be assessed.
The first real weather test day to which ESIAS-chem is applied is 15 April 2010, which was characterized by strong westnorth-westerly winds at Iceland. This is illustrated by the wind profile at the volcano (Fig. 4a) for the whole simulation length of 36 hours and the wind field over Europe at approx. 5 km height on 15 April 2010, 12 UTC (Fig. 4b). During this day, the polar front and the polar jet stream are located above Iceland, driving the volcanic ash to travel fast southeast towards 230 continental Europe, with wind speeds of up to 60 ms −1 at heights of 5-8 km.  thetical sub-Plinian eruption of the Eyjafjallajökull volcano on 29 April 2010 (Fig. 5). The same emission profile is taken as depicted in Fig. 2. This day is characterized by weak winds of approximately 10 ms −1 in the vicinity of the volcano, which is illustrated by Fig. 5. Thus, the emitted volcanic ash is only slowly transported. Additionally, the vertical wind shear on 29 April 2010, 12 UTC, is low, as indicated by a higher barotropicity above Iceland (Fig. 5c). The two dates were chosen because of their different wind patterns and the real eruption of the Eyjafjallajökull that occurred during these days. Hence, the identical 240 twin experiment provides an optimal case scenario for the application of ESIAS-chem to real volcanic eruptions.
In order to validate the results of the stochastic assimilation method, first the volcanic ash cloud of the ensemble mean is compared to the nature run. Two measures are considered, which provide information about the accuracy of the horizontal extent of the volcanic ash cloud and the distribution of volcanic ash within the cloud. The former is investigated by the pattern correlation coefficient (pcc, cf. Zidikheri et al., 2016) of the column mass loading of the nature run and the ensemble mean.

245
The latter is analyzed by the relative mean absolute error (RMAE) of the column mass loading.
The pcc is defined by (Zidikheri et al. (2016)) with va = va − va. Herein, the entries of the volcanic ash detection vector va for the model mean (subscript x) and observations (y) are equal to 1 if the grid column contains volcanic ash above 0.2 gm −2 , which is the detection limit of volcanic ash 250 column mass loading observations , and 0 otherwise. The mean volcanic ash detection va is calculated by va =< 1, va > / < 1, 1 >, where 1 denotes the vector with 1 on all entries and < ·, · > indicates the scalar product. The pattern correlation coefficient gives information about the compliance between the simulated and observed horizontal distribution of volcanic ash in the atmosphere.

255
If the volcanic ash clouds, indicated by the column mass loading, of the nature run and the ensemble mean perfectly coincide, the pattern correlation coefficient equals 1. If the volcanic ash clouds of the nature run and the ensemble mean are totally disjoint, the pattern correlation coefficient equals 0. It is noted that for this model setup with perfect meteorology, pcc < 1 indicates that the analysis contains volcanic ash either in model layers or at times, where no volcanic ash is emitted in the nature run.

260
The RMAE is defined by where N y is the number of grid columns in which volcanic ash column mass loading of the nature run y j and the ensemble mean x j are ≥ 0.2 gm −2 . The RMAE is as well calculated for volcanic ash concentrations, for which this threshold is ≥ 10 µgm −3 . In this case, y j corresponds to the volcanic ash concentration simulated by the nature run. The RMAE measures 265 the relative difference of column ash mass loads (or volcanic ash concentrations) between nature run and ensemble mean, averaged over all grid cells. Higher values of the RMAE are a result of different height-time-mass emission patterns between the nature run and the ensemble mean, given the assumed perfect meteorology used in this analysis.
The pattern correlation coefficient is shown in Fig. 6 for the two analysis days. The lines in Fig. 6 indicate results for different assimilation window lengths as illustrated by Fig. 3. Fig. 6 shows a constantly large pattern correlation coefficient > 0.95 270 after the artificial eruption terminated at 08 UTC for both analysis days, except for assimilation window lengths of 06 and 12 hours. By applying an assimilation window of 06 hours from the simulation start, the artificial volcanic eruption has not terminated, thus, the latest emissions from the nature run volcanic eruption are not considered in the analysis. This leads to a reduced pcc for the 06 hour assimilation window test case. The assimilation window of 12 hours ends four hours after the termination of the artificial volcanic eruption. The reduced pcc exhibits that this short assimilation window is not sufficient in Again, the pattern correlation coefficient does not account for deviations in the strength of volcanic ash column mass loading at locations in which the ensemble mean and the nature run differ in volcanic ash load.

280
Increasing the assimilation window length (i.e. taking later observations into account) increases the pattern correlation coefficient on both days. The analysis suggests that for the respective test cases an assimilation window of 18 hours, that is 10 hours after the artificial eruption terminated, is sufficient for ESIAS-chem to analyze the exact location of the volcanic ash cloud as observed from space leading to a pcc value that remains high (> 0.95) throughout the full analysis time period of 36 hours. Fig. 6 demonstrates that the assimilation system is able to accurately analyze the horizontal extend of the volcanic ash cloud.  test cases with longer assimilation windows (≥ 18 hours) show a RMAE of the order of 10 % for both days. These low values show the good performance of the analysis for these assimilation window lengths with respect to the nature run. In principle, Fig. 7 corroborates the same findings that are analyzed for the pattern correlation coefficient, i. e. increasing the assimilation 290 window length decreases the error of the ensemble mean. For both days and meteorological circulation patterns, an assimilation window of 18 hours is sufficient to reduce the RMAE to a value of approx. 10 % for column mass loading values above 0.2 gm −2 . On 15 April 2010 assimilation windows larger than 24 hours result in a slightly higher RMAE than the analysis using an assimilation window of 18 hours. This is a result of the convergence of volcanic ash in the upper troposphere south of Norway around 24 hours after the simulation has started (not shown). Thus, additional observations at later times do not 295 contribute significant information to the assimilation system. In summary, Fig. 7 proves that the assimilation system is able to analyze the distribution of volcanic ash column mass loading properly for weak and strong wind conditions.
The above analysis focuses on the comparison of the nature run and the ensemble mean with respect to column mass loading of volcanic ash. Thus, it does not provide any information about the vertical distribution of volcanic ash. The ability of ESIASchem to infer vertical profiles of volcanic ash is given in Fig. 8, which displays the relative mean absolute error of the volcanic 300 ash concentrations above 10 µgm −3 . The RMAE of the volcanic ash concentrations decreases by increasing the assimilation window length, which is especially visible for 29 April 2010. On both days, an assimilation window of only 6 hours results in a RMAE larger than 100 %. Therefore, this test case is not shown in Fig. 8. On average, the RMAE reduces to about 20 % on both days for assimilation windows larger than 18 hours, which shows the good performance of the ESIAS-chem analysis not only in terms of column mass loading but also in terms of the vertical distribution of the volcanic ash in the atmosphere. use. Fig. 9 and Fig. 10 display the profile of (a) the nature run emissions, (b) the ensemble mean emissions, (c) the relative error of the ensemble mean and (d) the relative ensemble standard deviation for the 24 hour assimilation window experiments on 15 April and 29 April 2010. Herein, x and y are the ensemble mean and nature run emissions, respectively, and σ x is the ensemble standard deviation.

315
The total nature run emissions on both days (4.25 · 10 8 tons and 4.30 · 10 8 tons on 15 April and 29 April, respectively) are well captured by the analyzed total emissions with mean emissions of 4.58 · 10 8 tons and 4.10 · 10 8 tons, respectively, and standard deviations of 3.67 · 10 7 tons and 3.47 · 10 7 tons. The relative error of the total emitted volcanic ash is 7.7 % and 4.7 %, respectively. On 15 April 2010, the analyzed emission profile of the ensemble mean shows the two explosive eruptions of the nature run emission profile with the correct height of the maximum emissions at the right time (Fig. 9b). Even though the 320 ensemble mean shows a vertically and temporally smoothed emission profile, the false emissions are low with respect to the maximum emissions. The relative error of the ensemble mean emissions is of the order of 10 %-20 % for most emission times and heights (cf. Fig. 9c) and therefore, the results are similar to the analysis presented before. The relative ensemble standard deviation is of the same order as the relative error of the ensemble mean emissions, indicating a reasonable ensemble spread.
The analyzed emission profile of the ensemble mean on 29 April 2010 (Fig. 10b) however shows strong deviations from the 325 nature run emission profile (Fig. 10a). Although the highest level emissions of the nature run emission profile in 8 km height  The analyzed emission profile on 29 April 2010 shows the limits of the ESIAS-chem approach. While the volcanic ash column mass loading have only low errors, the emission profile shows large deviations up to 60 unit% and more (Fig. 10c).

335
The ensemble standard deviation of the emission profile (Fig. 10d) Fig. 11b). Even though the volcanic ash in 4 km height in the center of the cross-section is covered from above by the elevated volcanic ash in 7 km height, the ensemble predicts a 50 % chance of volcanic ash at 4 km height. This is remarkably, since only vertically integrated observations of volcanic ash are assimilated.

360
The volcanic ash northeast of the center of the vertical cross-section (i. e. to the right in Fig. 11b) is predicted by only 20-30 % of the ensemble. The ensemble predicts this volcanic ash in this vertical column to be at a height of 6-7 km by a chance of more than 70 %. This may be because of the lack of vertical wind shear that prevents the distinction of volcanic ash emitted at different heights. Fig. 12a shows the probability of volcanic ash column mass loading exceeding 2 gm −2 as predicted by the ensemble on 365 30 April 2010, 12 UTC, i. e. 36 hours after the simulation start and 12 hours after the end of the assimilation window. Isolines of nature run's volcanic ash column mass loading for 0.5, 1, and 2 gm −2 are also given by blue lines. A vertical cross-section of the probability of volcanic ash concentration exceeding 2 mgm −3 along the red line in Fig. 12a is shown in Fig. 12b. Even

Conclusions
In this study, a new method for estimating volcanic ash emissions and its uncertainty from column mass loading observations is developed. This new method is realized by the atmospheric chemical part of the Ensemble for Stochastic Integration of Atmospheric Simulations (ESIAS-chem). On the one hand, the method comprises an ensemble-based particle smoother, which 375 extends the assimilation window to include the latest observations available, taken for the estimation of the emission profile for the whole assimilation window. On the other hand, the discrete-grid ensemble Nelder-Mead method (DENM) is developed in order to efficiently achieve a posterior ensemble representation of the optimum of the cost function. The particle smoother approach enables to use the latest observations for the emission within the whole assimilation window.
The presented investigation acts as a best case scenario for probabilistic volcanic ash assessments. The analysis is idealized 380 in different ways: Firstly, the uncertainties in meteorological fields, especially in winds, is neglected in the study. Secondly, the amount of observational data is exceptionally large, with observations of the full domain every 6 hours. The inclusion of observations of ash-free areas supports the systems ability to constraint the emissions to the correct times and heights. However, in a companion paper, ESIAS-chem is applied to real observations of volcanic ash. This study will also deal with the ability of the system to remove volcanic ash emissions at false height and times. Even though direct observations of volcanic ash were used in this study, ESIAS-chem is extremely flexible in terms of observational data. By estimating the emission profile the system is able to use other kind of observations not directly related to volcanic ash in the atmosphere, such as samples of tephra fall out, if available.
The flexibility of the system enables also to account for constraints on the emission profile. Those may result from radar or web cam observations (e. g., Arason et al., 2011). It is noted that the vertical and temporal resolution of the emissions may be 390 increased if it is suggested by information about the respective volcanic eruption. By the choice of an u-shaped emission profile, ESIAS-chem was challenged with highly variable volcanic ash emissions. By leaving the emission profile as unconstrained as possible, the system proves at least for strong wind conditions to estimate the emission profile even for this challenging situation.
The analysis showed that an assimilation window of 24 hours is sufficient in order to provide reliable forecasts of the vertically 395 resolved volcanic ash distribution in the atmosphere. The system verifies to capture high volcanic ash concentrations in the correct height with a high probability. Although the system lacks to estimate the true emission profile sufficiently well for weak wind conditions, the analysis of the probability of volcanic ash showed that its vertical distribution in the atmosphere is accurately predicted.
Besides volcanic ash eruptions, ESIAS-chem is applicable to a variety of emission scenarios, especially unexpected emission 400 events like forest fires and mineral dust events. Therefore, it provides a fast and efficient model for source term estimation including uncertainty representation. In principle, the method can be adapted to multi-source emission scenarios. The enhanced need for compute resources can partly be absorbed by a reduced resolution of the emission profile and will be in the focus of future work. Further, for the investigation of real volcanic eruptions, the uncertainties of meteorological variables, especially winds, need to be represented. It is noted that ESIAS-chem is flexible in integrating other modules and is applicable to other 405 atmospheric models as well. Competing interests. The authors declare that they have no competing interests.