Improving volcanic ash predictions with the HYSPLIT dispersion model by assimilating MODIS satellite retrievals

Currently NOAA’s National Weather Service (NWS) runs the HYS PLIT dispersion model with a unit mass release rate to predict the transport and dispersion of volcanic ash . The model predictions provide information for the Volcani c Ash Advisory Centers (VAAC) to issue advisories to meteorologi cal watch offices, area control centers, flight information c e ters, and others. This research aims provide quantitative foreca sts of ash distributions generated by objectively and optim ally estimating the volcanic ash source strengths, vertical distrib ution and temporal variations using an observation-modeli ng inversion 5 technique. In this top-down approach, a cost functional is d efined to mainly quantify the differences between model pred ictions and the satellite measurements of column integrated ash con centrations, weighted by the model and observation uncerta inties. Minimizing this cost functional by adjusting the sources pr ovides the volcanic ash emission estimates. As an example, M ODIS (MOderate Resolution Imaging Spectroradiometer) satelli te retrievals of the 2008 Kasatochi volcanic ash clouds are u sed to test the HYSPLIT volcanic ash inverse system. Because the sa tellite retrievals include the ash cloud top height but not t he 10 bottom height, three options for matching the model concent rations to the observed mass loadings are tested. Although t he emission estimates vary significantly with different optio ns the subsequent model predictions with the different rele ase estimates all show decent skill when evaluated against the una ssimilated satellite observations at later times. Among th e three options, integrating over three model layers yields slight y better results than integrating from the surface up to the volcanic ash cloud top or using a single model layer. Inverse tests also sh ow t at including the ash-free region to constrain the model is not 15 beneficial for the current case. In addition, extra constrai n s to the source terms can be given by explicitly enforcing “ o-ash” for the atmosphere columns above or below the observed ash cl oud top height. However, in this case such extra constraints re not helpful for the inverse modeling. It is also found that si multaneously assimilating observations at different time s produces better hindcasts than only assimilating the most recent obs ervations. 1 Atmos. Chem. Phys. Discuss., doi:10.5194/acp-2016-750, 2016 Manuscript under review for journal Atmos. Chem. Phys. Published: 13 October 2016 c © Author(s) 2016. CC-BY 3.0 License.


