Assimilation of MODIS Dark Target and Deep Blue observations in the dust aerosol component of NMMB-MONARCH version 1.0

Abstract

A data assimilation capability has been built for the NMMB-MONARCH chemical weather prediction system, with a focus on mineral dust, a prominent type of aerosol. An ensemble-based Kalman filter technique (namely the local ensemble transform Kalman filter – LETKF) has been utilized to optimally combine model background and satellite retrievals. Our implementation of the ensemble is based on known uncertainties in the physical parametrizations of the dust emission scheme. Experiments showed that MODIS AOD retrievals using the Dark Target algorithm can help NMMB-MONARCH to better characterize atmospheric dust. This is particularly true for the analysis of the dust outflow in the Sahel region and over the African Atlantic coast. The assimilation of MODIS AOD retrievals based on the Deep Blue algorithm has a further positive impact in the analysis downwind from the strongest dust sources of the Sahara and in the Arabian Peninsula. An analysis-initialized forecast performs better (lower forecast error and higher correlation with observations) than a standard forecast, with the exception of underestimating dust in the long-range Atlantic transport and degradation of the temporal evolution of dust in some regions after day 1. Particularly relevant is the improved forecast over the Sahara throughout the forecast range thanks to the assimilation of Deep Blue retrievals over areas not easily covered by other observational datasets.

The present study on mineral dust is a first step towards data assimilation with a complete aerosol prediction system that includes multiple aerosol species.

Introduction

Among the different aerosol species, mineral dust is one of the main components of the atmospheric aerosol loading and is of great interest for a variety of reasons. Mineral dust plays an important role in the earth's energy balance and has a relevant impact on economical activities, on the ecosystem, on health, as well as on weather and climate . The strong dust storms occurring near emission sources constitute a big hazard to air traffic and road transport as they can reduce the visibility down to few metres. However, dust does not affect only local economies: because of its transport over thousands of kilometres, it has an impact from local to global scales. Dust deposition provides nutrients to continental and marine ecosystems. African dust for example has a role in fertilization of the Amazon rainforest , while dust deposition over oceans has implications for sea biogeochemistry as the iron contained in the dust particles is a nutrient for phytoplankton, whose photosynthetic activity in turn affects the carbon cycle . Dust has health implications both close to and far from sources. For example, studies have shown the usefulness of dust aerosol climatologies in predicting part of the year-to-year variability of the seasonal incidence of meningitis in Niger , while particulate matter measurements taken in areas far from sources show that Saharan dust outbreaks have adverse effects of cardiovascular and respiratory conditions . Mineral particles perturb the earth-atmosphere's radiation budget through their interaction with the short-wave radiation, through scattering and absorption, and long-wave radiation, through absorption and re-emission. Due to this redistribution of the energy, dust aerosols can have a strong impact on atmospheric processes at short (weather) and long (climate) term periods, while they can affect atmospheric circulations at large spatial scales (e.g. Asian monsoons; ). Furthermore, this can generate feedback processes since changes in weather and climate can in turn lead to changes in the dust cycle.

Different types of ground-based (e.g., ) and space-borne (e.g., ) observations have been utilized to describe the variability of atmospheric dust. However, due to either insufficient spatial representativeness or accuracy, the spatio-temporal features of dust aerosols are not fully captured by the current observing system. Neither do models accurately describe atmospheric and surface dust concentrations . High uncertainties are also in our knowledge of the optical and micro-physical properties of dust, and in our representation of its vertical structure. The latter has implications for the radiation's budget and transport. On the other hand, an accurate quantification of dust's spatial and temporal distribution is key in correctly characterizing the effect that it has on the earth's energy balance, as well as in improving the skill in forecasting its concentrations in the atmosphere as well as in forecasting the weather .

Regional and global centres, predicting the most important aerosol species or dust only, participate in different model inter-comparison initiatives like the Aerosol Comparisons between Observations and Models (AeroCom; ) project, the International Cooperative for Aerosol Prediction (ICAP; ) initiative, and the WMO Sand and Dust Storm Warning Advisory and Assessment System (SDS-WAS; ). Multi-model ensemble spreads give an indication of large uncertainties in the modelling schemes and confirm the need for a better characterization of aerosols. Relatively recently because of these large uncertainties, the atmospheric composition community has begun to make use of data assimilation (DA) for better characterization and prediction of atmospheric constituents such as aerosols and trace gases . Though their dynamic is mainly driven by forcings such as emissions, recent studies showed that the usage of observations through data assimilation has improved significantly the accuracy of short-term forecast and the global monitoring of both aerosols and trace gases . Since the first experiments in the early 2000s, the assimilation of aerosol observations is now operational in some of the main aerosol forecasting centres . have highlighted in particular the importance of a combined assimilation of satellite products for aerosol forecast.

The Earth Sciences Department of the Barcelona Supercomputing Center (ES-BSC) is implementing a gas-aerosol module able to predict atmospheric composition at different spatial and temporal scales within the NMMB (Non-hydrostatic Multi-scale Model on the B grid; ) state-of-art meteorological model. This modelling system is known as the Multiscale Online Nonhydrostatic AtmospheRe CHemistry mode (NMMB-MONARCH). We report here on the extension of NMMB-MONARCH with a data assimilation functionality using satellite aerosol optical depth. NMMB-MONARCH version 1.0, as in where the model was previously named NMMB/BSC-CTM, considers dust only, but other aerosols are being implemented . The focus of this work on mineral dust is justified by the operational services provided by NMMB-MONARCH. This model provides an operational dust forecast for the Barcelona Dust Forecast Centre under an initiative of the World Meteorological Organization. It participates in the multi-model dust ensemble of the aforementioned ICAP initiative, providing daily global dust forecasts of up to 120 h. It also provides daily regional forecasts of up to 60 h to the WMO SDS-WAS system. Before this work, the system did not have an aerosol data assimilation capability and dust was produced uniquely from model estimated surface emission fluxes. The present study on mineral dust is a first step towards data assimilation with a complete aerosol prediction system that includes multiple aerosol species (not only dust but also sea salt, sulfate, and organics).

Previous studies of assimilation of dust aerosol only have been conducted for the Chinese Unified Atmospheric Chemistry Environment – Dust (CUACE/Dust) forecast system . These studies have used variational data assimilation techniques (3D-Var) which require, in their most practical implementation, pre-calculated and constant in time model error structures. Alternatively, ensemble-based techniques use flow-dependent model error amplitudes and structures which evolve during forecast and, in theory, should be able to capture better instabilities in the background flow . Dust AOD is currently assimilated at the UK Met Office with a hybrid variational data assimilation technique (hybrid 4D-Var).

In this work we present the coupling of NMMB-MONARCH with an ensemble-based technique known as local ensemble transform Kalman filter (LETKF; ). The LETKF scheme has been shown to be particularly suitable for the assimilation of aerosol information since it has been observed by , , and that aerosol fields have limited spatial correlations. Long-range transport of dust could be an exception to this. Since detailed studies of spatial correlation length scales for dust long-range transport are still missing in the literature, in this work we assume that what has been derived (limited spatial correlations) in general for aerosols is valid for dust. The utility of ensemble-based techniques for global aerosol simulations has been shown in previous studies (; and more recently ). The main novelty in our study is the creation of the ensemble, our implementation is based on known uncertainties in the physical parametrizations of the sophisticated dust emission scheme used by the NMMB-MONARCH model, as well as in the use of observations particularly relevant for dust applications, like MODIS Deep Blue.

