Accelerating volcanic ash data assimilation using a mask-state algorithm based on an ensemble Kalman filter : a case study with the LOTOS-EUROS model ( version 1 . 10 )

In this study, we investigate a strategy to accelerate the data assimilation (DA) algorithm. Based on evaluations of the computational time, the analysis step of the assimilation turns out to be the most expensive part. After a study of the characteristics of the ensemble ash state, we propose a mask-state algorithm which records the sparsity information of the full ensemble state matrix and transforms the full matrix into a relatively small one. This will reduce the computational cost in the analysis step. Experimental results show the mask-state algorithm significantly speeds up the analysis step. Subsequently, the total amount of computing time for volcanic ash DA is reduced to an acceptable level. The mask-state algorithm is generic and thus can be embedded in any ensemble-based DA framework. Moreover, ensemble-based DA with the mask-state algorithm is promising and flexible, because it implements exactly the standard DA without any approximation and it realizes the satisfying performance without any change in the full model.


Introduction
Volcanic ash erupted into atmospheres can lead to severe influences on aviation society (Gudmundsson et al., 2012). Turbine engines of airplanes are extremely threatened by ash ingestion (Casadevall, 1994). Thus, accurate real-time aviation advice is highly required during an explosive volcanic ash eruption (Eliasson et al., 2011). Using data assimilation (DA) to improve model forecast accuracy is a powerful approach (Lu et al., 2016a). Recently, ensemble-based DA (Evensen, 2003) has been evaluated as very useful for improving volcanic ash forecasts and regional aviation advice (Fu et al., 2016). It corrects volcanic ash concentrations by continuously assimilating observations. In Fu et al. (2016), real aircraft in situ measurements were assimilated using the ensemble Kalman filter (EnKF), which is the best known and most popular ensemble-based DA method. Based on the validation with independent data, ensemble-based DA was concluded as being powerful for improving the forecast accuracy.
However, to make the methodology efficient also in an operational (real-time) sense, the computational efforts must be acceptable. For volcanic ash DA problems, so far, no studies on the computational aspects have been reported in the literature. Actually, when large amounts of volcanic ash erupted into atmospheres, the computational speed of volcanic ash forecasts is just as important as the forecast accuracy (Zehner, 2010). For example, due to the lack of a fast and accurate forecast system, the sudden eruption of the Eyjafjallajökull volcano in Iceland from 14 April to 23 May 2010 caused an unprecedented closure of the European and North Atlantic airspace, resulting in a huge global economic loss of USD 5 billion (Oxford-Economics, 2010). Since then, research on fast and accurate volcanic ash fore-Published by Copernicus Publications on behalf of the European Geosciences Union.
casts has gained much attention, because it is needed to provide timely and accurate aviation advice for frequently operated commercial airplanes. It was shown that the accuracy of volcanic ash transport can be significantly improved by the DA system in Fu et al. (2016). Therefore, it is urgent to also consider the computational aspect, i.e., improving the computational speed of the volcanic ash DA system as quickly as possible. This is the main focus of this study.
Due to the computational complexity of ensemble-based algorithms and the large scale of dynamical applications, applying these methods usually introduces a large computational cost. This has been reported from the literature on different applications. For example, for operational weather forecasting with ensemble-based DA, Houtekamer et al. (2014) reported computational challenges at the Canadian Meteorological Center with an operational EnKF featuring 192 ensemble members, using a large 600 × 300 global horizontal grid and 74 vertical levels. An initialization requirement of over 7 × 10 10 values to specify each ensemble results in large computational efforts on the initialization and forecast steps in weather forecasting. For oil reservoir historymatching (Tavakoli et al., 2013), the reservoir simulation model usually has a large number of state variables; thus, the forecasts of an ensemble of simulation models are often time-consuming. Besides, when time-lapse seismic or dense reservoir data are available, the analysis step of assimilating these large observations becomes very time-consuming (Khairullah et al., 2013). Large computational requirements of ensemble-based DA have also been reported in ocean circulation models (Keppenne, 2000;Keppenne and Rienecker, 2002), tropospheric chemistry assimilation (Miyazaki et al., 2015), and many other applications.
To accelerate an ensemble-based DA system, the ensemble forecast step can first be parallelized because the propagation of different ensemble members is independent. Thus if a computer with a sufficiently large number of parallel processors is available, all the ensemble members can be simultaneously integrated. In the analysis stage, to calculate the Kalman gain and the ensemble error covariance matrix, all ensemble states must be combined together. In weather forecasting and oceanography sciences, Keppenne (2000), Keppenne and Rienecker (2002), and Houtekamer and Mitchell (2001) have reported using parallelization approaches to accelerate the expensive analysis stage. In reservoir history matching, a three-level parallelization has been proposed by Tavakoli et al. (2013); Khairullah et al. (2013) in recent years, to significantly reduce computational efforts of both forecast and analysis steps due to massive dense observations and large simulation models. The first parallelization level is to separately perform the ensemble simulations on different processors during the forecast step. This approach is usually quite efficient when a large ensemble size is used. However, the scale or model size of one reservoir simulation is constrained by the memory of a single processor. Thus, the second parallelization level is to perform one ensemble member simulation using a parallel reservoir model. These two levels do not deal with the analysis step, which collects all ensemble members to do computations usually on a single processor. Therefore, a third level of parallelization was implemented by Tavakoli et al. (2013) and Khairullah et al. (2013) by parallelizing matrix-vector multiplications in the analysis steps. Furthermore, some other approaches on accelerating ensemble-based DA systems have also been reported, such as GPU-based acceleration (Quinn and Abarbanel, 2011) in numerical weather prediction (NWP) and domain decomposition in atmospheric chemistry assimilation (Segers, 2002;Miyazaki et al., 2015). The observations used in an DA system can also be optimized with some preprocessing procedures, as reported by Houtekamer et al. (2014).
Although for other applications there were many efforts in dealing with large computational requirements in an ensemble-based DA system, most of them cannot be directly used to accelerate volcanic ash DA. This is because the acceleration algorithms are strongly dependent on specific problems, such as model complexity (high or low resolution), observation type (dense or sparse), or primary requirement (accuracy or speed). These factors determine, for a specific application, which part is the most time-consuming, and which part is intrinsically sequential. Thus, no unified approach for efficient acceleration of all the applications can be found. Although the successful approaches in other applications cannot be directly employed in volcanic ash forecasts, their success does stress the importance of designing a proper approach based on the computational analysis of a specific DA system. Therefore, the computational cost of our volcanic ash DA system will first be analyzed. Then, based on the computational analysis, we will investigate a strategy to accelerate the ensemble-based DA system for volcanic ash forecasts. This paper is organized as follows. Section 2 introduces the methodology of volcanic ash DA. Section 3 analyzes the computational cost of the conventional volcanic ash DA system. In Sect. 4, the mask-state algorithm (MS) is developed for acceleration. The comparison between MS and standard sparse matrix methods is presented in Sect. 5. The discussions on MS is in Sect. 6. Finally, the last section summarizes the concluding remarks of our research.
To simulate a volcanic ash plume, an atmospheric transport model is needed. In this paper, the LOTOS-EUROS (abbreviation of LOng Term Ozone Simulation -EURopean Operational Smog) model is used (Schaap et al., 2008) with model version 1.10 (http://www.lotos-euros.nl/). The LOTOS-EUROS model (Schaap et al., 2008) is an operational model focusing on nitrogen oxides, ozone, particulate matter, and volcanic ash. The model configurations for volcanic ash were discussed in detail by Fu et al. (2016). For volcanic ash simulation, the model is configured with a state vector of size 180 × 200 × 18 × 6 (the dimensions correspond to longitude, latitude, vertical level, and ash species), and the size of the model state is thus calculated as ∼ 10 6 .
The experiment in this study starts at t 0 (09:00 UTC, 18 May 2010 for this study) by considering an initial condition from a previous LOTOS-EUROS conventional model run (see Fig. 1a). In the second step (the forecast step) the model propagates the ensemble members from the time t k−1 to t k (k > 0, the time step is 10 min): The operator M k−1 describes the time evolution of the state which contains the ash concentrations in all model grid boxes. The state at the time t k has a distribution with the mean x f and the forecast error covariance matrix P f given by where L f represents the ensemble perturbation matrix. In this study, the forecast step is performed in parallel because of the natural/common parallelism of the independent ensemble propagation, which is a trivial approach when employing ensemble-based DA (Liang et al., 2009;Tavakoli et al., 2013;Khairullah et al., 2013). When the model propagates to 09:40 UTC, 18 May 2010, the volcanic ash state gets sequentially analyzed by the DA process by combining real aircraft in situ measurements of PM 10 and PM 2.5 concentrations until 11:10 UTC. The measurement route and values are demonstrated in Fig. 1b, c and the details are described in Weber et al. (2012) and Fu et al. (2016). The observational network at time t k is defined by the operator H k which maps the state vector x to the observational vector y by where y contains the aircraft measurements and v represents the observational error. H k selects the grid cell in x(k) that corresponds to the locations of the observation. When measurements are available, the ensemble members are updated in the analysis step using where K represents the Kalman gain, H is the observational matrix formed by the observational operator H , R represents the measurement error covariance matrix, and v j represents the realization out of the observation error distribution v. After the continuous assimilation ending at 11:10 UTC, the forecast at 12:00 UTC is illustrated in Fig. 1d, for which the forecast accuracy has been carefully evaluated as significantly improved compared to the case without DA (Fu et al., 2016). The EnKF with the above setups is abbreviated as "conventional EnKF" and used in this study for the computational evaluation. Note that in the study we do not use covariance localization as proposed by Hamill et al. (2001) for reducing spurious covariances. This is because although localization is possible, the ideal case is not to use it in order to have the correct covariances in a large (converged) ensemble. It is crucial for localization that when unphysical (spurious) covariances are eliminated, physical (correct) covariances can be well maintained (Petrie and Dance, 2010). If the "filtering length scale" for localization is too long (i.e., all the dynamical covariances are allowed), many of the spurious covariances may not be eliminated. If the length is too short, important physical dynamical covariances then may be lost together with the spurious ones. Therefore, essentially deciding on an accurate localization is a challenging subject (Riishojgaard, 1998;Kalnay et al., 2012), especially for accuracydemanding applications. Therefore, in this study we choose the ensemble size of 100 to guarantee the accuracy and avoid large spurious covariances.
3 Computational analysis for volcanic ash DA

Computational analysis of the total runtime
Ensemble-based DA is a useful approach to improve the forecast accuracy of volcanic ash transport. However, if it is timeconsuming, it cannot be taken as efficient due to the high requirement on speed for volcanic ash DA (see Sect. 1). Based on this consideration, we need to analyze the computational cost of a conventional volcanic ash DA system.
As introduced in Sect. 2, the total execution time of conventional EnKF comprises four parts, i.e., initialization, forecast, analysis, and other computational cost. The initialization time includes reading meteorological data, initializing model geographical and grid configurations, reading emission information, initializing stochastic observers for reading and transforming observations to the model grid, and initializing all the ensemble states and ensemble means. The fore- cast time is obtained from Eq. (1), while the analysis time corresponds to the computational sum from Eqs. (2) to (7). The other computational time includes script compiling, setting environment variables, and starting and finalizing DA algorithms.
The evaluation result of the conventional EnKF is shown in Table 1 (the middle column). It can be seen that the total computational time (4.36 h) is relatively large compared to the simulation window (3.0 h, i.e., from 09:00 to 12:00 UTC, 18 May 2010), which is too much in an operational sense. Therefore, in this study, we aim to accelerate the computation to within an acceptable runtime (i.e., requiring less runtime than the time period of the DA application).
It can also be observed from Table 1 that the main contribution to the total execution time is the analysis step. Compared to the initialization and forecast time, the analysis stage takes 72 % of the total runtime. Due to the expensive analysis step, although some approaches (such as MPI-parallel I/O Filgueira et al., 2014, domain decomposition Segers, 2002 can potentially accelerate the initialization and forecast step, the effect on the final acceleration of the total computational cost is little. Therefore, to get an acceptable computational time, the cost reduction in the analysis step is the target. One may wonder that since the number of observations is small, why does analysis take so much time? The large state vector seems to be left responsible for the problem. To know the exact reason, the detailed computational cost of the analysis step must be evaluated.

Cost estimation of all analysis procedures
We start with the formulations of the analysis step. The analysis step is represented by Eq. (7), which can be written in a full matrix format with Eq. (8), where the subscripts represent the matrix's dimensions. A f and A a represent the forecasted and analyzed ensemble state matrix, and are respectively built up from ξ f and ξ a with N ensembles. The measurement ensemble matrix Y is formed by an ensemble of y + v (see Eq. 7). H is the observational matrix, which is used to select state variables (at measurement locations) in the full ensemble state matrix corresponding to the measurement ensemble matrix Y. n is the number of model state variables in a three-dimensional (3-D) domain, i.e., ∼ 10 6 in this study (see Sect. 2). m is the number of measurements at one assimilation time, which depends on the measurement type. For aircraft in situ measurements used in this study (see Fig. 1c), two measurements are made at each time by one research flight, so that m is 2 here. N is the ensemble size and is taken as 100 in this study. As described in Eq.
(3), the ensemble perturbation matrix L f in EnKF can be re-written as where I is an N × N unit matrix and 1 is an N × N matrix with all elements equal to 1.
Here we explicitly express L f and HL f in the form of A f and O f , respectively. This is because in our volcanic ash DA system, A f and O f are two of the three inputs (another one is the measurement ensemble matrix Y for the analysis step). These are the three inputs used for actual computations in the analysis step. As shown in Fig. 2a, A f is obtained from the forecast step, and O f and Y are acquired from our stochastic observer module (see Fig. 2a) which is used for a volcanic ash transport model to integrate geophysical measurements. With the input Y, the measurement error covariance R, as introduced in Eq. (6), can then be computed with Based on previous definitions and Eqs.
(2) to (7), the analysis step can be reformulated as follows: Equation (11) shows how the analysis step is performed in a volcanic ash DA system. In order to accelerate the analysis step, the most time-consuming part must be reduced. each procedure in the analysis step. Considering that the state number n (∼ 10 6 ) is significantly larger than the measurement number m (m = 2 here) and the ensemble size N (N = 100), the most time-consuming procedure in the analysis step is thus the last one, that is A a = A f X with a computational cost of O(n N 2 ). Therefore, in our volcanic ash DA system, this part is the most time-consuming part in the analysis step.
for singular value decomposition (SVD) in our study is not timeconsuming, which is different from applications of reservoir history matching (Tavakoli et al., 2013;Khairullah et al., 2013). This is because of the SVD procedure costs O(m 3 ), and due to the measurement size on the order of the size of the state in those cases, the SVD procedure thus requires a huge computational cost for reservoir DA. Analysis in the previous section shows that A a = A f X is most expensive in the analysis step. Each column of A f is constructed from a forecasted ensemble state; thus, the dimension of A f is n × N. In each column, the element values correspond to volcanic ash concentrations in a 3-D domain. Figure 3 shows the coverage of all ensemble forecast states at a selected time, 10:00 UTC, 18 May 2010, without loss of generality. A common phenomenon can be observed: that is, only a part of the 3-D domain is filled with volcanic ash. The ash clouds only concentrate in a plume which is transported over time. This is because volcanic eruption is a fast and strong process. The advection dominates the transport, and the volcanic ash plume is transported with the wind. This is a particular characteristic of volcanic ash transport, Geosci. Model Dev., 10, 1751-1766, 2017 www.geosci-model-dev.net/10/1751/2017/ in contrast to other atmospheric-related applications such as ozone (Curier et al., 2012), SO 2 (Barbu et al., 2009), and CO 2 (Chatterjee et al., 2012). For those applications, the concentrations are everywhere in the domain, the emission sources are also everywhere, and observations are available throughout the domain too (especially for satellite data), whereas for application of volcanic ash transport, the source emission is only at the volcano; thus, usually only a limited domain is polluted by ash. As shown in Fig. 3, in the 3-D domain with a grid size of 3.888 × 10 6 , the number of grids in the area with volcanic ash is counted as 1.528 × 10 6 , whereas the number of no-ash grids is 2.36 × 10 6 . Note that shown in the figure are accumulated ash coverages of all ensemble states; thus, in the no-ash grids, there is no ash for any of the ensemble states. Thus a very large number of rows in A f are zero corresponding to the no-ash grids. These zero rows in A f have no contributions to A a = A f X, because a zero row in A f always results in a zero row in A a . Therefore, for the case of Fig. 3, 2/3 of the computations are redundant and can be avoided.
To realize this, one may think to limit the domain for the entire assimilation step; then, the number of zero rows certainly would be largely reduced. This is actually incorrect, because these zero rows are changing along with the transport of ash clouds, and are not constant at each analysis step. So the full domain must be considered and it should be adaptive (choose different zero rows according to different A f at different analysis times).

Derivation of the mask-state algorithm (MS)
Here we introduce item n noash to represent the number of zero rows in the ensemble state matrix A f , and use n ash to represent the number of other rows (also, n ash represents the grid size of ash plume). When computing A a = A f X, to avoid all the computations related to n noash rows with zero elements, the index of other n ash rows must first be decided. This index is meant to reduce the dimensions of A f . After getting a A a with a dimension of n ash × N, the index will be used again to reconstruct the full matrix A a with the dimension of n × N. Based on this idea, we propose a mask-state algorithm (MS) which deals with the time-consuming analysis update. MS includes five steps.
i. Compute ensemble mean stateĀ f . The mean statē A f n × 1 can be easily computed by averaging A f n × N along N columns. Due to all elements in A f n × N corresponding to ash concentrations, all elements in A f n × N are larger than or equal to zero. The index of non-zero rows inĀ f n × 1 is thus equivalent to that in A f n × N . The computational cost for this step is O(n N).
ii. Construct mask array z. Based on previously obtained A f n × 1 , we search the non-zero elements ofĀ f n × 1 and record the index into a mask array z n ash × 1 . With this strategy, we do not need to search the full matrix A f n × N and build an index matrix for storage. This is a benefit for saving memory. The computational cost for this step is O(n).
iii. Construct masked ensemble state matrix A f . Using the mask array z n ash × 1 obtained from step (ii), A f n ash × N can be constructed column by column according to Eq. (13), and the computational cost (overhead) for this step is O(n ash N).
iv. Compute A a by multiplying A f and X. Perform matrix computation A a n ash × N = A f n ash × N X N × N . This step is similar to A a = A f X, as described in Sect. 3.2, but the computational cost now becomes O(n ash N 2 ) instead of O(n N 2 ).
v. Construct analyzed ensemble state matrix A a . With the computed A a from step (iv) and the mask array z from step (ii), the final analyzed ensemble state matrix A a n × N can be constructed based on Eq. (14). The computational cost (overhead) for this step is O(n N).
A a (z(1 : n ash ), 1 : N ) = A a (1 : n ash , 1 : N ) According to the derivations of MS, the computational costs related to zero rows are avoided. Here the "zero rows" do not equal "zero elements". The former corresponds to the regions where there is no ash for all the ensemble members, while the latter also counts the no-ash regions specifically for some ensembles. Certainly the consideration of all "zero elements" can include all the sparsity information of the ensemble state matrix, but extra computations and memories must be spent on searching the full matrix A f n × N with a computational cost of O(nN ) and storing a mask-state matrix with dimensions of n × N. This is expensive compared to constructing the mask array in procedure (ii). Actually, after a careful check of the volcanic ash ensemble plumes, there is no "bad" ensemble which is really different from others. Although the concentration levels in ensemble members are distinct, the main direction and the occurrence to the grid cells are more or less the same. This means that the "zero rows" actually more or less equal "zero elements" but are much faster than the way with "zero elements", which confirms the suitability and advantage of procedure (ii). Probably when there are big meteorological uncertainties, the "zero elements" will be much larger than "zero rows". In this case, how to make use of the sparsity information in the ensemble state matrix will be considered in future.
Based on procedures of MS, the computational cost of A a = A f X can be reduced. However, without a careful evaluation, we cannot conclude MS is fast, because the algorithm also employs other procedures. If these procedures (i), (ii), (iii), and (v) are much cheaper than the main procedure (iv), MS can definitely speed up the analysis step, and vice versa. Now we analyze MS's computational cost, which can be summed as O(n N) + O(n) + O(n ash N) + O(n ash N 2 ) + O(n N), i.e., O(n N +n ash N 2 ). Thus, the computational overhead involved in transforming the full matrix to a small one (i.e., O(n ash N ) for procedure (iii)) has little effect on the total computation cost of MS (i.e., O(n N +n ash N 2 )). However, the computational overhead of transforming the small matrix to the full one (i.e., O(n N) for procedure (v)) does contribute a part, which cannot be ignored, to the total MS's computational cost. The computational cost without MS is O(n N 2 ).
The comparison between both costs (with and without MS, i.e., O(n N + n ash N 2 ) and O(n N )) indicates when the number of non-zero rows (n ash , i.e., the number of grids with ash) of the forecasted ensemble state matrix satisfies n ash < N−1 N n; then, MS can accelerate A a = A f X. Here, O(n N + n ash N 2 ) and O(n N) are on the same order when n ash < N−1 N n. The larger the difference between n ash and N−1 N n, the better the speedup can be achieved. According to this analysis, and the characteristic (e.g., n ash n approximately equals 1 3 in this case) of volcanic ash transport as described in Sect. 4.1, the relation is certainly satisfied and is actually n ash N −1 N n (significantly smaller) for our study. Therefore, for our volcanic ash DA system, with MS, the computational cost for the time-consuming part A a = A f X is O(n ash N 2 ), which is much reduced compared to O(n N 2 ) with conventional computations.
The relation n ash < N−1 N n indicates whether we would have speedup by the MS method; actually, it can be extended to Eq. (15), which explicitly specifies the expected amount of speedup (S ms ) of A a = A f X by the MS algorithm. In this case study, N is taken at 100 and n ash n ≈ 1 3 , so S ms is approximately 3.0. According to Amdahl's law (Amdahl, 1967), the total computational speedup (S total ) by MS can be predicted by Eq. (16), where p ms is the proportion of the computational cost of A a = A f X in the overall DA computations. It has been evaluated that the computational cost of A a = A f X dominates the analysis step (see Fig. 2b); thus, the proportion of the computational cost of A a = A f X approximates the proportion of the analysis step in the total DA computations (i.e., p ms ≈ 72 % in this case, as described in Sect. 3.1). Therefore, based on Eq. (16), the maximum ("ideal") computational speedup can be predicted to be 1 1−p ms (i.e., ≈ 3.57 for this case study) when S ms approximates infinity. However, this is not the actual speedup because S ms is in fact specified by Eq. (15). Based on the discussions above, S total can therefore be estimated by Eq. (15) at ≈ 2.0 in this case.

Experimental results
Analysis of the algorithmic complexity of MS shows that MS is an efficient approach to reduce the computational cost of Geosci. Model Dev., 10, 1751-1766, 2017 www.geosci-model-dev.net/10/1751/2017/ the time-consuming A a = A f X. Now MS will be applied in the real volcanic ash DA system, to investigate whether in practice it can speed up the analysis step well. We perform MS in the conventional EnKF, which means initialization and forecast steps are all computed as the conventional EnKF. The only difference between MS-EnKF and conventional EnKF is that in the former MS is employed for the analysis step, and in the latter is the standard analysis step. The result and related specifications are shown in Table 1. As introduced in Sect. 2, the forecast step has been configured with the conventional parallelization; thus, N+2 (102 here) cores are actually used (one core for the DA algorithm, the other N+1 cores for the parallel forecast of N ensemble members and one ensemble mean). It can be seen from Table 1 that MS indeed largely accelerates the analysis step (as expected, by a factor of about 3.0 for this study), which confirms the theoretical cost evaluation. The detailed experimental time for each step of MS is shown in Table 2. As expected, the densedense matrix multiplication in step (iv) takes the largest part (i.e., 0.8474 h for this case study) of the total computational time (0.87 h) of MS. However, step (iv) has been a big improvement compared to the case without MS (3.14 h; see Table 1), which is because the computational time for the other steps (e.g., steps (i-iii) cost only 0.0156 h to reduce the size of the ensemble state matrix) is little and ignorable. Note that the total computational time of A a = A f X with MS (i.e., 0.87 h in Table 2) is not exactly equal to the computational time of the MS-EnKF analysis procedures (i.e., 0.88 h in Table 1). The subtraction (i.e., 0.01 h) corresponds to the summed computational time of all the other analysis procedures (i.e., procedures 1-8) except for A a = A f X (see Fig. 2b and Table 3).
MS is now experimentally proven as efficient to significantly reduce the computational time for the analysis step during volcanic ash DA. Note that it can also be observed that the computational time for the "other" parts in Table 1 (such as operations for setting environmental variables, starting and finalizing DA algorithms, as mentioned in Sect. 3.1) is slightly reduced by the MS method (i.e., 0.03 h in this case). This is because in the conventional EnKF, the ensem-ble mean stateĀ f is calculated in the "other" parts as an output to finalize the DA algorithms, while in MS-EnKF, the calculations ofĀ f are needed and directly involved in the "Analysis" part.
The result shows that, benefitting from the success of a reduced analysis step, the overall computational cost indeed gets significantly reduced. The total execution time is 1.95 h, which is less than the simulation window of 3 h (09:00-12:00 UTC, 18 May 2010). This result satisfies our goal to accelerate the computation to an acceptable runtime (i.e., requires less runtime than the time period of the DA application). Therefore, aviation advice based on the MS-EnKF can be provided as not only accurate, but also sufficiently fast. Note that the result (1.95 h) is obtained after the volcanic ash is transported to continental Europe. If the assimilation is performed in the starting phase of volcanic ash eruption (when aircraft measurements are available), a more significant acceleration would be obtained. This is because in this case the volcanic ash is only transported in an area near to the volcano; thus, the number of no-ash grid cells will take a large proportion (much higher than 2/3 for this case study) of the full domain.
There is another interesting point. According to Fig. 3, the ash grids comprise 39.3 % of the total grids. Thus, the minimum computing time by using MS to utilize this model's characteristic should be ≈ 1.234 h (i.e., 0.393 × 3.14 h). However, the experimental result shows that the computational time goes down to 0.88 h (see Table 1). One reason for this time decrease is that when the size of the matrix is reduced, the memory access cost also goes down (e.g., through better cache usages). Another possible reason is that the ash grid number actually decreases with time (not always taking 39.3 % of the total grid number), due to ash sedimentation and deposition processes (Fu et al., 2015).
Note that in this study we only perform the commonly used ensemble parallelization for the forecast step (already efficient compared to the expensive analysis step) but do not choose model-based parallelization (e.g., tracer or domain decomposition). As specified in Table 1, no parallelization is implemented on the six tracers. This is because due to the important aggregation process (Folch et al., 2010), there are big dependencies between different ash components and thus it does not make much sense to parallelize them. As for domain-decomposed parallelization (Segers, 2002), it is not efficient for our application. This is because volcanic ash is special in the sense that the model is only doing computations in a small part of the domain (i.e., there are no data in a rather large part of domain), and this active part is continuously changing. Thus, a fixed domain decomposition is not very useful here because of the changing plume position. In this sense, some advanced approach such as adaptive domain-decomposed parallelization (Lin et al., 1998) should be adopted to achieve additional acceleration to the volcanic ash forecast stage. This is an interesting subject for future application, when a more complicated model is employed, only ensemble parallelization may be not enough for the forecast stage.
5 Comparison between MS and standard sparse matrix methods

Issues related to the generation of CSR-based arrays
According to Sect. 4, MS has proven to be capable of solving the computational issue of A a = A f X. Motivated by the model's characteristics, MS was proposed from an application's perspective and achieved a good result by managing the irregular sparsity in our complicated volcanic ash DA system. The main reason why MS is efficient is that the sparsity of A f can be well utilized by MS. In Sect. 4, we only performed the comparison between MS and the case of full storage dense matrices. However, the problem abstracted here (A f X) is actually a sparse-dense matrix multiplication (SDMM) problem, since A f is sparse and X is dense (see Eq. 12 for X). Thus, one may wonder what the result would be if the comparison of MS is made to more standard sparse matrix methods, such as compressed sparse row (CSR)-based methods (Saad, 2003;Bank and Douglas, 1993), which are commonly used for sparse matrix vector/matrix multiplication.
Before we make the comparison, we need to first address the intrinsic problem when considering standard sparse matrix methods in EnKF for A a = A f X. The issue is that it is not possible to directly generate a sparse storage format (e.g., CSR) of A f without first generating the full matrix A f . This is mainly because A f comes from the model-driven ensemble forecast step, where each A f column corresponds to one member of the ensemble. During model forecast, we know there are indeed no-ash grids. However, it is not certain where the plume is exactly after one forecasting time step. This is highly dependent on the weather conditions and the model processes (e.g., advection and diffusion for horizontal grids, sedimentation and deposition for vertical grids). Thus, a fixed and wide domain is usually needed by the model to avoid complications, resulting in the generation of the full storage of A f (to be used in A a = A f X). Therefore, if we want to implement a CSR storage format for the sparse matrix A f , we must first generate the full storage A f from the ensemble forecast step, and then we generate the three CSR arrays based on A f . Generating CSR arrays is usually much more expensive (computationally) than a single sparse matrix-vector multiplication (SpMV). Thus, if we generate CSR arrays for only performing one-time SpMV, it would be meaningless from HPC's point of view. Fortunately, this is not the case for A a = A f X (i.e., SDMM), which can actually be considered as N-times SpMV. (Here, X has N columns, and one SpMV means the multiplication of A f by one column of X.) Thus, CSR-based SDMM might also be a candidate in reducing the computation time of A a = A f X. It remains interesting to compare the performance of CSR-based SDMM and MS in dealing with A a = A f X for our study case.

Result of CSR-based SDMM
To implement CSR-based SDMM for A a = A f X, the three CSR arrays for A f (denoted val, col_idx, and row_ptr in this study) need first to be generated. The array val of size n val stores non-zero values of A f , where n val ≈ n ash N . (In this study case, n = 3.888 × 10 6 , n ash ≈ 1 3 n, and N = 100.) The array col_idx of the same size n val stores the column index of the non-zeros. The array row_ptr saves the start and end pointers of the non-zeros of the rows in A f . The size of row_ptr is n + 1.
After the above three CSR arrays are generated, CSRbased SpMV can be performed for multiplying A f by a column vector v in X (see Algorithm 1 in Fig. 4a). With that, Algorithm 2 (Fig. 4b) can be implemented by looping Algo-Geosci. Model Dev., 10, 1751-1766, 2017 www.geosci-model-dev.net/10/1751/2017/ rithm 1 for N times to obtain A a = A f X. The experimental result of CSR-based SDMM is shown in Table 4, where all the environmental conditions (such as the DA system, the programming environment) are the same as the case of MS. This gives a fair comparison between CSR-based SDMM and MS. In addition, for a pure algorithmic comparison with the serial MS, here the CSR-based SDMM is also performed in a serial case.
From Table 4, we can first confirm that the computational time (i.e., 0.0407 h) for the generation of the three CSRbased arrays (val, col_idx, and row_ptr to represent the sparse matrix A f ) indeed takes more time than the computational time of one CSR-based SpMV (i.e., 0.0117 h). Thus, there is little value in performing sub-step (i) (see Table 4) if only one SpMV (i.e., sub-step (ii)) is needed. However, to get all N (i.e., 100) columns of A a , the sub-step (ii) is looped for N times, resulting in an ignorable impact of sub-step (i) on the total computational time (i.e., 1.21 h) of CSR-based SDMM.
The result of CSR-based SDMM also shows that the standard sparse matrix methods can reduce the computational time of A a = A f X, by comparing with the conventional way in Table 3. However, it can also be observed that the computational time of CSR-based SDMM is larger than MS (i.e., 1.21 h versus 0.87 h in Table 3). Thus, although application of sparse matrix multiplication methods is positive, it is still slower than MS on our problem.

Comparison between CSR-based SDMM and MS
In the CSR-based SDMM, only non-zero elements in A f participate in the multiplication between A f and X; thus, redundant computation (related to zero elements in A f ) is avoided. So the computation time of A f X is reduced with CSR-based SDMM. In the following, we analyze the performance difference between CSR-based SDMM and MS.
Firstly, from the programming's perspective, in CSRbased SDMM, the loop number for the rows of A f is from 1 to n (see Fig. 4a), while the corresponding loop number in MS is from 1 to n ash (see step (iv) of MS in Sect. 4.2, n ash ≈ (1/3)n). Although only non-zero elements are used in the multiplication in CSR-based SDMM, the length of the outer loop is still n (much larger than n ash ), which is the essential reason that MS is faster than CSR-based SDMM. Note that as discussed in Sect. 4.1, there are many zero rows in A f ; thus, CSR-based SDMM actually does nothing when it comes to a zero row, but still needs to execute the loop. Within each loop number, it has to check the information from row_ptr (size n + 1), where the value corresponding to a zero row is usually set to be the value in row_ptr corresponding to the first subsequent non-zero row.
Secondly, with respect to the algorithm, CSR-based SDMM utilizes the sparsity of A f by its generation of three CSR arrays, while MS not only utilizes the sparsity information of the sparse matrix A f , but also utilizes the consistency of ensemble forecasts; that is, ensemble forecasted states are not consistent in values but usually consistent in non-zero locations. This is a typical property in ensemble-based DA, resulting in N ensemble plumes being different in concentration values but having similar transport directions/shapes (see Sect. 4.2). Thus, most of the zero elements in A f are actually in zero rows of A f for an EnKF application, which leads to a small number of non-zero rows (n ash ) compared to the full number of rows (n) of A f . Therefore, only considering n ash rows in A f n ash × N X N × N (see step (iv) of MS in Sect. 4.2) is more advantageous for an EnKF application than considering all n rows in CSR-based SDMM. Based on the above analysis, MS can be considered a specific sparse matrix method, which typically works for ensemble-based DA applications.
It is useful to apply standard sparse matrix methods (e.g., CSR-based SDMM) for our assimilation application. The accelerated analysis step by CSR-based SDMM (1.22 h; see Table 3) also reduces the total computational time (i.e., 2.29 h; see Table 1 for the computational time of initialization, forecast, and others) to an acceptable level (i.e., less than 3 h for our case study). In practice, due to the better performance of MS than CSR-based SDMM, we will use MS as a better choice for assimilation applications. In addition, we do not only intend to present MS, but also intend to reveal which part is the most time-consuming part for plume-type assimilation of in situ observations. 6 Discussions on MS

Applicability
For volcanic ash forecasts, only a relatively small domain is polluted compared to the full 3-D domain, so that MS can work efficiently. Using MS is also applicable for many other DA problems, where the domain is not fully polluted by the species. It does not matter what the emission looks like and whether the releases are short-or long-lived species. Given an assimilation problem, the only restriction for MS to gain an acceleration is whether the whole domain is fully polluted or partly polluted. The assimilation problems where MS can achieve the acceleration effect on the computations of A a = A f X include all the volcanic-related ash/gas assimilations, e.g., assimilation of satellite data/LIDAR data/in situ data; (sand/desert) dust-storm-related assimilation; tornadorelated assimilation; assimilation of exploding nuclear plants or factories; chemicals or oils leaking into seas; global (forecast) fire assimilation; and assimilation of environmental pollutant transport, e.g., severe smog. In addition, for DA applications (e.g., ozone, SO 2 ) where pollutants spread over the whole domain, usually the focus is only on the high concentrations, and a threshold can be set to ignore the very low values without losing the necessary assimilation accuracy. In this case, MS can also lead to a potential acceleration since many very low concentrations can be explicitly truncated to be zeros. It has been analyzed that when the number of non-zero rows (n ash , i.e., the number of ash grids in a 3-D domain) of A f satisfies n ash < n, MS can work faster than standard EnKF. For volcanic ash application, because n ash is much less than n, the acceleration is quite large. Hence, in this case, we propose to embed MS in all ensemble-based DA methods because it is fast and the implementation using MS is exact to the standard ensemble-based methods; i.e., it does not introduce any approximation in view of MS procedures. Actually this proposal can be extended to all real applications, even if the condition is not satisfied. This is because in this case the computational cost of MS for A a = A f X becomes O(n N 2 ), which is the same as that of using the standard assimilation (shown in Fig. 2b). Therefore, if the state numbers are equal to or close to the total number of grid points in the domain, the added computational cost of using MS is very small (negligible), so that the computational time with MS is almost the same as the time of using the standard approach, whereas when the condition n ash < n is satisfied, MS will accelerate the analysis step. Thus MS is generic and can be directly used in any ensemble-based DA, and this acceleration can be automatically realized for some potential applications, without spending time investigating whether the condition is satisfied. In a real (or operational) 3-D DA system, MS can be easily included; i.e., we only need to invoke the MS module when computing A a = A f X, without any other change to the current framework. Note that MS is applicable in ensemblebased DA but not in variational-based DA. This is because in a variational-based DA system, the minimization of a cost function is mostly operated within several/many continuous time steps (Lu et al., 2016b; thus, it is convenient to always use the full (i.e., non-masked) domain to represent different state matrices (corresponding to different time steps in variational-based DA).
As stated in Eq. (15), the speedup of the MS method is approximately the inverse of n ash n . So far there are no statistical data on the value of n ash n . Considering the problem of volcanic ash transport, there is one emission point (at the volcano); all the ashes in atmospheres are transported by the directional wind drive from the same source point. Thus volcanic ash cloud is actually transported in a shape of a plume, which in general does not cover the full but only a small part of the 3-D domain. At the start phase of a volcanic ash eruption, n ash n is much smaller than 1.0 (started from 0). During transport over a long time (1.5 months for this case study), n ash n increases to approximately 1/3. Therefore, the speedup of MS in ensemble-based volcanic ash DA will be significant.

MS and localization
Based on the formulation of MS, one may think it can be taken as a localization approach (Hamill et al., 2001). There is indeed a similarity between MS and the localization approach, in a sense that when computing A a = A f X, both get rid of a large number of cells and only do computations related to the selected grids. These two algorithms are however functionally different. This is because the localization approach is meant for reducing spurious covariances outside a local region which is built up around the measurement; thus, the results with and without localization approaches are different, while MS is developed for the acceleration purpose. The masked region is discontinuous and independent of locations of measurement, but dependent on the model domain. Thus, there is no difference in the assimilation results between using MS and without using it. Therefore, based on the functional difference, MS cannot be taken as a localization approach.
In this study, we do not employ the localization strategy in the analysis step, because we use a rather large ensemble size of 100 to guarantee the accuracy, as introduced in Sect. 2. But for some applications (e.g., ozone, CO 2 , sulfur dioxide), especially when assimilating satellite data, localization is a necessary approach and has been widely used in reducing spurious covariances (Barbu et al., 2009;Chatterjee et al., 2012;Curier et al., 2012). In these cases, because the localization approach forces the analysis only to update the state within a localization region, one may think that localization could replace MS and that there would be no significance in employing MS. Actually this is not correct. We explain the reason as follows.

1998, 2001) given by
The Schur product f•P f in Eq. (17) is defined by the elementwise multiplication of the covariance matrix P f and a localization matrix f. f is defined based on the distance between two locations; thus, it is dependent on the domain and needs information on the full ensemble state locations. In this way, f • P f can contain more zeros than P f , but the dimensions are not changed, so that the computations related to f • P f are actually not reduced. Therefore, we can understand the localization approach in the analysis step as that the states within and outside a local region are both updated with increments, but just the increments outside the region are zero (which seems like not updating). This is also the reason why the localization approach is not meant for acceleration, but only for reducing spurious covariances. Now it is clear that localization cannot replace MS. Actually both can be performed together in dealing with the time-consuming part A a = A f X. The localization approach can first transfer A f to a localized matrix with more zero rows. Then MS can be used to accelerate the multiplication of the localized matrix and X. In this way, MS is expected to accelerate A a = A f X with a high speedup rate, because the computational cost of more zero rows in the localized ensemble state matrix is avoided.

MS and parallelization
Motivated by the model's physics, the implementation MS currently is for the serial case. This implementation has reduced the computation time to an acceptable time (i.e., the simulation time is less than the period of forecast in realworld time). It is however interesting to discuss the potential of parallelization of the dense-dense matrix multiplication ( A a n ash × N = A f n ash × N X N × N ) in step (iv) of the algorithm (see Sect. 4.2 and Table 2). The related matrix multiplication can be easily parallelized on multiple processors. Optimization and evaluation on the parallelized MS will be considered in future. For the current case study, the computational time (3.13 h; see Table 3) for an "ideal" reduction by parallelization of MS is not much larger than the acceleration (already) gained by MS (2.26 h, subtraction between 3.13 h and 0.87 h; see Table 3). Therefore, from the application's perspective, further acceleration by parallelization is not required.
Alternatively, one may also consider to (1) directly parallelize the expensive matrix multiplication of A a n × N = A f n × N X N × N , without first performing MS, or (2) implement CSR-based SDMM (see Sect. 5) with parallelization. Both are possible alternative approaches to accelerate the expensive matrix multiplication. The first approach can be implemented by a user's own designed parallelization, or by utilizing scaLAPACK (https://www.netlib.org/scalapack/, where the main function is "pdgemm"). The second approach can be realized by using some general parallel sparse-dense matrix multiplication methods (e.g., sending each column of X and three CSR arrays of A f to each processor to calculate each column of A a ) or using a good parallel algebra library like PeTSC (https://www.mcs.anl.gov/petsc/) which allows users to specify own orderings and comes with machine optimized parallel matrix-matrix multiplication operations. However, given the fact that MS can also be parallelized using similar ways or the same libraries, it is fair to not consider parallelization for all cases (i.e., using MS, not using MS, using CSR-based SDMM). Actually, the parallelization in MS could be performed much more easily than other approaches in dealing with A a = A f X, because the dense-dense matrix multiplication (parallelization in step (iv) of MS) is easier to parallelize than the sparse-dense matrix multiplication (direct parallelization for A a = A f X or parallelized CSR-based SDMM).
In this paper, for the current usage, we keep the possibility of parallelization open, because a serial MS has been efficient already.

Conclusions
In this study, based on evaluations of the computational cost of volcanic ash DA, the analysis step turned out to be very expensive. Although some potential approaches can accelerate the initialization and forecast steps, there would be no notable improvement to the total computational time due to the dominant analysis step. Therefore, to get an acceptable computational cost, the key is to efficiently reduce the execution time of the analysis step.
After a detailed evaluation of various parts of the analysis stage, the most time-consuming part was revealed. The mask-state algorithm (MS) was developed based on a study of the characteristics of the ensemble ash states. The algorithm transforms the full ensemble state matrix into a relatively small matrix using a constructed mask array. Subsequently, the computation of the analysis step was sufficiently reduced. MS is developed as a generic approach; thus, it can be embedded in all ensemble-based DA implementations. The extra computational cost of the algorithm is small and usually negligible.
The conventional ensemble-based DA with MS is shown to successfully reduce the total computational time to an acceptable level, i.e., less than the time period of the assimilation application. Consequently, timely and accurate volcanic ash forecasts can be provided for aviation advice. This approach is flexible. It boosts the performance without considering any model-based parallelization such as domain or component decomposition. Thus, when a parallel model is available, MS can easily be combined with the model to gain a further speedup. It implements exactly the standard DA without any approximation and with easy configurations, so that it can be used to accelerate the standard DA in a wide range of applications.
In this case study with the LOTOS-EUROS model (version 1.10), after the parallelization is performed for the forecast step of EnKF assimilation, the analysis step takes 72 % of the total runtime, which means the analysis step is the bottleneck. This case might not be general for all ash forecasts, as the computational cost for initialization and forecast greatly depends on the forecast model that is used. For the current development, it makes sense to use the LOTOS-EUROS model, because the model has been configured and evaluated in (Fu et al., 2015) by comparison with other famous models (e.g., NAME, Jones et al., 2007, andWRF-Chem, Webley et al., 2012) in simulating volcanic ash transport. However, if a more expensive ash forecasting model is used, then the bottleneck would be the forecast step. In this case, the forecast step should be the goal for acceleration, and probably a parallel model or adaptive domain decomposition (as discussed in Sect. 4.3) needs to be employed together with the parallel ensemble forecasts.
The use of in situ measurements is one important reason why MS works perfectly. For each analysis step, the number of measurements are quite small, and the procedure of the singular value decomposition (SVD) costs little. However, in some applications when many measurements are assimilated (e.g., satellite-based data Fu et al., 2017 or seismicbased data Khairullah et al., 2013), and the number of measurements is on the same order as the number of state variables, the most time-consuming part will be the SVD. In these cases, the contributions of MS will be limited. The reduction of the total computing time using MS is therefore less significant; an effective acceleration algorithm for the analysis step must be used and should consider the computationally expensive SVD in the first place.