Introduction
Large amounts of ash particles are produced during violent volcanic eruptions.After the initial ejection momentum carrying them upwards, ash particles rise buoyantly into the atmosphere.Then, volcanic ash travels away from the volcano following the atmospheric flow.Fine ash particles may remain in the atmosphere for days to weeks or longer and can travel thousands of miles away from the source (Rose and Durant, 2009).They have severe adverse impacts on the aviation industry, human and animal health, agriculture, buildings, and other infrastructure (Prata and Tupper, 2009;Gordeev and Girina, 2014;Wilson et al., 2011;Horwell and Baxter, 2006;Wilson et al., 2012).To help prepare for and mitigate such impacts, it is important to not only monitor but also to forecast the volcanic ash transport and dispersion.
Starting with a memorandum of understanding (MOU) signed between the United States National Oceanic and Atmospheric Administration (NOAA) and the Federal Aviation Administration (FAA) in December 1988, the NOAA Air Resources Laboratory (ARL) developed a volcanic ash forecast transport and dispersion (VAFTAD) model for emergency response, focusing on hazards to aircraft flight operations (Heffter and Stunder, 1993).Currently, the NOAA National Weather Service (NWS) runs the HYSPLIT dispersion model (Draxler and Hess, 1997;Stein et al., 2015a) with a unit mass release rate to qualitatively predict the transport and dispersion of volcanic ash.The model predictions provide important information for the Volcanic Ash Advisory Centers (VAAC) to issue advisories to meteorological watch offices, area control centers, flight information centers, and others.
In order to quantitatively predict volcanic ash, realistic source parameters need to be assigned to the volcanic ash transport and dispersion models.Mastin et al. (2009) compiled a list of eruptions with well-constrained source parameters.They found that a mass fraction of debris finer than 63 µm (m63) could vary by nearly 2 orders of magnitude between small basaltic eruptions (∼ 0.01) and large silicic ones (> 0.5).Default source parameters were assigned to the world's more than 1500 volcanoes.They may be used for ash cloud modeling when few observations are available in the event of an eruption.
With the advancement of remote-sensing techniques, satellites have played an important role in detecting and monitoring volcanic ash clouds (Seftor et al., 1997;Ellrod et al., 2003;Pergola et al., 2004).An automated volcanic ash cloud detection system has been developed and continuously improved (Pavolonis et al., 2006(Pavolonis et al., , 2013(Pavolonis et al., , 2015a, b), b).In addition to detecting and monitoring ash clouds, satellite measurements allow many ash cloud characteristics to be quantified.For instance, Wen and Rose (1994) used two-band data from the NOAA Advanced Very High Resolution Radiometer (AVHRR) to retrieve the total mass of a volcanic ash cloud from the eruption of Crater Peak/Mount Spurr on 19 August 1992 in Alaska.Using multispectral satellite data from the AVHRR/2 and ATSR-2 instruments, Prata and Grant (2001) provided a quantitative analysis of several properties of the Mount Ruapehu ash cloud in New Zealand, including mass loading, cloud height, ash cloud thickness, and particle radius.The quantified ash cloud parameters can be directly inserted into transport and dispersion models as "virtual sources" far from the vent.Wilkins et al. (2014Wilkins et al. ( , 2016) ) applied this technique to the eruption of Eyjafjallajökull in 2010 using infrared (IR) satellite imagery and the NAME model.It was also applied by Crawford et al. (2016) to the 2008 Kasatochi eruption using the HYSPLIT model.
Under a general data assimilation and inverse modeling framework, satellite measurements can be used to constrain the model and estimate emission parameters using various techniques.For instance, Stohl et al. (2011) applied an inversion scheme to the Eyjafjallajökull eruption using a Lagrangian dispersion model with satellite data and demonstrated the effectiveness of the method to yield better quantitative volcanic ash predictions.Schmehl et al. (2012) proposed a variational technique that uses a genetic algorithm (GA) to assimilate satellite data to determine emission rates and steering winds.A HYSPLIT inverse system based on a four-dimensional variational data assimilation approach has been built and successfully applied to estimate the cesium-137 releases from the Fukushima Daiichi Nuclear Power Plant accident in 2011 (Chai et al., 2015).The present work further develops the HYSPLIT inversion system to estimate the time-and height-resolved volcanic ash emission rate by assimilating satellite observations of volcanic ash.The system is tested with the 2008 Kasatochi eruption using the satellite retrievals from passive IR sensors.
The paper is organized as follows.Section 2 describes the satellite observations of volcanic ash, the HYSPLIT model and configuration, and the inverse modeling methodology.Section 3 presents emission inversion results, and Sect. 4 discusses the corresponding volcanic ash forecasts with the estimated source terms.A summary is given in Sect. 5.

Satellite observations
The volcanic ash observations are based on MODIS retrievals from the Terra and Aqua satellites.They include ash mass loading, cloud top height, and effective particle radius.Pavolonis et al. (2013Pavolonis et al. ( , 2015a, b) , b) described the details of the retrieval methodology and how the ash cloud observations are derived from the retrieved parameters, such as radiative temperature and emissivity.Here, volcanic ash observations from the 2008 Kasatochi eruption at five different instances are utilized.The observations were projected onto a latitudelongitude grid with a resolution of 0.05 • latitude and 0.1 • longitude.Figure 1 shows the volcanic ash mass loadings and ash cloud top heights of five granules.Each granule contains 6 minutes of data covering an area of approximately 1500 km along the orbit and 1650 km wide.Note that the satellite observations outside the shown domain are discarded.As the discarded data are mostly located upwind of the volcano vent, they are not expected to provide useful information to estimate the source strength.The places where satellite retrievals did not detect the existence of ash show zero mass loading.It will be shown later that such ash-free regions may be used along with the observed ash cloud to constrain the dispersion model.Note that the ash-free regions do not apply to regions with missing ash mass loadings due to meteorological clouds or other reasons.Table 1 shows the observation time and the number of grid cells with and without    ash detected for each granule.The clear regions dominate the satellite observations.Integrated mass loadings based on the satellite data are also listed in Table 1.They decrease from 9.68 ×10 8 kg for the first granule (G1) to 3.25 ×10 8 kg for the last granule (G5).This probably reflects the gradual loss of the total volcanic ash mass due to deposition.Note that the total mass is likely slightly underestimated for the second granule (G2), where the satellite lost sight of the eastern edge of the ash cloud.

HYSPLIT model configuration
In this study, volcanic ash transport and dispersion are modeled using the HYSPLIT model (Draxler andHess, 1997, 1998;Stein et al., 2015a).A large number of threedimensional Lagrangian particles are released from the source location and passively follow the wind afterward.A random component based on local stability is added to the mean advection velocity in each of the three-dimensional wind component directions to simulate the dispersion.Ash concentrations are computed by summing each particle's mass as it passes over a concentration grid cell and dividing the result by the cell's volume.
Both the NOAA Global Data Assimilation System (GDAS) (Kleist et al., 2009) and the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim global atmospheric reanalysis (Dee et al., 2011) were used as inputs for HYSPLIT.The basic information of the two data sets is listed in Table 2.The concentration grid is set at 0.05 • of resolution for latitude and 0.1 • for longitude with a vertical spacing of 2 km extending from the surface to 20 km.
A total of 290 independent HYSPLIT simulations were run with a unit emission rate released from all possible combinations of 29 different hours from 19:00 UTC on 7 August 2008 to 23:00 UTC on 8 August 2008 and 10 different 2000 m layers.Note that at the first layer, particles are released from the top of the vent from 300 m above sea level to 2000 m.At other layers, particle releases are uniformly distributed throughout the layer at the center of the grid as a line source.In each simulation, particles of four different sizes are released as different pollutants with different fall speeds according to Stokes's law (Heffter and Stunder, 1993).At all release time and height combinations, the contributions to the total mass are assumed to be constant at 0.8, 6.8, 25.4, and 67.0 % for particle sizes of 0.6, 2.0, 6.0, and 20.0 µm, respectively.The same particle size distribution was originally used in the NOAA ARL VAFTAD model (Heffter and Stunder, 1993).Webley et al. (2009) evaluated the sensitivity of the grain size distribution on the modeled ash cloud and found that this predefined distribution is sufficient for HYSPLIT volcanic ash simulation.MODIS effective particle radii (r eff ) are retrieved to describe the ash particle size distributions.However, radii greater than 15-20 µm are not retrieved since the retrievals cannot be performed reliably when r eff exceeds 15 µm (Pavolonis et al., 2013).

Model diagnostic
As shown in Sect.2.1, satellite observations provide ash mass loadings and ash cloud top heights after detecting ash.There are several options to construct the model counterparts for observed ash cloud mass loadings.Three different model diagnostic choices are tested here.In the first option, model volcanic ash concentrations from the ground or sea level up to the model layer where the observed cloud top height resides are integrated to calculate the mass loadings in the model simulation.In the second option, the single model layer where the retrieved cloud top height resides is used to construct the mass loadings.However, the retrieved cloud top heights are associated with uncertainties.Pavolonis et al. (2013) showed that the retrieved cloud top height had a low bias of 0.77 km relative to lidar.Crawford et al. (2016) compared MODIS cloud top height retrievals with CALIOP vertical profiles for the same event.In general, the MODIS top heights agree well with the top aerosol level indicated by CALIOP profiles but can be off by several kilometers.When CALIOP shows two levels of ash, the MODIS top height falls between them.In addition, the cloud top height retrievals typically lie in the middle of thick ash cloud layers rather than at the top (Pavolonis et al., 2013).To compensate for such uncertainties in ash cloud top height position, the third option is designed to integrate model volcanic ash concentrations over three model layers, i.e., from one layer below to one layer above the cloud top layer.
When ash is not detected, the grid cells are flagged as clear or ash-free.This is equivalent to zero mass loading for the entire atmospheric column at such a location.In this case, the model counterpart is obtained by integrating simulated concentrations from the surface to the domain top.Constraining the model simulation with these zero-value observations is expected to help remove spurious sources from which the transport and dispersion will likely generate additional ash clouds that are not observed.
At locations where ash is detected, the observations can be further exploited to provide additional constraints.As ash cloud top heights are provided along with the mass loadings, they indicate that no ash is above the cloud top.How- ever, no information can be inferred for the region below the cloud top.As a result, each ash cell actually generates two pieces of information.Besides the observed volcanic ash cloud mass loadings mentioned earlier, a clear atmospheric column above the cloud top is the other implicit piece of information that can be used in emission inversion.Similar to using zero-value observations at ash-free locations, the integrated mass loadings above the ash cloud top may also be used to filter out unlikely sources.When the "observed" ash cloud is assumed to be limited to one single model layer or three layers, it is also possible to add no-ashbelow-cloud constraints in the inverse modeling.Although such constraints are based on an assumption that is not always true, it will be tested nonetheless.
In addition to detected ash and clear cells, another scenario exists when satellite observations cannot provide positive or negative answers for ash detection, e.g., due to meteorological cloud obstruction.In such a case, no useful information can be used to constrain the model.For the 2008 Kasatochi eruption, overlying meteorological clouds were nearly absent, and valid observations appear across the satellite swaths.

Transfer coefficient matrix (TCM)
A transfer coefficient matrix (TCM) of 290 columns can be generated using all or a subset of the re-gridded MODIS observations listed in Table 1 and the results of the 290 HYS-PLIT simulations with unit emission.A transfer coefficient in the TCM is essentially the mass loadings at an observation point that the row represents, resulting from a dispersion run with a unit emission that the column indicates.
Figure 2 shows the two-dimensional transfer coefficient matrices averaged over all ash pixels for five granules.As a transfer coefficient corresponds to the mass loadings resulting from a unit ash release rate, integrating over more model layers would produce larger transfer coefficients.It is clearly seen that the single layer option, shown as the middle column in Fig. 2, has the averaged TCMs with the lowest values.Figure 2 also shows that integrating from the surface up to the ash cloud top layer generally results in TCMs with the largest values among the three options.As the option to add over three layers (right column in Fig. 2) includes a layer above the cloud top layer that is not included in option one, transfer coefficients at the upper layers may have larger values.Note that a block of zero transfer coefficients after 10:00 UTC on 8 August appear for G1.Ash releases after the observation time of G1 do not affect the G1 observations.In addition, releases need time to travel to the observed location.Figure 2 shows that, as expected, the averaged transfer coefficients tend to be smaller for later observations due to dispersion.The averaged TCMs using ECMWF meteorological data (not shown) are similar to the GDAS results shown here.

Emission inversion
Following a general top-down approach, the unknown emission terms are obtained by searching for the emissions that would provide the model predictions that most closely match the observations.In the current application with the known volcano location, the emission rates vary with time and release heights.With the potential emission time period divided into 29 hourly intervals and the release heights separated into 10 vertical layers, the discretized two-dimensional unknown emission has 290 components to be determined.
Similar to Chai et al. (2015), the unknown releases can be solved by minimizing a cost functional that integrates the differences between the model predictions and the observations, the deviations of the final solution from the first guess (a priori), and other relevant information written into penalty terms (Daley, 1991).For the current application, the cost functional F is defined as where q ij is the discretized two-dimensional emission over M = 29 h and N = 10 layers.q b ij is the first guess or a priori estimate, and σ 2 ij is the corresponding error variance.Note that we assume that the uncertainties of the release at each time-height are independent of each other so that only the diagonal term σ 2 ij of the typical a priori error covariance matrix appears in Eq. (1).For emission points at which the release generates no simulated ash corresponding to any of the assimilated observations, the first guesses remain unchanged.To avoid unrealistic release rates for such emission points, we chose a small constant emission rate

Average transfer coefficient (h m )
-2 Figure 2. Averaged TCMs using GDAS meteorological data with three different options for calculating the model mass loadings (column 1, integrating from the surface to the observed cloud top; column 2, calculated for a single layer where the observed cloud top height resides; column 3, integrating over three layers centered at the observed cloud top layer).Rows 1-5 (from top to bottom) correspond to observations G1-G5.
of 10 4 g h −1 (≈ 2.8 × 10 −3 kg s −1 ) at all hours and layers as the first guess.Large uncertainties are given in the following tests to reflect the fact that little was known for the mass emission rates.a h m and a o m are the mass loadings simulated by HYSPLIT and retrieved by MODIS, respectively.The observations here refer to both the volcanic ash mass loadings for the ash cloud and the zero values for the ashfree regions, which are later included as extra constraints in Sect.4.3.Zero mass loadings also include those calculated over the atmospheric columns above or below ash clouds as discussed earlier in Sect.2.4. 2 m includes the variances of the observational and representative errors.For simplicity, 2 m are referred to as observational errors hereafter and are assumed to be uncorrelated.Dubuisson et al. (2014) studied the remote sensing of volcanic ash plumes from SEVIRI, MODIS, and IASI instruments.The total uncertainty in MODIS mass loading resulted from errors in the input atmospheric parameters, such as ash layer altitude, particle size distribution, and particle composition, and was estimated to be ∼ 50 %.Their inter-comparison among six satellite configurations shows a standard deviation of 0.3 g m −2 for the mean mass loading estimates.In this study, the observational errors are estimated using m = 0.50 × a o m + 0.3 g m −2 .No smoothness penalty term is included in the cost functional because of the abrupt nature of the volcanic eruptions.A large-scale boundconstrained limited-memory quasi-Newton code, L-BFGS-B (Zhu et al., 1997), is used to minimize the cost functional F defined in Eq. ( 1).The maximum number of cost functional evaluations is set as 250 for cases in Sect. 3 and 2500 for those in Sect. 4. To ensure non-negative q ij solutions from the optimization, q ij is converted to ln(q ij ) as input to the L-BFGS-B routine.An alternative to this is enforcing the q ij ≥ 0 with lower bounds enabled by the L-BFGS-B routine.As they solve the same mathematical problem, these two options are expected to arrive at the same results with enough iterations.Chai et al. (2015) provide a detailed discussion on the conversion of control and metric variables.Although they showed that using logarithmic concentration differences in the cost functional performed better than directly using concentration differences in their application, the logarithmic conversion on the metric variable a m is not beneficial for the current application.This is because the range of the volcanic ash mass loadings here is much smaller than that of the cesium-137 air concentrations encountered in their application.In addition, the utilization of zero mass loadings in many ash-free regions prohibits using ln(a o m ).In this study, the mass loadings are directly compared in the cost functional without logarithmic conversion.

Emission estimates
The emission estimates obtained by minimizing the cost functional F introduced in Eq. ( 1) highly depend on the uncertainties given to the a priori and observations.Sensitivity tests are first performed by changing the magnitudes of the a priori error variances, while the observational error estimation is fixed.Chai et al. (2015) demonstrated that the emission inversion results were not sensitive to the first guess of the emissions when large uncertainties are presumed.
In the sensitivity tests, ash cloud data at G1 and G2 are assimilated.Note that the zero mass loading values for ashfree regions are not used here.Large a priori error variances are presumed with σ ij ≈ 10 12 g h −1 (≈ 2.8 × 10 5 kg s −1 ) and σ ij ≈ 10 16 g h −1 (≈ 2.8 × 10 9 kg s −1 ).In these cases, the HYSPLIT simulated mass loadings were calculated by integrating from the surface to the observed ash cloud top heights at the ash cells.Figure 3 shows that the emission inversion results are slightly different from each other when the a priori errors are assumed differently, as expected.However, similar patterns are apparent for both cases with the different a priori error variances.A peak release greater than 5000 kg s −1 is observed at 04:00 UTC on 8 August 2008 at the 6-8 km layer for both cases.This demonstrates that the emission estimates are most decided by the satellite data when a priori errors are assumed to be large enough.Note that a larger a priori term with smaller a priori error variances in Eq. ( 1) typically helps the minimization procedure in emission inversion.Since the results using the two a priori errors are similar, the a priori error variances are set as σ ij ≈ 10 12 g h −1 (≈ 2.8 × 10 5 kg s −1 ) in the following tests.Waythomas et al. (2010) characterize the eruption by three major explosive events and two smaller events.Events one and two started at 22:01 UTC on 7 August and 01:50 UTC on 8 August, respectively.These two events reached 14 km and produced water-rich but ash-poor clouds.Event three happened at 04:35 UTC on 8 August.It generated an ashrich cloud that rose up to 18 km.About 16 h of continuous ash emission was punctuated with events four and five at 07:12 UTC and 11:42 UTC on 8 August.
Figure 4 shows the emission estimates using all three options in calculating model mass loadings.The zero values for ash-free regions are not used here.The emission results are significantly different with different options.For cases in which the model counterparts of the satellite mass loadings are obtained by integrating from the surface to the cloud top, the ash releases started at 01:00 UTC on 8 August 2008 from the 8-10 km layer.The emissions lasted for 4 h and extended to multiple layers, reaching up to the 14-16 km layer and down to the 4-6 km layer.After 1 h without ash, moderate volcanic ash releases continued for 6 h until 12:00 UTC on 8 August, mainly between 8 and 16 km.A small ash emission of less than 80 kg s −1 is seen at the 12-14 km layer starting at 15:00 UTC for 1 h.If the model mass loadings are obtained by only considering a single layer where the cloud top height resides, the resulting release source terms are limited to layers between 12 and 16 km.The ash releases started at 03:00 UTC on 8 August and lasted for 3 h before resuming again 2 h later.With emissions on and off for the next 2 h Figure 3. Volcanic ash release results with different a priori error estimations (top, σ ij ≈ 2.8 × 10 5 kg s −1 ; bottom, σ ij ≈ 2.8 × 10 9 kg s −1 ).The TCMs for the emission inverse were generated using HYSPLIT runs with GDAS meteorological data.Only the ash cells from the satellite data at G1 and G2 are used in the emission inverse.The model counterparts are obtained by integrating from the surface to the ash cloud top heights at ash cells.
at the 14-16 km layer, the ash release continued for 6 h and peaked between 14:00 and 15:00 UTC on 8 August at the 12-14 km layer.There is also an isolated emission point at the 14-16 km layer starting at 23:00 UTC on 8 August and lasting for 1 h.In the last case, for which the model mass loadings are calculated by integrating over three layers centered at the cloud top layer, the ash releases are drastically different from the first two cases.The ash releases start much earlier, at 20:00 UTC on 7 August, and the release heights are within the 14-18 km range.The release then extended to more layers, but the main sources went lower.This lasted for 13 h before stopping at 09:00 UTC on 8 August.A second spurt of ash release started at 11:00 UTC from the 14-16 km layer and remained above 12 km before pausing again 5 h later.Several weaker ash releases are found between the 14 and 18 km layers at later times, from 19:00 UTC on 8 August to 00:00 UTC on 9 August.The three emission estimates in Fig. 4 do not reproduce the eruption as described by Waythomas et al. (2010), but they manage to capture some characteristics of the eruption.Without information on the vertical profiles of the ash cloud, how the mass loadings are interpreted greatly affects the release estimates, as shown by the drastic differences between the estimates in Fig. 4. Thus, it is difficult to generate reliable and accurate actual volcanic ash emission estimates if the ash cloud vertical structures are undetermined.However, it will be shown later that such emission estimates can still help to improve ash cloud forecasts.

Ash predictions with top-down emission estimates
A series of tests was performed to find the best inverse modeling setup.In Sect.4.1, the evaluation metrics are described.In Sect.4.2, the choices for calculating the model counterparts of the satellite mass loadings are compared.In Sect.4.3, whether to use the ash-free region to constrain the model is investigated.In Sect.4.4, the effect of keeping older observations when newer observations become available is discussed.

Evaluation metrics
For the model evaluation, total column mass loadings are constructed by integrating predicted concentrations from the surface to the domain top.They are used to compare the satellite observations in each granule shown in Fig. 1, including both ash and clear points.The aim is to use the total column mass loadings instead of any of the options described in Sect.2.4 to provide a fair comparison among the three options by avoiding the complexities associated with the vertical structures of the volcanic ash cloud.Note that Crawford et al. (2016) excluded masses below 2 km when integrating the model results to obtain the mass loadings because the satellite retrieval is less sensitive to low-level ash.Such exclusion may improve the evaluation statistics, but it is not expected to affect the inter-comparison between different model runs.Mean bias (MB), fractional bias (FB), root mean square error (RMSE), normalized RMSE (NRMSE), and the Pearson correlation coefficient (R) are calculated.FB and NRMSE are scaled by the average of model and observation means.In addition, the critical success index (CSI) defined below is calculated for ash detection: A threshold of 0.1 g m −2 , the approximate lower limit of the MODIS satellite data set, is used to categorize ash existence for both the model predictions and the observations.N Hit , N FalseAlarm , and N Miss denote the numbers of grid points where ash is predicted and observed, where ash is predicted but not observed, and where ash is observed but not predicted by the model, respectively.
Following Draxler (2006), the Kolmogorov-Smirnov parameter (KSP) and the "rank" are calculated.The KSP measures the largest difference between the cumulative distribution functions of the model predictions and the observations.As shown in Eq. ( 3), the rank adds up four components that all range from 0 to 1.The larger rank values indicate a better overall performance of the model results: (3)

Model mass loadings
The HYSPLIT predictions using the estimated source terms after assimilating the G1 and G2 observations are evaluated against the satellite observations of G2, G3, G4, and G5.Note that the zero mass loadings for ash-free regions are not used here.The three options to calculate the model ash mass loadings discussed earlier are employed in the inverse modeling.The statistics are listed in Table 3.
Comparing against the G2 observation, Table 3 shows that integrating over three model layers yields (option M1) slightly better results based on most statistics.This is true for cases with both GDAS and ECMWF meteorological fields.The advantage of the M1 option is not apparent when comparing against other observations.Based on rank, the ECMWF cases are better than the GDAS cases against G2, but the ranks for the ECMWF cases deteriorate faster with time and become worse than the GDAS cases when the model output is compared to the G4 and G5 observations.The model predictions have the best statistics compared against G4 than against the other satellite granules (G2, G3, and G5).The case with GDAS meteorological fields and the three-layer mass loading option M1 has the best rank of 3.02 (FB = 0.04, R = 0.72, CSI = 0.62, KSP = 0.10).If only G2 observations were assimilated, the model performance would be expected to peak when compared against G2.However, as both the G1 and G2 observations are assimilated, this is no longer true.The effect of assimilating different observations will be discussed later in Sect.4.4.Table 3 shows that the model tends to underestimate the ash mass loadings of G2 and G3 and then mostly overestimate the ash mass loadings of G4 and G5.It results in the best FB against G4 for GDAS cases and the best FB against G3 for ECMWF cases as the FB signs change.Since the volcanic ash will disperse with time, the average mass loadings get smaller.This is reflected in a basic trend of decreasing RMSEs with time, although the NRMSEs slightly increase.
While different evaluation metrics may not always agree with each other, the overall performance parameter rank provides a simplified way to compare the model results.Only the ranks are listed and used to compare the model predictions hereafter.Using HYSPLIT ensembles, Stein et al. (2015b) estimated the uncertainties of the rank as 0.08, 0.08, 0.09, 0.08, 0.11, and 0.07 for six different tracer releases.The uncertainties of the rank for the current application could vary, but they are not expected to be too different.

Extra constraints
As discussed in Sect.2.3, ash-free regions indicate zero mass loadings for the entire atmospheric columns.Cloud top heights can also be used to enforce ash-free atmospheric columns above volcanic ash clouds.In addition, ash-free at-Table 3. Evaluation statistics against the G2, G3, G4, and G5 observations for cases with different ways to calculate model mass loadings.G1 and G2 are assimilated for all cases listed here.MET is meteorological input; OBS indicates satellite observations used for evaluation; ML is mass loading; MA is integrating from the surface to the cloud top; M0 is calculated for a single layer where the cloud top height resides; M1 is integrating over three layers centered at the cloud top layer; MB is the mean bias; FB is the fractional bias; RMSE is the root mean square error; NRMSE is normalized RMSE; R is the Pearson correlation coefficient; CSI is the critical success index; and KSP is the Kolmogorov-Smirnov parameter.Rank is defined in Eq. (3).mospheric columns below the ash cloud may be assumed if an ash cloud thickness is estimated.Note that the term "above or below ash cloud" is in relation to the chosen model's cloud diagnostic.For instance, if the M1 option is chosen, aboveand below-ash-cloud constraints are enforced over the model layers outside the three ash layers.Whether such extra constraints benefit the inverse modeling is tested here using the 22 inverse cases listed in Table 4.The ranks evaluated against G2-G5 are listed.It is found that when the additional constraints of including the clear pixels outside the ash cloud are used, the ranks decrease.This holds true against G2-G4 for all three mass loading calculation options and for both sets of meteorological data.Two exceptions are found against G5 for the ECMWF cases with the M0 and M1 options, in which ranks increase from 2.17 to 2.32 and 2.28 to 2.38, respectively.Enforcing the extra constraints of the ash-free regions makes the inversion results very sensitive to transport errors since the HYSPLIT simulated ash plume outside the MODIS ash cloud starts to affect the emission inversion results.Table 4 shows that the emission inversion with extra constraints of clear pixels using ECMWF data performs better than using GDAS data, except in a single case with the MA option against G4.

MET
Adding the extra constraints of a clear column above the ash cloud generally causes a decrease in rank.An exception is the ECMWF case with the M1 option (three model layers used for mass loading calculation) in which the extra "top" constraint results in marginally better predictions evaluated against G5 (a rank of 2.39 versus 2.38).It is found that the ECMWF cases perform better than all of their GDAS counterparts after adding the top constraints.When the constraints of a clear column below the ash cloud are further added for the M0 and M1 options, the ranks decrease significantly, especially for the M0 option in which a single model layer is used to construct the model mass loadings.Clearly, model and observation uncertainties have to be carefully addressed Table 4. Ranks of the inverse tests with various extra constraints evaluated against the G2, G3, G4, and G5 observations (OBS).For mass loading (ML), MA is integrating from the surface to the cloud top; M0 is calculated for a single layer where the cloud top height resides; and M1 is integrating over three layers centered at the cloud top layer.For extra zero observation constraints, H is with clear pixels; T is with a clear column above the ash cloud; and B is with a clear column below the ash cloud.Ash cells are assimilated in all inverse cases.The satellite data at both G1 and G2 are used for all cases listed here.to take advantage of the extra constraints in order to benefit the emission inversion.This requires further investigation in future studies.

Older observations
As newer observations become available, whether to include the older observations in the assimilation remains a question.Table 5 lists the statistics for 10 cases evaluated against granules 2-5 using both GDAS and ECMWF fields.In the inverse modeling, only ash pixels were used, and the model mass loadings are calculated by integrating over three layers centered at the cloud top layer (M1 option).It is found that assimilating G2 and G1 yields greater ranks when comparing against G3 and G4 observations than assimilating G2 alone.At G5, there is little difference between the two strategies.Note that assimilating G2 alone helps to generate better statistics against the same observations than assimilating G1 and G2 at the same time, although this does not help the G3 forecasts about 12 h later.After G3 is available, three strategies to utilize the available observations G1, G2, and G3 are tested.The results show that assimilating G2 along with G3 observations achieves better forecasts at G4 and G5 moments than assimilating only G3.It is also found that including G1 in the assimilation does not make much difference.Again, the assimilation of G3 alone results in a closer match between the model predictions and the G3 observations, but the forecasts at later times are worse than when the earlier observations are also assimilated.
Figure 5 shows the comparison between the MODIS observations and the HYSPLIT simulations using the estimated source terms obtained by assimilating G1, G2, and G3 with both GDAS and ECMWF meteorological fields, listed as the last two cases in Table 5.The simulated ash clouds corresponding to G1 are narrower than the satellite observations, and the mass loading values are underestimated.Crawford et al. (2016) found that cylindrical source terms performed better than the line sources assumed here.Waythomas et al. (2010) showed that the source area was quite broad with a width of about 75 km from 06:00 UTC to 10:00 UTC on 8 August 2008.Inverse modeling with cylindrical sources will be investigated in the future.The HYSPLIT simulations with both meteorological fields agree well with granules G2 and G3, and this is reflected by the high rank vales (Table 5).This is expected, as the same observations were assimilated to obtain the ash release rates.Against G4, the model results capture the ash cloud locations and magnitudes very well for both cases.The case with GDAS input appears to have similar mass loading values as the observations, while the ECMWF case has a narrow ring inside the main ash cloud with higher values than the MODIS observations.In addition, the ECMWF case shows two tails, while the GDAS case has only one tail resembling the MODIS observations.Both cases show tapering shapes of the tails, which appear different from the satellite view.Against the later observations of G5, the HYSPLIT simulations start to deviate from the MODIS, as indicated by the lower rank.Both the GDAS and ECMWF simulations capture the ash cloud at similar locations as observed by the satellite but show smoother structures.It is speculated that meteorological fields with higher spatial and temporal resolutions might be able to improve the ash cloud predictions.  1 for detail)."+" shows the location of the Kasatochi volcano (52.1714 • N, 175.5183 • W).The white areas indicate regions outside the satellite granules for MODIS observations.For HYSPLIT simulations, the white areas indicate zero mass loadings in order to reveal the ash cloud boundaries.The ash release rates for the HYSPLIT simulations were obtained by assimilating granules G1, G2, and G3.In the inverse modeling, only ash pixels were used, and the model mass loadings are calculated by integrating over three layers centered at the cloud top layer.
Table 5. Ranks against G2-G5 for HYSPLIT simulations after assimilating various combinations of observation inputs.The model counterparts of the satellite mass loadings are calculated using the M1 option, i.e., integrating over three layers centered at the cloud top layer.Only the ash cells are assimilated for all the inverse cases listed here."()" indicates that the observations have been assimilated.There were several lidar observations of the Kasatochi ash cloud provided by the CALIPSO satellite (Winker et al., 2010;Kristiansen et al., 2010;Crawford et al., 2016).The HYSPLIT simulations shown in Fig. 5 are also compared against the 532 nm backscatter vertical profiles along the three CALIPSO overpasses coincident with G1, G4, and G5.The comparisons reveal that both the GDAS and ECMWF simulations captured the main ash cloud features at approximately the same location and altitude.