The NMMB-MONARCH chemical weather prediction system is described in more detail in Sect. , with emphasis on its dust module. A description of the data assimilation scheme and of the assimilated observations follows respectively in Sects.  and . We report then in Sect.  about the characteristics of the simulations that we have run, in Sect.  about the evaluation methodology that we have followed, and in Sect.  about the evaluation results of our simulation experiments. The final section concludes the paper with a summary of this development, the main results achieved, and future perspectives.

The NMMB-MONARCH model and its mineral dust component

The ES-BSC is implementing a new gas-aerosol module within the NMMB meteorological model from the Unites States National Centers for Environmental Prediction (NCEP). The new modelling system is known as NMMB-MONARCH where it was previously named NMMB/BSC-CTM, and is developed in collaboration with NCEP and other research institutions. The chemistry (aerosols included) and meteorology are fully online integrated. NMMB-MONARCH is able to work with a wide range of spatial scales thanks to its unified non-hydrostatic dynamical core, keeping consistent parametrizations at different spatial and temporal scales. Furthermore, the dynamical core and the coupled modules are computationally highly efficient, satisfying current and projected operational requirements. The rest of this section will briefly describe some characteristics of the dust component of NMMB-MONARCH, with particular focus on the emission scheme.

The dust emission scheme implemented in NMMB-MONARCH follows the empirical relationship of and according to which the vertical dust flux is proportional to the horizontal sand flux. The horizontal to vertical flux ratio reflects the availability of dust in four soil populations (clay, silt, fine/medium sand, and coarse sand) . The horizontal sand flux is modelled as the flux of the saltating particles $H$ simulated according to and is proportional to the third power of the wind friction velocity. A soil moisture-dependent threshold on the friction velocity determines the velocity above which the soil particles begin to move in horizontal saltation flux. This threshold is dynamically estimated according to soil characteristics. Soil moisture effects are included following through the soil moisture correction parameter included in the expression for the threshold on the friction velocity. A sectional approach is used for the transport of dust particles, i.e. the dust size distribution is decomposed into small size bins. More exactly, dust is modelled using eight dust size bins varying from 0.1 to 10 microns, and, within each transport bin, dust is assumed to have a time-invariant lognormal distribution . The total vertical flux mass is distributed among the dust transport bins according to a specific dust distribution at sources. NMMB-MONARCH uses a distribution over sources derived from which assumes that the vertical dust flux is size distributed according to three lognormal background source modes. More explicitly, the dust vertical mass flux $F_b$ (kg s${}^{-\mathrm{1}}$ m${}^{-\mathrm{2}}$) in a given transport bin $b$ at each grid cell is given by $$F_b=CS(\mathrm{1}-V)\mathit{\alpha }H\sum _{i=\mathrm{0}}^{\mathrm{3}}{m}_i{M}_{i,b}b=\mathrm{1},\mathrm{\dots },\mathrm{8},$$ where $S$ is a source erodibility factor defined on bare ground surfaces, representing the probability of having accumulated sediments in the given grid cell that are potential dust sources; $(\mathrm{1}-V)$ is the grid's fraction of bare soil; $\mathit{\alpha }$ (m${}^{-\mathrm{1}}$) is the horizontal to vertical flux ratio calculated for four soil population classes (clay, silt, fine/medium sand, and coarse sand); $H$ (kg s${}^{-\mathrm{1}}$ m${}^{-\mathrm{1}})$ is the horizontal sand flux; $M_{i,k}$ is the mass fraction of background source mode $i$ carried in a transport bin $k$ calculated following , and weighted by a specific background source mode coefficient $m_i$; and $C$ is a global tuning factor empirically set to 0.768, which is meant to compensate for the uncertainty associated with the different components of $F_k$. More details about the above formulation of dust emission can be found in .

The mineral dust module has been extensively evaluated in studies at global and regional scales , showing that its evaluation scores lie in the upper range of the AEROCOM model evaluation performance scores. However, these evaluation efforts confirmed, similarly to other modelling systems, different sources of uncertainty in the NMMB-MONARCH dust modelling.

The data assimilation scheme

We have coupled NMMB-MONARCH with the LETKF scheme with four-dimensional extension as described in , in order to estimate optimal initial conditions for our dust model. The overall scheme implements an iterative approach consisting in a forecast step and state estimation step. The state estimation step combines information from mineral dust observations and a prior first guess, or background from our model. A short-term forecast is used as background information. The background incorporates information from past observations; therefore, the analysis is estimated using both current and past observations. LETKF is a development of the ensemble-based transform Kalman filter (ETKF; ) and of the local ensemble Kalman filter (LEKF; ), and is particularly suited to high-performance computing applications. A very attractive feature of an ensemble-based technique is the use of a flow-dependent background error covariance, which is derived from the ensemble of model states at the assimilation time, and evolves during forecast. At any given time in fact the state estimate is represented by an ensemble of system states and its uncertainty is derived from the ensemble. LETKF has the advantageous feature that it applies localization, i.e. it performs the analysis locally (at each grid point only observations within a certain distance are assimilated), allowing the global analysis to explore a much higher-dimensional space than the subspace spanned by the ensemble (whose dimensionality is limited by the number of ensemble members). Localization also reduces the effect of spurious long-range covariances, effectively eliminating them after a given distance. This is particularly suitable for the assimilation of aerosol information since, as mentioned in the introduction, it has been observed that aerosol fields have limited spatial correlations ($\sim$ 100 km). have already shown the positive impact of assimilating aerosol ground station observations using a LETKF assimilation system for the SPRINTARS model, while used it to assimilate CALIOP vertical profiles in the MASINGAR model and used it to ingest MODIS observations in the NICAM-SPRINTARS model.