Summary
An inverse system based on HYSPLIT has been developed to solve the effective volcanic ash release rates as a function of time and height by assimilating satellite mass loadings and ash cloud top heights.The Kasatochi eruption in 2008 was used as an example to test and evaluate the current top-down system with both GDAS and ECMWF meteorological fields.
When quantifying the differences between the model predictions and the satellite observations, the model counterparts can be calculated differently using the threedimensional model concentration results because the observed ash cloud bases are unknown.Three options to construct the model mass loadings are tested for this inverse system, integrating volcanic ash concentrations from the surface up to the cloud top height or only using one or three model layers.It is found that the emission estimates vary significantly with different options.However, all of the predictions with the different estimated release rates show decent skill when evaluated against the unassimilated satellite observations at later times.The option of integrating over three model layers yields slightly better results than integrating from the surface up to the cloud top or using a single model layer.
The extra constraints of enforcing zero mass loading in the ash-free regions are tested with the inverse system.The model predictions using the emission estimates generated with such extra constraints are worse than those using the emission estimates generated by only assimilating the ash pixels.Additional "no-ash" constraints for the atmosphere columns above or below the observed ash cloud top height are found to further deteriorate the subsequent model predictions using the top-down emission estimates.
Assimilating multiple granules at different times proves to be beneficial.As new observations become available, the effect of 1-day old observations becomes marginal, but assimilating mass loadings from the most recent observations and those from about 12 h earlier yields better results than only assimilating the most recent observations.
The spatial and temporal resolutions of the meteorological fields may need improvement for future studies.The line source assumed here can be replaced by more realistic cylindrical sources in the future.A simple particle size distribution with four different particle sizes is used at all release heights and times.With the MODIS effective radius available, a more realistic way to represent the particle size distribution can be explored.