Here we discuss the basic concepts behind the LETKF algorithm; a more detailed description of the scheme can be found in . Consider a state vector $\mathit{x}$ of the dynamic variables of a system (for our application these are dust mass mixing ratios). The mean analysis increment at a grid point is estimated as a linear combination of the background ensemble perturbations ${\mathbf{X}}^{\text{b}}$: $${\stackrel{\mathrm{‾}}{\mathit{x}}}^{\text{a}}={\stackrel{\mathrm{‾}}{\mathit{x}}}^{\text{b}}+{\mathbf{X}}^{\text{b}}\mathbf{w},$$ where we use the superscripts “a” and “b” to denote respectively the analysis and background state vector, and where the $i$th column of the matrix ${\mathbf{X}}^{\text{b}}$ is ${\mathit{x}}^{\text{b}\left(i\right)}-{\stackrel{\mathrm{‾}}{\mathit{x}}}^{\text{b}}$, $\mathit{\{}i=\mathrm{1},\mathrm{2},\mathrm{\dots },k\mathit{\}}$ with $k$ ensemble members, i.e. the difference between the $i$th ensemble forecast ${\mathit{x}}^{\text{b}\left(i\right)}$ and the ensemble forecast mean ${\stackrel{\mathrm{‾}}{\mathit{x}}}^{\text{b}}$. $\mathbf{w}$ is termed the “weight” matrix specifying what linear combination of the background ensemble perturbations is added to the background mean to obtain the analysis ensemble. The “weight” matrix is given by $$\mathbf{w}=[{\mathbf{Y}}^{\text{b}}{\mathbf{R}}^{-\mathrm{1}}{\mathbf{Y}}^{\text{b}}+(k-\mathrm{1}\left)\mathbf{I}{]}^{-\mathrm{1}}{\mathbf{Y}}^{\text{b}}{\mathbf{R}}^{-\mathrm{1}}\right({\mathit{y}}^{\text{o}}-{\stackrel{\mathrm{‾}}{\mathit{y}}}^{\text{b}}),$$ where ${\mathbf{Y}}^{\text{b}}$ is the background ensemble perturbation matrix in observation space (or background observation ensemble perturbation matrix), $\mathbf{R}$ is the observation error covariance matrix which we assume diagonal, $\mathbf{I}$ is the identity matrix, ${\mathit{y}}^{\text{o}}$ is the vector of observations, and ${\stackrel{\mathrm{‾}}{\mathit{y}}}^{\text{b}}$ is the mean background observation ensemble. The background observation ensemble is obtained by applying the observation operator $h(\cdot )$ to the ensemble forecast members ${\mathit{x}}^{b\left(i\right)}$, i.e. ${\mathit{y}}^{b\left(i\right)}=h\left({\mathit{x}}^{b\left(i\right)}\right)$.

LETKF uses R-localization, i.e. the localization is performed in the observation error covariance matrix, making the influence of an observation on the analysis decay gradually toward zero as the distance from the analysis location increases. To achieve this, the observation error is divided by a distance-dependent function that decays to zero with increasing distance: $e^{-\frac{{\text{dist}}^{\mathrm{2}}}{l^{\mathrm{2}}}}$, where “dist” is the distance in the grid space between an observation and the model grid in which the analysis is calculated, and $l$ is the horizontal localization factor.

Ensemble perturbations

We run the data assimilation scheme under an imperfect model scenario assumption: each ensemble member is run with a different perturbation of uncertain model parameters in the dust emission scheme. In dust modelling, the emission source term is a particularly large contributor to model error . In the case of NMMB-MONARCH one of the components of the uncertainty in the emission term has been identified for example in the vertical flux distribution at sources .

The model ensemble is created by perturbing the vertical flux of dust in each of the eight dust bins. As described in Sect. , NMMB-MONARCH follows a sectional approach for dust, i.e. the size distribution is decomposed into small size bins that from bin 1 to bin 8 go from 0.1 to 10 $\mathrm{\mu }$m with division intervals at 0.18, 0.3, 0.6, 1, 1.8, 3, and 6 $\mathrm{\mu }$m. This is equivalent to perturbing the total vertical flux as well as its size distribution at sources. The perturbations are extracted imposing some physical constraint: correlated noise is used across the bins so that noise correlation decreases with increased difference of the normalized cubic radius among the bins; the noise is applied multiplicatively and has mean 1 and standard deviation of 30 % of the unperturbed value in each bin; and the final distribution is unimodal. Figure  shows how the vertical flux is perturbed in our ensemble simulations. Additionally, we have perturbed the threshold friction velocity for dust emission by extracting a multiplicative random factor from a normal distribution with mean 1 and spread 0.4. This considers the uncertainty of the model with respect to both surface winds and soil humidity. At low resolution, model surface winds are typically underestimated over dust sources. Also, the model uses the scheme of to calculate the increase in the threshold friction velocity with soil humidity, which is typically overestimated in arid regions . The spin-up period for the ensemble ensures that perturbations applied at the sources propagate everywhere in the globe. For this reason at this first stage of development of our ensemble system we did not consider a combined meteorology and source perturbation necessary. The structure of our source perturbations, for both types of perturbations, is temporally and spatially constant.

Distribution of the mass vertical flux at sources across the eight dust transport bins for the different ensemble members in different colours, where the bin sizes from bin 1 to bin 8 go from 0.1 to 10 $\mathrm{\mu }$m with division intervals at 0.18, 0.3, 0.6, 1, 1.8, 3, and 6 $\mathrm{\mu }$m. The distribution derived from , and used in the standard forecast, is the dashed red line, with horizontal bars indicating the standard deviation of the noise used to create the perturbations. The mean of the ensemble perturbations is the dash-dotted line.

[Figure omitted. See PDF]

Observation operator

Our state vector is the dust mass mixing ratios. Therefore an observation operator is needed to map the ensemble mean state vector into the observation space. The simulated AOD at wavelength $\mathit{\lambda }$ is calculated at a given observation location according to the following linear operator: $${\text{AOD}}_{\mathit{\lambda }}=\sum _{b=\mathrm{1}}^{\mathrm{8}}{\displaystyle \frac{\mathrm{3}}{\mathrm{4}{\mathit{\rho }}_b{r}_b}}{M}_b{Q}_{\mathit{\lambda }b},$$ where ${\mathit{\rho }}_{\text{b}}$ (kg m${}^{-\mathrm{3}}$) is the particle mass density, $r_{\text{b}}$ (m) is the effective radius, $M_{\text{b}}$ (kg m${}^{-\mathrm{2}}$) is the dust column mass loading calculated from each dust bin mixing ratio, and $Q_{\mathit{\lambda }b}$ is the extinction efficiency factor which is calculated for using the Mie scattering theory assuming dust spherical, non-soluble particles, and, within a bin, a lognormal distribution for dust with a geometric radius of 0.2986 $\mathrm{\mu }$m and a standard deviation of 2.0.

When using in the state vector the total mass mixing ratio, as we will explain in Sect. , an ensemble averaged extinction efficiency is calculated as in as an average of the extinction efficiency of the individual bins weighted by the bin mixing ratios.

Hereafter, when we will use the term AOD without specifying the wavelength, we imply that we refer to aerosol optical depth at a wavelength of 550 nm, which is the most commonly reported value in the literature.

Observational data

MODIS Dark Target and OMI

We consider for assimilation the MODIS Level 3 AOD product produced by the US NRL and the University of North Dakota (hereafter called NRL MODIS). The NRL MODIS product is produced for the purpose of assimilation into aerosol transport models , post-processing the Level $\mathrm{2}$ MODIS Dark Target product from the so-called Collection 5 , and is available over both land and ocean. The MODIS Level $\mathrm{2}$ product is an average of the 1 km by 1 km retrievals (at nominal resolution) generated by the Dark Target algorithm applied to top-of-atmosphere reflectances observed by the MODIS sensor onboard NASA polar-orbiting satellites Terra and Aqua. The NRL MODIS Level $\mathrm{3}$ product is filtered and corrected in order to eliminate outliers and gross systematic biases, spatially aggregated into a 1${}^{\circ }$ by 1${}^{\circ }$ mesh in order to avoid the assimilation of sub-grid features, and an error is estimated for each observation. The product is generated every 6 h at 00:00, 06:00, 12:00, and 18:00 UTC and is based on MODIS Level $\mathrm{2}$ observations in a 6 h interval around those times. The retrieval errors estimated by NRL/University of North Dakota were used for the observation errors. They include the instrumental error variance and the spatial representation error variance. Following the approach in , we assume uncorrelated observation errors. These observations pertain to the total atmospheric aerosol column, not just to dust AOD. The selection of observations in dust-dominated conditions is performed using retrievals of the Ångström exponent (AE) from the original MODIS Level $\mathrm{3}$ product , for information about the size of the particles, and using retrievals of the Aerosol Absorbing Index (AAI) from the Ozone Monitoring Instrument (OMI) sensor , for information about the absorption characteristics of the particles. Ångström exponent (AE) values are based on quality-assurance-weighted 470 and 660 nm optical depths over land, and 550 and 865 nm optical depths over sea. Observations are selected when daily MODIS Aqua or Terra products contain a value for AE $<$ 0.75 and daily OMI products contain a value for AAI $>$ 1.5. Figure  shows an example of the NRL MODIS Level $\mathrm{3}$ product for a day of May 2007 after the filter for dust-dominated conditions is applied.

Aerosol optical depth (top) and its associated observation error (bottom) for 10 May 2007 for the NRL MODIS Level $\mathrm{3}$ product after the application of a filter for dust-dominated conditions.

[Figure omitted. See PDF]

MODIS Deep Blue

The MODIS Dark Target product does not provide information over very bright reflective surfaces, including deserts, as the retrieval algorithm assumes low surface albedo. We consider the assimilation of the MODIS Deep Blue Level $\mathrm{3}$ daily AOD product from Collection $\mathrm{6}$ whose algorithm retrieves AOD also over bright arid land surfaces, such as deserts. The Collection $\mathrm{6}$ product covers all cloud-free and snow-free surfaces, and can be potentially very useful for mineral dust applications as it is able to provide an observational constraint close to dust sources. The Deep Blue algorithm uses top-of-atmosphere reflectances at $\mathrm{412}$ and $\mathrm{470}$ nm. In the presence of heavy dust load the reflectance at 650 nm is also used. The algorithm exploits the fact that, over most surfaces, a darker surface and stronger aerosol signal is seen in the blue wavelength range than at longer wavelengths. The quality of the MODIS Deep Blue AOD product is improved in Collection $\mathrm{6}$ compared to Collection $\mathrm{5}$, as the work of , based on Level $\mathrm{2}$ retrievals, showed. Similar findings, for the northern African and Middle East deserts, were reported by , who used Level $\mathrm{3}$ retrievals over the period 2002–2014.

We have applied to this product the same filter for dust-dominated conditions described in Sect. . In addition we have masked out Level $\mathrm{3}$ retrievals obtained with less than $\mathrm{30}$ Level $\mathrm{2}$ retrievals, since showed that the agreement between MODIS-AERONET is improved when the sub-pixel spatial representativeness is increased. The MODIS Deep Blue observations are not corrected for possible systematic biases; however, we are aware that for future applications we should address any possible bias in the product. It is important to notice that the Level $\mathrm{3}$ product is an aggregation of Level $\mathrm{2}$ retrievals that is produced using the highest-quality retrievals (i.e. retrievals with quality-assurance flag value $\mathrm{3}$). Furthermore, we have applied a quality control to all the assimilated observations based on normalized first-guess departures. As a proxy for the normalization factor, we have used the standard deviation of first-guess departures.

Aerosol optical depth (top) and its associated observation error (bottom) for 10 May 2007 for the MODIS Deep Blue Collection $\mathrm{6}$ Level $\mathrm{3}$ product after the application of a filter for dust-dominated conditions.

[Figure omitted. See PDF]

Characteristics of the simulation runs.

Experiment	Ensemble	Dust initial conditions	Spin-up	Dust initial conditions
name	configuration	at 00:00 UTC on day 1	period	at 00:00 UTC after day 1
Control	No	Cold start	1 month	FC $+$ 24
				from previous day run
ENS-free-run	Yes	Warm start	11 days	FC $+$ 24
		from Control		of the individual members
				from previous day run
DA-NRL	Yes	Warm start	None	Analysis at 00:00 UTC
		from ENS-Free-run		of the individual members
				from previous day DA cycle
DA-NRL-DB	Yes	Warm start	None	Analysis at 00:00 UTC
		from ENS-Free-run		of the individual members
				from previous day DA cycle
AN-initialized	No	Warm start	None	Ensemble mean analysis
		from Control		from DA-NRL-DB

A study by shows that the highest-quality data have an absolute uncertainty of approximately (0.086 $+$ 0.56 AOD${}_{\mathrm{550}}$) $/$ AMF, where AMF is the geometric air mass factor with a typical AMF value of 2.8. We have used this quantification of the uncertainty for the Level $\mathrm{3}$ product. Furthermore, we have defined the representation component of the error as the standard deviation of the values used in the aggregated product. Although a more accurate treatment for the representation error could be envisaged following for example the approach of , we deem small the impact that our approximation has on the analysis. Figure  shows an example of the MODIS Deep Blue Collection $\mathrm{6}$ Level $\mathrm{3}$ product for a day of May 2007 after the filter for dust-dominated conditions is applied.

The number of MODIS Deep Blue and Dark Target observations used over the experimental period is shown in Fig. .

AERONET

For validation purposes we have used observations from the ground-based stations of the global Aerosol Robotic Network (AERONET; ) of direct-sun photometers. These observations have not been assimilated in our test simulations. In particular, we have used their retrievals of column-integrated aerosol optical depth from direct-sun photometric measurements. The retrievals are obtained observing the extinction of direct solar radiation due to the presence of aerosols in the atmosphere. For this reason AERONET retrievals are not available under cloudy sky conditions and during night-time. These observations suffer from a relatively sparse spatial coverage but are very valuable for validation purposes as their uncertainty in these retrievals is estimated to be between $\mathrm{0.01}$ and $\mathrm{0.02}$ . Several studies have in fact used the AERONET data for validation purposes, or for the correction of biases in satellite measurements . We considered cloud-screened and quality-assured (Level $\mathrm{2.0}$) direct-sun AOD retrievals between $\mathrm{440}$ and $\mathrm{870}$ nm. AERONET AOD at $\mathrm{550}$ nm was obtained using the Ångström law.

Number of NRL MODIS and MODIS Deep Blue Level 3 observations assimilated between May and August 2007.

[Figure omitted. See PDF]

Numerical simulation set-up

We have run a set of different experiments (listed in Table ): a control experiment to produce a 5-day forecast (hereafter called the Control experiment) with the same operational configuration (but at a coarser resolution) and version that provides daily global forecast to the aforementioned ICAP multi-model ensemble, and which is initialized for dust from the previous day 24 h forecast (FC $+$ 24). Assimilation experiments were run with NRL MODIS AOD (hereafter called the DA-NRL experiment) and with NRL MODIS AOD and MODIS Deep Blue AOD (hereafter called the DA-NRL-DB experiment) with a pre-processing to the observations as described in Sect. . Additionally, we have also run free ensemble simulations without assimilating any observation (hereafter called ENS-free-run). We have also run a 5-day forecast experiment initialized from the analysis produced by the DA-NRL-DB experiment (hereafter called the AN-initialized experiment) in order to evaluate the impact of the analysis on the forecast. The Control experiment was run for 5 months in the spring/summer period of 2007 (from 1 April to 31 August 2007) starting from a cold start for dust and with a spin-up period of 1 month (April 2007) which is excluded from the analysis of the results. Also, the ensemble is spun up before data assimilation is applied.