Figure 1 .
Figure1.MODIS volcanic ash mass loadings (left) and ash cloud top height (right) listed from top to bottom by observation time (see Table1for detail)."+" shows the location of the Kasatochi volcano (52.1714 • N, 175.5183 • W).Note that the satellite observations to the left of the map domain are not used in this paper.

Figure 4 .
Figure 4. Volcanic ash release estimates with different options in the model mass loading calculation: (top) integrating from the surface to the cloud top (same as Fig. 3, top), (center) calculated for a single layer where the cloud top height resides, and (bottom) integrating over three layers centered at the cloud top layer.

Figure 5 .
Figure5.Volcanic ash mass loadings from MODIS (left) and HYSPLIT simulations with GDAS (center) and ECMWF (right) from top to bottom by observation time (see Table1for detail)."+" shows the location of the Kasatochi volcano (52.1714 • N, 175.5183 • W).The white areas indicate regions outside the satellite granules for MODIS observations.For HYSPLIT simulations, the white areas indicate zero mass loadings in order to reveal the ash cloud boundaries.The ash release rates for the HYSPLIT simulations were obtained by assimilating granules G1, G2, and G3.In the inverse modeling, only ash pixels were used, and the model mass loadings are calculated by integrating over three layers centered at the cloud top layer.

Table 1 for
detail)."+" shows the location of the Kasatochi volcano (52.1714 • N, 175.5183 • W).Note that the satellite observations to the left of the map domain are not used in this paper.

Table 1 .
Description of MODIS ash cloud observations."Ash cells" and "clear cells" show the number of grid cells with and without detected ash, respectively.The total mass is obtained by integrating mass loadings over the observed region.

Table 2 .
Description of GDAS and ECMWF meteorological data.