We use a 24 h assimilation window and observations are considered for assimilation at four time slots within the window, at 00:00, 06:00, 12:00, and 18:00 UTC. The system uses as first guess a 1-day forecast with output every $\mathrm{6}$ h. Simulated observation and background departures are calculated at each time slot. The time slots are exactly the ones corresponding to the times in which NRL MODIS AOD observations are available. We are using the LETKF implementation with a four-dimensional extension as described in . The state vector comprises the mixing ratio at all the time slots considered and so does the observation AOD vector. Background observation means ${\stackrel{\mathrm{‾}}{\mathbf{y}}}_j$ and perturbation matrices ${\mathbf{Y}}_j$ are formed at the various time slots $j$ when the observations are available. They are then vertically concatenated to form a combined background observation mean $\stackrel{\mathrm{‾}}{\mathbf{y}}$ and perturbation matrix $\mathbf{Y}$. $\stackrel{\mathrm{‾}}{\mathbf{y}}$ and $\mathbf{Y}$ are used for the calculation of a weight matrix $\mathbf{w}$ using the standard LETKF, i.e. we calculate a single $\mathbf{w}$ based on all innovations throughout the day. This same $\mathbf{w}$ is then applied to the state vector at different times throughout the assimilation window.

We have tuned different aspects of the data assimilation system, including testing the number of ensemble members, different perturbations of the ensemble, and a different state vector for the control variables. Using $\mathrm{24}$ ensemble members did not produce a significant impact on the dust analysis compared to the use of $\mathrm{12}$ ensemble members. This could be explained by our setting of a localization factor which makes the influence of an observation on the analysis decay gradually toward zero as the distance from the analysis location increases. We have set the horizontal localization factor to the value $\mathrm{1}$ in all the data assimilation experiments. This means that after two grid points the localization function is very close to zero. The value chosen is in the range of the ones used in previous studies such as and . Covariance localization in fact effectively eliminates background spatial correlations after a certain distance, and might have solved possible sampling errors introduced by the low dimensionality of the 12-member ensemble compared to the $\mathrm{24}$-member ensemble. We also apply vertical localization following the approach of localizing the error covariance vertically for radiance assimilation. The observation error is divided by the square of the model AOD normalized sensitivity function.

We have tested the usage of different perturbations of the dust emission scheme: a perturbation of the mass vertical flux per dust bin, or a perturbation of both the mass vertical flux and the threshold on the wind friction velocity. As we show in the next section, the latter configuration was deemed better as it spans a larger space of possible system states.

We have tested two different options for the state vector: a control variable consisting of the mixing ratio of eight individual dust bins or the total dust mixing ratio defined as the sum of the eight dust bins at each grid point and for all the vertical levels. In the latter case the mixing ratios for the individual dust bin after data assimilation are determined from the background, i.e. from their relative fractions before assimilation. The observation operator is calculated using the original mixing ratio following the approach for the observation operator in . The tests that we have performed show that representing individually the bins in the state vector does not have any significant impact on the analysis, while it increases the computational cost of the assimilation compared to using the total mixing ratio. Moreover, the use of a bulk approach is common in systems assimilating total AOD values as the observations are not able to fully constrain the individual bin profiles. We should note that this choice of state vector makes still meaningful the creation of the ensemble perturbing the vertical flux for the individual bins, as this allows us to express in the background the uncertainty in the size distribution at sources, and to span a larger space of possible system states.

In the next section we show the results of assimilating NRL MODIS NRL and MODIS Deep Blue observations using 12 ensemble members obtained by perturbing the mass vertical flux per bin at sources together with the threshold on the wind friction velocity, as described in Sect. , and using the total dust mixing ratio as the analysis variable in the state vector. All simulations were run on a global domain with 40 hybrid pressure-$\mathit{\sigma }$ layers, 5 hPa top pressure, and a horizontal resolution of $\mathrm{2.8}$${}^{\circ }$ by $\mathrm{2}$${}^{\circ }$. The NCEP final analysis at $\mathrm{1}$${}^{\circ }$ by $\mathrm{1}$${}^{\circ }$ at 00:00 UTC was used to initialize the meteorology at every forecast run.

Methodology for the evaluation of the simulations

The evaluation of the simulations is done in three stages: (a) an internal check of the data assimilation system; (b) evaluation of the analysis using as reference independent observations; (c) evaluation of a 5-day forecast with and without analysis initialization using as reference independent observations.

Regional domains and respective groups of AERONET stations used for validation purposes.

Regional domain (short name)	AERONET stations
Long Atlantic transport (LongAtl)	La_Parguera, White_Sands_HELSTF, Univ_of_Houston
Short Atlantic transport (ShortAtl)	Capo_Verde, Dakar, La_Laguna
Sub-Sahel (SubSahel)	IER_Cinzana, Banizoumbou, Ilorin, Agoufou
Sahara (Sahara)	Tamanrasset_INM
Extended Mediterranean (ExtMediter)	Saada, FORTH_CRETE, Lecce_University, Rome_Tor_Vergata
	Villefranche, Avignon, Evora, Barcelona, Granada
Middle East (MiddleEast)	SEDE_BOKER, Solar Village, Hamim
Central Asia (CenAsia)	–
East Asia (EastAsia)	–

Map of AERONET stations and of the different regional domains used for validation purposes. The regional domains are indicated with different colours: Long Atlantic transport (LongAtl) in blue, Short Atlantic transport (ShortAtl) in red, Sub-Sahel (SubSahel) in orange, Sahara (Sahara) in green, Extended Mediterranean (ExtMediter) in yellow, Middle East (MiddleEast) in pink, Central Asia (CenAsia) in granada, and East Asia (EastAsia) in cyan.

[Figure omitted. See PDF]

The consistency of the data assimilation system is checked through considerations of statistics of the ensemble, of simulation departures from assimilated observations, and of the temporal mean of assimilation increments. The ensemble mean and the coefficient of variation for the ensemble are calculated with and without data assimilation. The coefficient of variation is defined as the ratio of the standard deviation of the ensemble to the ensemble mean. Additionally, statistics for first-guess (FG) and analysis (AN) departures are calculated, where departures are defined as the difference between assimilated observations and simulations (first-guess or analysis), while mean increments are defined as the temporal mean of differences between analysis and first guess at the different time slots of the assimilation window.

Dust optical depth averaged for the month of May 2007 for the Control (top left), ENS-free-run (top right), DA-NRL (centre left), and DA-NRL-DB (centre right) experiments, and dust optical depth difference between the DA-NRL (bottom left), DA-NRL-DB (bottom right), and ENS-free-run experiment.

[Figure omitted. See PDF]

The evaluation of analysis and forecast with respect to independent observations are performed in terms of statistics of model field errors $e_i$ from observations, where $e_i=m_i-o_i$, with index $i$ indicating an instance of observation $o_i$ and where $m_i$ is the model field in observation space, bi-linearly interpolated at the observation location. We consider the root mean square error (RMSE), the mean error (BIAS), the standard deviation of the error (SD), the fraction gross error (FRGE), and the correlation coefficient (CORR) of the model AOD compared to either quality-assured (Level 2.0) AERONET or satellite retrievals. The FRGE $=\frac{\mathrm{2}}{n}{\sum }_{i=\mathrm{0}}^n\left|\frac{o_i-m_i}{o_i+m_i}\right|$ is added to the most widely used set of statistics for the error as it behaves symmetrically with respect to underestimation and overestimation without emphasizing the outliers, and is normalized to the sum of observation and simulation values. The SD of the error, though it can be derived from the other statistics, is also reported so as to make more explicit the changes in the bias-free mean square error and aid the interpretation of the evaluation results. The above set of evaluation statistics are calculated for measurements from individual ground-based stations, groups of stations, regional domains observed by satellite sensors, and globally.

For AERONET AOD measurements dust-dominated conditions are identified using the approach of as follows: AOD is classified as “Dust” AOD when the associated AE $<$ 0.75; we set “Dust” AOD to 0 when the associated AE $>$ 1.3; we identify a mixed aerosol type when the associated 0.75 $<$ AE $<$ 1.3. The latter values are excluded from the validation. We use the AERONET AOD value closest to the model time step and within a $\pm$30 min interval. For satellite AOD retrievals we use the set of satellite observations quality-controlled and filtered for dust-dominated conditions used in the assimilation step. We use these satellite observations to validate uniquely the forecast range following the assimilation window. We show the forecast evaluation statistics corresponding to measurements and simulations at 12:00 UTC only, so that they refer to an approximately equal number of pairs of observations and model simulated values at each forecast lead time that we are considering. Hence a smaller number of AERONET observations (at 12:00 UTC only) are used to verify the forecast compared to the ones used in the evaluation of the analysis.

We have identified eight regions of interest for the validation purposes in our study period, namely Long Atlantic transport (LongAtl), Short Atlantic transport (ShortAtl), Sub-Sahel (SubSahel), Sahara (Sahara), Extended Mediterranean (ExtMediter), Middle East (MiddleEast), Central Asia (CenAsia), and East Asia (EastAsia). These names do not necessary correspond to the conventional names of exact geographical locations but are meant to identify regional domains in a convenient way according to dust intrusions and to group observational stations. Most of the regional domains contain ground-based stations with a minimum number of observations during the study period (stations with fewer than 30 “Dust” observations are discarded), with the exception of Central and East Asia. The ground-based stations are listed in Table  and shown in the map in Fig.  together with regional domains used for the validation of the experiments against either ground-based or satellite observations.

Coefficient of variation for the month of May 2007 for the ENS-free-run (top), DA-NRL (centre), and DA-NRL-DB (bottom) experiments, when the ensemble is created perturbing the emitted mass vertical flux for each dust bin and the threshold on the friction velocity generating dust horizontal flux.

[Figure omitted. See PDF]

Coefficient of variation for the month of May 2007 for the ENS-free-run (left) and DA-NRL-DB (right) experiments, when the ensemble is created perturbing the emitted mass vertical flux for each dust bin.

[Figure omitted. See PDF]

Evaluation results

Ensemble, departure, and increment statistics

We compare here the dust fields in the Control, ENS-free-run, DA-NRL, and DA-NRL-DB experiments in terms of mean values and, when applicable, ensemble spread. Figure  shows the dust AOD values averaged over a month of the study period for the four above experiments, and the difference in AOD between the data assimilation experiments and the ENS-free-run. By visual inspection it can be noticed that the ensemble mean of the ENS-free-run experiment compares well with the Control experiment, which suggests that the ensemble perturbations are altering only to a small extent the model mean state as defined by a standard run. The analysis clearly shows conspicuous changes in the dust field compared to the Control experiment or the ENS-free-run. Figure  shows the coefficient of variation for AOD in the

Statistics of departures of first guess and analysis from assimilated observations, calculated for May to August 2007.

Experiment (departures)	Observations	BIAS	RMSE	CORR	FRGE	SD
DA-NRL (FG)	NRL	$-\mathrm{0.074}$	0.37	0.59	0.66	0.36
DA-NRL (AN)	NRL	$-\mathrm{0.118}$	0.27	0.75	0.54	0.24
DA-NRL-DB (FG)	NRL	$-\mathrm{0.160}$	0.35	0.58	0.70	0.31
DA-NRL-DB (AN)	NRL	$-\mathrm{0.169}$	0.29	0.72	0.61	0.24
DA-NRL-DB (FG)	DB	$-\mathrm{0.001}$	0.35	0.40	0.49	0.35
DA-NRL-DB (AN)	DB	$-\mathrm{0.075}$	0.23	0.64	0.35	0.22

ENS-free-run and the data assimilation experiments. Data assimilation clearly lowers the values of the coefficient of variation in the regions where observations are present, with values lower for the DA-NRL-DB than for the DA-NRL experiment, which indicates a reduction of the ensemble spread due to the assimilated observations. The high values of the coefficient of variation in the Southern Hemisphere, with or without data assimilation, are due to the perturbation of the dust sources present in the southern part of the globe. These values are not negligible due to differences among the ensemble members normalized to small dust AOD values. The ensemble of Fig.  (and Fig. ) is created by perturbing the emitted mass vertical flux for each dust bin and the threshold on the friction velocity generating dust horizontal flux. Creating the ensemble without perturbing the threshold on the friction velocity produces a reduced spread. See Fig.  for this second configuration of the ensemble with a coefficient of variation for the ENS-free-run in the left panel and for the experiment with data assimilation in the right panel. Perturbing the threshold on the friction velocity has an impact on the spread also outside source regions because, as explained earlier, the spin-up period for the ensemble ensures that perturbations applied at the sources propagate everywhere. Furthermore, this ensemble configuration better represents model uncertainty since the ratio of the prior total spread (the square root of the sum of the ensemble background variance and the observation error variance) to the prior RMSE (of the ensemble background against NRL MODIS and MODIS Deep Blue observations) is closer to $\mathrm{1}$ compared to when no perturbation is applied to the threshold on the friction velocity. It should be noted, however, that this chosen ensemble configuration under-represents uncertainty since it has a prior total spread smaller than the RMSE (ratio equal to 0.82). Other better perturbations should be tested for a future implementation since an underrepresentation of the background uncertainty might translate into giving a lower weight to the observations with respect to the background.

Binned scatter plots of the counts of the logarithm of assimilated observations and first guess (left plot) and analysis (right plot) for the DA-NRL experiment (top row) and DA-NRL-DB experiment (central and bottom rows), calculated for May to August 2007. A logarithmic scale is used for the counts.

[Figure omitted. See PDF]

Mean dust AOD analysis increments for May to August 2007 at 12:00 UTC for the DA-NRL experiment (left) and for the DA-NRL-DB experiment (right).

[Figure omitted. See PDF]

Time series of AOD values for May and August 2007 in La Parguera (top left), Dakar (top right), Ilorin (centre left), Tamanrasset INM (centre right), Lecce University (bottom left), and Hamim (bottom right) for the Control (blue), DA-NRL (green), and DA-NRL-DB (red) experiments, for MODIS AOD (NRL and DB; magenta circles), and for AERONET AOD (black triangles) in dust-dominated conditions. Analysis values are used for the data assimilation experiments.

[Figure omitted. See PDF]

Maps of validation statistics: BIAS, RMSE, CORR, and FRGE for the Control (left), DA-NRL (centre), and DA-NRL-DB (right) experiments calculated against AERONET AOD for a selection of stations providing observations during the study period (May to August 2007). Maps of the observation counts used for validation are shown in the bottom row.

[Figure omitted. See PDF]

We evaluate in the rest of this section the assimilation experiments in terms of statistics of the departures of the analysis and first guess from the assimilated satellite observations. Figure  shows for May to August 2007 first-guess dust AOD (in the left panels) and analysis dust AOD (in the right panels) versus observations for the DA-NRL and DA-NRL-DB experiments. The departure statistics with respect to the two sets of observations that we have assimilated are in Table . In both experiments a smaller scatter and a higher correlation coefficient for the analysis indicate that the assimilation improves the agreement with observations and hence a positive sanity check of the data assimilation system. The asymmetric behaviour of all the analysis scatter plots suggests that the system is more successful in correcting too high AOD values than correcting too low AOD values, which could be due to the fact that usually we have larger observation errors and a smaller ensemble spread for low AOD values. The BIAS is significantly smaller than the RMSE and the RMSE improves in the analysis over the forecast. The issue of a higher BIAS in the analysis departures compared to the first-guess departures has been identified in other assimilation systems (see , Sect. 4) and might be attributed to the fact that AOD is a positive definite variable, as this provides a deviation from the Gaussianity condition in the prior which is assumed in the analysis step. A solution to this problem worth investigating in the future would consist in applying a transformation of the state variables into new variables which present Gaussian features, a procedure known as Gaussian anamorphosis .

Figure  shows global maps of mean dust AOD analysis increments, i.e. the monthly averaged difference between analysis and short-term forecast respectively in the case in which only NRL MODIS AOD observations are assimilated and in the case in which also MODIS Deep Blue AOD observations are assimilated. Both experiments show non-zero systematic increments which are to be interpreted as systematic corrections that these sets of observations are making, in particular removing mass close to sources and, to a lesser extent, adding mass in the outflow. The spatial distribution of the increments highlights the role that MODIS Deep Blue observations play in particular over the Sahara dust sources. There are some regions where the two data assimilation experiments show opposite increments. This could be due to unresolved conflicting biases between the two types of MODIS retrievals.

Validation of the analysis

We perform in this section a validation of the dust fields simulated either with or without data assimilation through a comparison with observations from ground-based stations that have not been assimilated for May to August 2007. We calculate the statistics for individual stations and for groups of stations. Figure  shows the time series of dust AOD values for May to August 2007 for the Control experiment (blue), for the analysis of the DA-NRL (green) and of the DA-NRL-DB (red) experiment, and for AERONET observations in dust-dominated conditions (black) at six locations within the different regional domains of Fig. , which are in the proximity of dust sources (Tamanrasset in Algeria), affected by short-range dust transport (Dakar in Senegal, Ilorin in Nigeria, and Hamim in the United Arab Emirates), or affected by long-range dust transport in Europe (Lecce in Italy), and across the Atlantic (La Parguera in Puerto Rico). For reference, the MODIS AOD observations from the assimilated dataset (NRL and Deep Blue), which are at the closest distance and within a $\mathrm{2}$${}^{\circ }$ radius from the location of the AERONET station, are also included in the time series (magenta circles). Note, however, that these latter observations are not an independent reference for validation of the analyses, nor are they entirely representative of the observational constraint used to calculate the analysis in the given station location. The time series show an overestimation in the Control experiment of the optical depth near the sources, and to a smaller extent in the transport, which clearly suggests that the model tends to overestimate dust emissions. The current calibration for model version 1.0 has the shortcoming of accurately capturing long-range transport at the expense of an overestimation over the sources. This overestimation is reduced with data assimilation. By a first eyeball inspection, the AOD simulation variance is reduced by data assimilation and is more in accordance with the AOD observation variance.

BIAS, RMSE, CORR, FRGE, and SD for the Control experiment, for the experiment assimilating MODIS NRL observations (DA-NRL), and for the experiment assimilating MODIS NRL and MODIS Deep Blue observations (DA-NRL-DB) calculated against AERONET observations for all the stations in Fig. . The dust mean AOD for the observations used for validation during the experiment period is also reported.

[Figure omitted. See PDF]

Maps in Fig.  show results of validation statistics calculated for the full study period at each AERONET station for the three experiments performed. These maps allow us to appreciate the strongest features of the three simulations at individual AERONET stations and how those stations are representative of the regional domains that we have identified. The Control experiment shows that the strongest BIAS and highest RMSE are in the sub-Sahel region. The BIAS indicates that the model systematically overpredicts AOD in that region. The highest FRGEs are in the long transport over the Atlantic or Europe, as expected in areas of low AOD values. The correlation between model and observation values is in general lower near source areas than in outflow regions. This could be due to the too coarse model resolution not able to follow as well as the observations the dynamic of the dust field near source areas. The assimilation of MODIS NRL observations decreases some of the strongest biases, in particular in the dust outflow regions in the Sahel and over the African Atlantic coast, which is reflected in a reduced FRGE and RMSE, and is associated with improved correlation. The assimilation of the MODIS Deep Blue observations additionally to the NRL MODIS observations is of further benefit: it reduces the BIAS and RMSE downwind from the strongest dust sources of the Sahara. It is also relevant to notice that the additional assimilation of MODIS Deep Blue observations improves the correlation over the above areas and in the Arabian Peninsula.

BIAS, RMSE, CORR, FRGE, and SD for the Control experiment, the DA-NRL experiments, and the DA-NRL-DB experiment calculated against AERONET observations for groups of stations within the regional domains in Fig. . The dust mean AOD for the observations used for validation during the experiment period is also reported.

[Figure omitted. See PDF]

BIAS, RMSE, CORR, and FRGE for the forecast at $\mathrm{12}$, $\mathrm{36}$, $\mathrm{60}$, $\mathrm{84}$, and $\mathrm{108}$ h of the Control (blue) and AN-initialized (red) experiments, i.e. the experiment initialized with the DA-NRL-DB analysis, calculated against AERONET observations (left) and against global satellite retrievals, both NRL MODIS and MODIS Deep Blue, (right) filtered for dust-dominated conditions. The AERONET stations are the ones in Fig. . The dust mean AOD for the observations used to validate the $\mathrm{12}$ h forecast during the experiment period is also reported.

[Figure omitted. See PDF]

As the right panel of Fig.  but for the different regional domains of Fig. .

[Figure omitted. See PDF]

The chart plots for the validation statistics calculated for all the AERONET stations considered (hereafter called global statistics) and for stations grouped according to regional domains of interest are respectively in Figs.  and . Global statistics show that assimilation produces in general a better representation of dust concentrations in the atmosphere, and that the assimilation of Deep Blue retrievals has a positive impact over the assimilation of Dark Target retrievals only.

When considering the regional domains, the assimilation of NRL MODIS AOD has a positive impact on the quality of the analysis everywhere, with the only exception of a slight increase in RMSE in the Middle East region. This positive impact is more pronounced in the short Atlantic transport and in the sub-Sahel region. The additional assimilation of MODIS Deep Blue AOD has a considerable positive impact in the Sahara, sub-Sahel, and Middle East regions, and is neutral or slightly detrimental in the rest of the transport, in particular in the long-range Atlantic transport. The correlations for the global domain and for all the regional domains are highly statistically significant with the exception of the Sahara region (in the Control and DA-NRL experiments only), where the number of observations is smaller than other regional domains.

It should be noted, however, when interpreting the above statistics that the validation against AERONET observations may introduce significant errors when comparing a global model grid box against a point observation .

Validation of the forecast

We have validated the forecast up to 5 days ahead initialized at 00:00 UTC from either the Control experiment or an analysis (from DA-NRL-DB). We have calculated for May to August 2007 the errors for the forecast at 12, 36, 60, 84, and 108 h (hereafter indicated as FC $+$ 12, FC $+$ 36, FC $+$ 60, FC $+$ 84, and FC $+$ 108) with respect to either AERONET observations or satellite observations. As mentioned when describing our evaluation methodology, we use as a reference the set of satellite observations from the Dark Target and Deep Blue algorithm ingested in the assimilation step, i.e quality-controlled and filtered for dust-dominated conditions. They are used only to validate the forecast range following the assimilation window. As expected, all the validation statistics worsen with increased forecast step in both experiments (see Fig.  for global statistics). The impact of initializing the model with a dust analysis is positive on the first day. The analysis produces a better forecast in terms of BIAS and RMSE (and also SD of the error) up to FC $+$ 108, and a better correlation on the first day. The correlation is slightly lower from FC $+$ 36 onwards. The conclusions drawn by validating against AERONET or satellite observations are equivalent. Results calculated for regional domains (Fig. ) show that the Control experiment tends to overestimate AOD everywhere with the exception of Central and East Asia. This suggests an overestimation in particular of the Sahara emissions which is consistent with the bias found in the analysis and which is maintained during the forecast. The correlations for the global domain and for all the regional domains, at all forecast lead times, are highly statistically significant. Initializing the 00:00 UTC forecast with the DA-NRL-DB dust analysis reduces the overestimation compared to satellite retrievals on the first day of the forecast consistently with the improvement observed in the analysis in the previous section. However, this produces an underestimation of AOD in the long-range Atlantic transport during all the forecast lead times, which, because of the relatively small AOD values in that area, is reflected in particular in the FGRE. Although there is an overestimation of AOD, there is a better agreement of the temporal evolution in that region. The underestimation of AOD in the Atlantic transport might be due to too strong a deposition which affects in particular the long-range transport, and in the standard run is compensated for by an overestimation over the sources. As said earlier, a shortcoming of the current model calibration is to capture well the long-range transport at the expense of an overestimation over the sources, which data assimilation reduces. To identify the exact cause of it will require, however, further investigation together with a better adjustment of the current model parameters. With the exception of this underestimation of AOD across the Atlantic, all the error statistics and correlation coefficients are improved in the first day of the forecast in all the regional domains. The error of the analysis-initialized forecast is lower also in the rest of the forecast range (up to 5 days), though, after day 1, the correlation with satellite observations in some regions (SubSahel and ShortAtl) is lower for the analysis-initialized forecast than for a standard forecast. It is particularly relevant to notice that the dust forecast over the Sahara is improved for all the statistics and throughout the forecast range.

Conclusions

We have developed a data assimilation system for the NMMB-MONARCH model version 1.0, which considers dust only, while other aerosols are being implemented. We have coupled NMMB-MONARCH with an ensemble-based data assimilation technique known as LETKF. For this purpose we have created a forecast ensemble based on known uncertainties in the physical parametrizations of the mineral dust emission scheme. We have processed satellite aerosol optical depth retrievals for assimilation with a dust filter. Due to the presence of other aerosols in the selection of dust-dominated conditions, uncertainties might have been introduced into our assimilation process. It should be noted however that the identification of dust-dominated conditions is performed in this study as a proof of concept to demonstrate the potential of using data assimilation in NMMB-MONARCH, and will not be strictly necessary in a future model upgrade including all the major aerosol species. Still, efforts towards aerosol speciation could continue to be pursued when assimilating information about total aerosol optical properties. In this respect, operational centres currently rely merely on model background to distribute assimilation increments among the different aerosol species.

Assimilation experiments showed that aerosol optical depth retrieved with the Dark Target algorithm can help NMMB-MONARCH to better characterize atmospheric dust. This is particularly true for the analysis of the dust outflow in the Sahel region and over the African Atlantic coast. The additional assimilation of Deep Blue retrievals has a further positive impact in the analysis downwind from the strongest dust sources of the Sahara and in the Arabian Peninsula.

An analysis-initialized forecast performs better (lower forecast error and higher correlation) than a standard forecast everywhere on the first day of the forecast. The only exception to this is an underestimation of the forecast of AOD in the long-range Atlantic transport. The error of the analysis-initialized forecast is lower also in the rest of the forecast range (up to 5 days), though, after day 1, in sub-Sahel and short Atlantic transport the temporal evolution of dust is less in agreement with independent observations, compared to a standard forecast. Particularly relevant is the improved forecast over the Sahara throughout the forecast range thanks to the assimilation of Deep Blue retrievals over areas not easily covered by other observational datasets. To the best of our knowledge, this is the first study quantifying the benefit of assimilating MODIS Deep Blue from Collection 6 specifically for mineral dust simulations. This product is currently operationally assimilated by the UK Met Office, who consider only Deep Blue observations over desert, and by the European Centre for Medium-Range Weather Forecasts.

In our future implementation of the forecast ensemble, we plan to exploit spatial patterns of variation in model parameter uncertainty, for example source-dependent uncertainties, as well as uncertainties in the deposition term. A better representation of uncertainties in dust emission flux inherently will help the representation of uncertainties in other parts of the dust cycle. A recent study by shows that, for their system, combined meteorology and aerosol source ensembles are necessary to produce sufficient spread in outflow regions. Notwithstanding that their conclusion might be system-dependent, we will take into account their results in our future studies.

Copies of the code and data are readily available upon request from the corresponding authors.

The authors declare that they have no conflict of interest.

Acknowledgements

This work was funded by the SEV-2011-00067 grant of the Severo Ochoa Program awarded by the Spanish Government, the CGL-2013-46736-R grant of the Spanish Ministry of Economy and Competitiveness, and the ACTRIS Research Infrastructure Project of the European Union's Horizon 2020 research and innovation programme under grant agreement no. 654169. The authors thank all the Principal Investigators and their staff for establishing and maintaining the AERONET sites, NRL/University of North Dakota for the MODIS AOD L3 product, and the MODIS and OMI mission scientists and associated NASA personnel for the production of the AOD, AAI, and AE data used in this investigation. The authors thankfully acknowledge the computer resources at MareNostrum and the technical support provided by the Barcelona Supercomputing Center (RES-AECT-2015-1-0007). They also thank Francesco Benincasa for his technical support. Carlos Pérez García-Pando acknowledges long-term support from the AXA Research Fund, as well as the support received through the Ramón y Cajal programme (grant RYC-2015-18690) of the Spanish Ministry of Economy and Competitiveness. Comments from two anonymous reviewers are gratefully acknowledged. Edited by: A. Lauer Reviewed by: two anonymous referees