Data-driven aeolian dust emission scheme for climate modelling evaluated with EMAC 2.55.2

Aeolian dust has significant impacts on climate, public health, infrastructure and ecosystems. Assessing dust concentrations and the impacts is challenging because the emissions depend on many environmental factors and can vary greatly with meteorological conditions. We present a data-driven aeolian dust scheme that combines machine learning components and physical equations to predict atmospheric dust concentrations and quantify the sources. The numerical scheme was trained to reproduce dust aerosol optical depth retrievals by the Infrared Atmospheric Sounding Interferometer on board the MetOp-A satellite. The input parameters included meteorological variables from the fifth-generation atmospheric reanalysis of the European Centre for Medium-Range Weather Forecasts. The trained dust scheme can be applied as an emission submodel to be used in climate and Earth system models, which is reproducibly derived from observational data so that a priori assumptions and manual parameter tuning can be largely avoided. We compared the trained emission submodel to a state-of-the-art emission parameterisation, showing that it substantially improves the representation of aeolian dust in the global atmospheric chemistry–climate model EMAC.

1 Introduction

Aeolian dust is one of the most abundant aerosol components worldwide and substantially affects the Earth system in many ways. In contrast to sea salt, the only other component contributing a comparable fraction to the total aerosol loading, aeolian dust is emitted over land, where accordingly the highest dust loads occur and severely interfere with human health and activities. The generally high and often exceptional particulate matter concentrations in the vicinity of dust sources both cause acute health problems (Goudie, 2014) and contribute to the long-term exposure with associated health risks and excess mortality (Lelieveld et al., 2019b, a). Reduced visibility during dust events can interrupt road and air traffic, cause accidents, and reduce solar electricity production, the latter aggravated by deposited dust, resulting in high economic costs (Middleton et al., 2019). Since the particles interact with radiation and clouds, atmospheric dust significantly affects weather and climate. Deposited on snow and ice, mineral dust reduces the surface albedo and accelerates glacier melting (Di Mauro et al., 2019; Francis et al., 2022). On the other hand, deposited dust particles are an important source of mineral nutrients both on land and in the ocean (Bristow et al., 2010). Although aeolian dust predominantly originates from natural sources, it can be affected by anthropogenic climate change and land use (Klingmüller et al., 2016); moreover, atmospheric dust interacts with anthropogenic air pollution (Klingmüller et al., 2019, 2020). Consequently, aeolian dust is not a purely natural phenomenon.

These are only some of the diverse aspects which constitute the great relevance of aeolian dust in Earth science and drive the demand for accurate representations in atmospheric models. However, meeting this demand is challenging, as can be seen, for example, from the large uncertainty range of model estimates for the global annual dust emission (Huneeus et al., 2011). Aeolian dust emissions depend on many environmental factors. Considering that some of those have been identified only recently (Rodriguez Caballero et al., 2022), presumably not all relevant factors and their impacts are known or well quantified. Temporally changing environmental factors have to be considered, which can vary within seconds (wind and turbulence), seasons (vegetation), years (land cover) and geological timescales (soil composition). Quantitative observations of dust emissions are difficult, and thus a direct comparison of modelled emission fluxes with observations, especially on global scale, is not feasible. Therefore, dust emission models can only be validated indirectly based on atmospheric dust concentration observations or scarce dust deposition measurements, which requires that atmospheric processing and transport are taken into account.

We present a novel approach to address these challenges combining two developments that have become main drivers of Earth system science, the increasing availability of satellite observations and the advancing machine learning technology (Reichstein et al., 2019; Bauer et al., 2021; Eyring et al., 2021; Irrgang et al., 2021). We used the resulting hybrid dust emission scheme to substantially improve the representation of aeolian dust in the global atmospheric chemistry–climate model EMAC (Jöckel et al., 2006), demonstrating a valuable synergy of machine learning and physical process-based modelling.

This article is structured as follows. The data sets are described in Sect. 2. The architecture of the data-driven dust model is presented in Sect. 3. Its training and evaluation are described in Sects. 4 and 5, respectively. Section 6 documents the EMAC setup used in Sect. 7 to demonstrate the benefits of the newly derived dust emissions. Conclusions are drawn in Sect. 8.

2 Data

Our selection of input variables follows that of common online dust emissions schemes (e.g. Astitha et al., 2012; Klingmüller et al., 2018). Data for several input variables were taken from the fifth-generation European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalysis (ERA5), including the single-level variables of surface friction velocity, total precipitation, snow depth, volumetric soil water in the topmost soil layer, leaf area index (LAI) of low vegetation, LAI of high vegetation and surface geopotential (Hersbach et al., 2018b) and the height dependent variables of the u-component of wind, v-component of wind and vertical velocity (Hersbach et al., 2018a). The latter variables where considered at four pressure levels, 650, 750, 850 and 950 hPa, to represent the atmosphere up to about 4 km altitude. Mineral dust can reach higher altitudes, but since most of the dust mass remains within these layers, this approximation is a reasonable compromise between a realistic representation of dust transport and an acceptable computational burden. For all time-dependent variables, hourly data were considered. Horizontally, we re-gridded the data to a regular longitude–latitude grid at 1^∘ resolution.

The ERA5 data were complemented by the static soil clay fraction distribution from the SOILPOP30 database (Nickovic, 2011) and the annual Moderate-resolution Imaging Spectroradiometer (MODIS) collection 6 International Geosphere–Biosphere Programme (IGBP) land surface class fractions (Friedl and Sulla-Menashe, 2015), and both of these are re-gridded to the same 1^∘ grid.

As the output variable representing the atmospheric dust distribution we consider the dust aerosol optical depth (DAOD), which is closely related to the vertically integrated amount of dust and can be retrieved remotely from satellites so that data sets with consistent global coverage are available. We used version 8 of the 10 µm DAOD product developed at the Université libre de Bruxelles (ULB) based on observations by the Infrared Atmospheric Sounding Interferometer (IASI) on board the MetOp-A satellite (Clarisse et al., 2019). In the infrared, coarse particles scatter radiation more efficiently than submicrometre particles. Therefore, observations at a wavelength of 10 µm are most sensitive to mineral dust and sea salt and are affected little by other aerosols, resulting in a particularly reliable DAOD retrieval. We considered all years that are completely covered by the data set, in total 12 years from 2008 to 2019. Years with incomplete coverage and the continuation of the product based on retrievals from the MetOp-B satellite were omitted. The level 2 data were split into hourly data by rounding the observation times to the nearest full hour and gridded to the 1^∘ grid of the input data by averaging all retrievals covered by each grid cell. Grid cell values below the 0.05 % percentile and above the 99.95 % percentile were omitted.

Towards the poles, fewer DAOD retrievals are obtained, the grid cell geometry increasingly deviates from a rectangular shape that we assume in our dust scheme and the Courant–Friedrichs–Lewy criterion gains importance, while only a small fraction of the global dust is transported to or emitted from these regions. Therefore, we restricted the domain of all data sets to latitudes between 60^∘ S and 60^∘ N.

We assume the aerosol extinction contributing to the DAOD to be proportional to the dust concentration. Accordingly, the total amount of dust over each grid cell is proportional to the DAOD and the horizontal area of the grid cell. By multiplication with the latitude-dependent cell area, we converted the DAOD to a quantity proportional to the total amount of dust over each grid cell, which after normalisation corresponds to the vertical integral of the dust representation in our model.

Surface friction velocity, soil moisture, snow depth, LAI, geopotential and clay fraction were normalised by subtracting their mean and division by their standard deviation. Dust amount and precipitation were only divided by their standard deviation to preserve their sign. This normalisation maps all input variables and associated gradients into a similar and relatively small range to support the model parameter optimisation routine during the training process.

The 12 years of DAOD data and the corresponding input data were split into a validation period (2008), a test period (2009) and a training period (2010 to 2019). While we used the validation period to monitor the training procedure and validate the model development, the test period was exclusively used to evaluate the final trained model. Only the training period was used for model training. The order of the periods was chosen such that reliable emission data for EMAC was readily available for the test period, but at the same time the training period could be easily extended in the course of the model development. Moreover, the validation period before the test period serves as spin-up period of the EMAC simulations.

3 The data-driven dust model

The data-driven dust model combines multiple modules representing different processes, including dust emission, transport and removal. Figure 1 provides an outline of the model architecture. The processes are applied iteratively once every hour, synchronous with the hourly time step of the input and output data. Consistent with the input wind fields, the instantaneous dust concentration is implemented as three-dimensional array (d_ijk), where $i=\mathrm{1},\mathrm{\dots },n_{\text{long}}$, $j=\mathrm{1},\mathrm{\dots },n_{\text{lat}}$ and $k=\mathrm{1},\mathrm{\dots },n_{\text{lev}}$, corresponding to a grid with 1^∘ horizontal spacing (n_long=360, n_lat=120) and n_lev=4 vertical layers. All $n_{\text{long}}\times n_{\text{lat}}\times n_{\text{lev}}$ array elements d_ijk are initialised with 0. We used n_t=168 iterations, corresponding to a time interval of 168 h, i.e. 7 d.

Figure 1 Outline of the dust model. It combines multiple modules representing different processes and effects (small grey boxes), each with trainable parameters. The dust emissions are represented by three submodels (large grey box).

Figure 2 Predicted vs. observed hourly grid cell DAOD values during the test period. The colours represent the density of data points. Correlation coefficient (r), root-mean-square error (RMSE), mean absolute error (MAE) and the number of data pairs (n) are indicated.

The centrepiece of the model is the emission module, which itself consists of three components. At each horizontal grid cell (i, j), an interim emission flux s_ij is calculated by a stack of densely connected neural layers using the six input variables of surface friction velocity, LAI of low vegetation, LAI of high vegetation, soil moisture, snow cover and clay fraction. This generic deep neural network architecture is computationally efficient but still very powerful for various classification and, as in our application, regression problems.

The stack consists of two hidden layers of 64 units each with leaky rectified linear unit (ReLU) activations, $\text{LeakyReLU}\left(x\right)=max(\mathrm{0},x)+\mathrm{0.01}min(\mathrm{0},x)$, and a final output unit with softplus activation, $\text{Softplus}\left(x\right)=\ln (\mathrm{1}+\exp (x\left)\right)$, to ensure positive output. The non-zero gradient of the leaky ReLU function at negative values avoids vanishing gradients that can hinder the model optimisation process when using the ordinary ReLU activation function, which is zero at negative values. Within the layers, each unit calculates the weighted sum of its input variables, where the weights are trainable parameters, adds a trainable bias and applies the activation function before passing the result to the next layer. In the densely connected stack, all units in one layer are connected to all units in the next layer. The emission module is translation invariant in both space and time because identical trainable weights and biases are used for all locations (i,j) and time steps. Therefore, spatial and temporal variations of the modelled emissions are based solely on changing input parameters.

The effect of the topography is considered by multiplying the interim emission flux with a topography-dependent factor β_topo,ij, similarly to the use of topographic source functions in existing dust emission schemes (Ginoux et al., 2001; Klingmüller et al., 2018),

$$\begin{array}{}\text{(1)}& s_{ij}\leftarrow {\mathit{\beta }}_{\text{topo},ij}{s}_{ij}.\end{array}$$

Unlike the intermediate emission flux, which depends only on the local input parameter values, the topography effect considers the surrounding terrain to account for, e.g. the accumulation of sediments in valleys and depressions. In deep learning, such dependency on surrounding values is commonly addressed by convolutional layers that consider all values from grid cells defined by the convolution kernel and apply trainable weights and biases that are translation invariant. The topography factor β_topo,ij is computed from the geopotential height field by a stack of convolutional layers with three by three grid cell kernels. The stack consists of three hidden layers with four output channels each and leaky ReLU activations and a one channel output layer with sigmoid activation, $\text{sig}\left(x\right)=(\mathrm{1}+\exp (-x){)}^{-\mathrm{1}}$, to limit the output to the interval (0, 1). This architecture computes the topography factor for each grid cell based on the geopotential in the surrounding nine by nine grid cells, i.e. a 9^∘ by 9^∘ region, similar to the 10^∘ by 10^∘ region used by Ginoux et al. (2001).

The size of the convolutional network and the densely connected neural network, i.e. the number of layers and number of units or channels per layer, were adjusted based on experiments with different settings. The model should be large enough to provide sufficient capacity to learn the relevant relationships but not larger than necessary to avoid overfitting and reduce computational requirements. Model performance is not very sensitive to limited changes in model size but could be improved by systematic optimisation of the latter.

The land cover effect is also considered by multiplying the interim emission flux with a corresponding factor, generalising the concept of the land-cover-class-based emission mask used by Klingmüller et al. (2018) for example. The land cover factor is represented by a single linear layer with 16 inputs representing the 16 IGBP land cover classes. The layer multiplies the fraction f_ijl of grid cell (i,j) that is covered by class l with a trainable, class-specific emission efficiency β_lc,l ranging from 0 to 1 and sums over all classes to obtain the average emission efficiency of the grid cell ${\mathit{\beta }}_{\text{lc},ij}={\sum }_{l=\mathrm{1}}^{\mathrm{16}}{f}_{ijl}{\mathit{\beta }}_{\text{lc},l}$. The interim emission flux is multiplied by this emission efficiency to obtain the final emission flux,

$$\begin{array}{}\text{(2)}& s_{ij}\leftarrow {\mathit{\beta }}_{\text{lc},ij}{s}_{ij}.\end{array}$$

According to this emission flux, the dust concentration in the lowest vertical layer (k=n_lev) of the dust array is increased,

$$\begin{array}{}\text{(3)}& d_{ij{n}_{\text{lev}}}\leftarrow d_{ij{n}_{\text{lev}}}+s_{ij}{A}_{ij},\end{array}$$

where A_ij is the grid cell area, before the dust array is processed further by the transport module.

The transport module successively applies vertical diffusion, vertical wind transport, north–south transport and east–west transport. The vertical diffusion considers that some fraction of the dust is exchanged vertically between grid cells independent of the large-scale vertical velocity of the air masses, e.g. by convection or sedimentation, and is parameterised by two trainable parameters, the total fraction of dust leaving the grid cell and the fraction thereof moving upwards. The transport by wind in each of the three spatial dimensions is calculated by considering the overlap of the grid cells translated according to the wind speed and the time step length of 1 h with the neighbouring cells and has no trainable parameters. The wind speed components are limited so that the Courant number C≤1. In the vertical transport calculation, dry deposition is considered by allowing a fraction of the dust transported onto the surface to remain on the surface, the rest is retained in the lowest dust layer. The fraction represented by the latter is defined by a trainable parameter. However, during training this parameter typically quickly converges to 1, indicating that the current implementation of our approach is not sensitive enough to identify dry deposition. The reason is likely the limited atmospheric residence time discussed below, representing an unrealistic additional removal process. Moreover, a fraction of the deposition that always occurs co-located with emissions might be represented by a reduced emission flux.

We consider wet deposition as the dominant removal process. Every hourly time step the dust concentration in all vertical layers is reduced by a fraction proportional to the local total precipitation p_ij,

$$\begin{array}{}\text{(4)}& d_{ijk}\leftarrow (\mathrm{1}-{\mathit{\beta }}_{\mathrm{p}}{p}_{ij})d_{ijk},\end{array}$$

where the proportionality factor β_p is a trainable parameter. This approximation considers the correlation between precipitation and dust removal as the dominant statistical relationship with the strongest influence on dust source parameters during training. More detailed parameterisations could improve future versions of the model. In addition, the finite number of iterations considered by the model implies a corresponding maximum atmospheric residence time for the dust. Here, the use of n_t=168 iterations corresponds to a maximum residence time of 1 week. In reality, dust particles can reside much longer in the atmosphere, but because the contribution of those particles to the global DAOD distribution is comparatively small, most DAOD features and the dominant dust sources can be identified regardless of this restriction. However, a noticeable bias of the DAOD predictions is expected in regions remote from dust sources.

The three-dimensional dust field is summed along the vertical dimension to obtain the total amount of dust in the vertical column. Finally, de-normalisation and division by the cell area yield the DAOD prediction.

Since by design the dust emission flux is positive and neither the transport nor the removal processes can alter the sign of the dust concentration in each grid cell, the amount of dust and the DAOD predictions are guaranteed to be positive, which is desired for a realistic model. On the other hand, the noise range of the IASI retrievals includes negative values that crucially contribute to a realistic overall DAOD level. Even though these retrievals cannot be reproduced by the model, they are included in the training process to avoid a positive bias in the DAOD predictions. In combination with the expected negative bias in remote regions, this a priori limits the agreement of predicted and observed DAOD. However, for our main purpose, the extraction of enhanced dust emissions for climate simulations, it is sufficient, as demonstrated in the following.

We implemented the model using the PyTorch machine learning framework (Paszke et al., 2019). This framework is very flexible and allows the seamless combination of the different modules described above, including standard deep learning models such as densely connected neural layers and convolutional networks, classical regression models, and iterative physical processing as used in the transport module, into a single hybrid model. The autograd engine, one of PyTorch's core features, automatically differentiates this complex model with respect to the trainable parameters, and the resulting gradients are used by the optimisation routines during model training.

Figure 3 Correlation coefficient of the observed and predicted DAOD time series for each grid cell during the test period.

Figure 4 Daily correlation coefficient of the observed and predicted DAOD during the test period.

4 Training

During training, the predicted amount of dust in each vertical column is compared to the observed DAOD (the latter normalised and weighted by the grid cell area as described above) by means of a loss function.

In addition to the discrepancy between predictions and observations, the training procedure should also consider reasonable a priori assumptions: The surface friction velocity initiates and sustains the horizontal flux of saltating particles and is therefore the driving force of mineral dust emissions. Consequently, the dust emission flux generally increases with increasing surface friction velocity. In contrast, vegetation, soil moisture and snow cover reduce soil erodibility and saltation and thus the emission flux. To incorporate this information into the emission module we imposed a knowledge-based regularisation by penalising violations of these rules similar to Gupta et al. (2019). For this purpose, before each optimisation step during training, we computed the partial derivatives of the emission module with respect to the corresponding input variables for n_s=10⁶ randomly sampled input parameters. The values consistent with the rules were set to zero, i.e. positive derivatives with respect to the surface friction velocity and negative derivatives with respect to the high and low vegetation LAI, soil moisture and snow depth. The mean of the absolute value of the thus-adjusted derivatives, divided by the mean of the n_s=10⁶ emission values, was added to the mean squared error (MSE) loss function f_MSE to obtain the total loss function,

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{f}_{\text{loss}}=& f_{\text{MSE}}+\left({\displaystyle \frac{\mathrm{1}}{\mathrm{5}{n}_{\mathrm{s}}}}\sum _{l=\mathrm{1}}^{n_{\mathrm{s}}}\sum _{m=\mathrm{1}}^{\mathrm{5}}\text{ReLU}\left({\mathit{\sigma }}_m{\displaystyle \frac{\partial s_l}{\partial x_m}}\right)\right)\\ & \times {\left({\displaystyle \frac{\mathrm{1}}{n_{\mathrm{s}}}}\sum _{l=\mathrm{1}}^{n_{\mathrm{s}}}{s}_l\right)}^{-\mathrm{1}},\end{array}\end{array}$$

where (x_m) are the input variables and ${\mathit{\sigma }}_m=-\mathrm{1}$ if x_m is the surface friction velocity, otherwise σ_m=1. The total loss was then minimised using the Adam optimiser (Kingma and Ba, 2014).

To generalise the model and prevent overfitting during training, we regularised the stack of densely connected layers in the emission module by applying a 5 % dropout rate to the hidden layers (Hinton et al., 2012). The combination of machine learning components and numerical implementations of physical processes such as atmospheric transport in a single model adds a novel aspect to the overfitting problem. For pure machine learning models, analysing the validation data set is sufficient to detect overfitting. In contrast, when combined with numerics, the machine learning algorithm can reduce the validation loss further by adapting to the numerical approximations (e.g. maximum transport period and discretisation), which is then not necessarily based on a more realistic representation of the underlying process. To avoid this type of overfitting, we trained the model using the whole training data set only once and did not repeat the training for multiple epochs. Tests using more epochs showed that the training loss can indeed be reduced further, meaning that the DAOD predictions become more realistic; however, the resulting emissions did not clearly improve the EMAC simulations.

5 Trained model evaluation

Figure 2 shows the hourly grid cell DAOD values predicted by the trained model during the test period (year 2009) vs. the corresponding IASI retrievals. In the range typically observed over the dust belt between West Africa and East Asia (10 µm DAOD ≳ 0.05), the predictions are realistic. Towards smaller values, the model seems to underestimate the observations. This is partly attributed to the expected underestimation in remote areas due to the aforementioned limitation of the transport time. The large scatter of the observed values around DAOD = 0 corresponds to a maximum of the density distribution at about 0.01 on the logarithmic scale. Accordingly, observations in this range are consistent with 0 and the low model predictions. The overall good performance of the model is reflected by a high Pearson correlation coefficient of 0.502, calculated using all 14 371 211 prediction–observation pairs for the test data set.

Figure 5 Annual mean 10 µm DAOD during the year 2009 observed by IASI (a), predicted by the trained dust model (b), and simulated by EMAC with reference (c) and trained (d) mineral emissions.

Figure 6 The importance of the individual input variables, estimated by their impact on the model skill score S. For each variable, the impact was obtained by comparison with a model evaluation replacing the variable with random or, in case of the topography, constant values.

The temporal correlation coefficients of the observed and predicted hourly DAOD values within each grid cell (i,j), $r_{ij}=r({\text{DAOD}}_{\text{obs},ij},{\text{DAOD}}_{\text{pred},ij})$, are typically greater than 0.5 over the regions affected by desert dust (Fig. 3). In regions remote from dust sources, where mineral dust is of minor relevance, the low, randomly varying dust concentrations result in correlation coefficients close to zero, e.g. over large parts of the Pacific and Southern oceans. The negative correlation coefficients over some regions, notably west of Australia over the southern Indian Ocean, result from the transport time limit in our model. They occur when enhanced DAOD values are caused by dust reaching the region after more than 1 week of transport over non-emitting regions including the ocean, while the model predicts clean air with low DAOD values. The daily correlation coefficient, $r\left(t\right)=r\left({\text{DAOD}}_{\text{obs}}\right(t),{\text{DAOD}}_{\text{pred}}(t\left)\right)$, quantifying the similarity of the spatial pattern of predictions and observations on day t, also reaches or exceeds 0.5 during the dust season in Northern Hemispheric spring and summer (Fig. 4). During winter, the decline in dust activity and thus the absence of distinct spatial patterns results in lower correlation coefficients.

Figure 7 The modelled emission flux as a function of surface friction velocity (a, without soil moisture) and soil moisture (b, with surface friction velocity 0.7 m s⁻¹). The flux was computed assuming the absence of vegetation and snow and an average clay fraction (20.6 %). The topography and land cover factors were set to unity and the normalisation described in Sect. 6 was applied.

The distribution of the predicted annual mean DAOD (Fig. 5b) agrees well with the IASI result (Fig. 5a). The most notable discrepancy is an overestimation by the predictions over the Sahara. Here, the mean values of the IASI retrievals are significantly affected by small or even negative values between the major dust outbreaks, where the model usually predicts small but non-zero values which accumulate to a larger annual average.

We conclude that the model generally reproduces the main features of the global dust distribution. This makes it suitable for the application in climate models, as discussed in the next sections. Additionally, valuable information can be obtained directly from analysing the model.

Table 1 Fractional contribution of major source regions to the global mineral dust emissions.

To study the relative contribution of different source regions to the global dust emissions, we compute hourly dust emission fluxes by evaluating the dust emission submodel (light grey box in Fig. 1) from the trained model alone. Accumulating the emission fluxes from the 9 major dust source regions considered by Kok et al. (2021) over the test period (year 2009) yields the fractional contributions to the global dust emissions shown in Table 1. These contributions are well within the range of the models considered by Kok et al. (2021), with a comparably large contribution from the Middle East and Central Asia.

To estimate the importance of the individual input variables, we made use of the skill score S defined by Taylor (2001),

$$\begin{array}{}\text{(6)}& S={\displaystyle \frac{\mathrm{4}(\mathrm{1}+r{)}^{\mathrm{4}}}{({\mathit{\sigma }}_{\mathrm{1}}/{\mathit{\sigma }}_{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{2}}/{\mathit{\sigma }}_{\mathrm{1}}{)}^{\mathrm{2}}(\mathrm{1}+r_{\mathrm{0}}{)}^{\mathrm{4}}}},\end{array}$$

where r is the correlation coefficient and σ₁ and σ₂ are the standard deviations of modelled and observed values, respectively. As maximum attainable correlation coefficient we simply used r₀=1 since we focus on relative changes of the skill score (a more accurate estimate r₀<1 would result in higher skill scores).

Figure 8 Mean aeolian dust emission flux during the test period.

Figure 9 Annual mean 550 nm AOD during the year 2009 observed by MODIS (a) and simulated by EMAC with reference (b) and trained (c) mineral dust emissions.

We computed the DAOD during the test period multiple times, each time replacing the values of one input variable with random values sampled from the probability distribution of the same variable, or, in the case of topography, using its mean value (695 m) as a constant input (random values would correspond to an unrealistic terrain). Subtracting the resulting skill scores from those obtained using the full input data yields the impact ΔS of each variable, which indicates the importance of the variable to the model predictions. Figure 6 identifies soil moisture and surface friction velocity as important input parameters. Due to the slightly different methodology, the exact value obtained for the topography should not be directly compared with the other values, but the use of topographic data clearly improves the model predictions. In contrast, the effect of snow is small because conditions where snow suppresses emissions that are not already suppressed by other input variables are rare. We assume that the importance of the variables for the model predictions is to some extent representative of their real importance for dust emissions. This assumption neglects the interdependence of the input variables as well as the effect of model approximations and needs to be verified by detailed analysis in future studies. Nevertheless, the importance of soil moisture in our model is remarkable and consistent with the conclusions of Klingmüller et al. (2016).

Similarly, the dependence of modelled emissions on input parameters should indicate realistic relationships. As an example, the dependence of emissions on the two local input variables with the greatest influence on the skill score, surface friction velocity and soil moisture, was assessed (Fig. 7). Under ideal conditions, i.e. without soil moisture, vegetation and snow, emissions increase only slowly with surface friction velocity until they increase more rapidly from about 0.2 m s⁻¹. As soil moisture increases, emissions decrease until they vanish at a volumetric soil moisture of about 0.3. Like the parameter importance, this analysis will be addressed in detail in a separate study.

6 EMAC setup

The output of the emission submodel corresponds to the change of the DAOD (weighted by the grid cell area and the standard deviation of the DAOD). Avoiding assumptions on the size distribution and optical properties of the mineral dust particles and consistent with our simplifying assumption that dust extinction and concentration are proportional, we converted the hourly emission flux output to a mass concentration flux in $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ using a single tuning parameter chosen such that the global annual emission for 2008 amounts to 4.3 Gt, which is the upper limit of the range reported by Huneeus et al. (2011). This choice represents a relatively high global total emission. However, since our model does not use any a priori emission mask, regions with locally very low emissions also contribute to the global total emissions, which significantly enhances the value due to the large area of those regions. Our choice for the total emission value is supported by the realistic DAOD levels we obtain over the dust belt. The mean of the resulting emission flux during the test period in 2009 is shown in Fig. 8.

Figure 10 Location of the AERONET stations used for the model evaluation.

The resulting dust emission flux data were used in a validation simulation with the ECHAM/MESSy atmospheric chemistry–climate model (EMAC) (Jöckel et al., 2006). We employed a recent EMAC version based on release 2.55.2. A list of the activated EMAC submodels including references is provided in Table A1 in Appendix A.

The total dust emission flux was distributed over the accumulation and coarse modes consistently with the online dust emission scheme, where 94.7 % of the dust is emitted into the coarse mode. The modes are represented by log-normal distributions with fixed geometric standard deviations (σ_g=1.59 for the accumulation and σ_g=2 for the coarse mode). The count median diameter of each mode can vary with a fixed threshold between accumulation and course mode at 1.4 µm. Accordingly, the percentage of fine particles in the emissions is comparable to the 4.3 % (95 % confidence interval 3.5 % to 5.7 %) found by Kok et al. (2017) for particles smaller than 2 µm. There is no strict upper limit for the size of coarse particles, but super-coarse particles with diameters exceeding 20 µm are only represented by the upper tail of the log-normal coarse mode.

The simulation covered 3 years from 2008 to 2010, where the year 2008 is used as the model spin-up period. The results for the years 2009 and 2010 are used to evaluate the performance of the trained dust emissions during test and training periods, respectively.

Horizontally, the setup used a Gaussian T63 grid with approximately 1.9^∘ latitude–longitude spacing. The simulation assimilated observational data by nudging temperature, vorticity and divergence above the boundary layer to meteorological reanalysis data of the ECMWF and by using the sea surface temperature from the same data set.

For comparison, we performed a reference simulation with the same settings but using the online dust emission scheme of Klingmüller et al. (2018) instead of the trained emissions.

7 EMAC evaluation

The aerosol optical depth (AOD) at 10 µm wavelength is dominated by coarse aerosol particles, which in turn are dominated by mineral dust and sea salt particles. Accordingly, we obtain the EMAC DAOD by weighting the 10 µm AOD with the volume fraction of dust in the local dust (DU, 2.65 g cm⁻³) and sea salt (SS, 2.17 g cm⁻³) mixture, DAOD(10 µm) = AOD(10 µm) $V_{\text{DU}}/(V_{\text{DU}}$ + V_SS).

Figure 11 Skill score S and correlation coefficient r comparing the daily AOD simulated by EMAC with AERONET observations. The AERONET data were interpolated to the 550 nm wavelength of the model output using the Ångström exponent. The red, green and grey bars depict the differences between the results for the trained emissions and for the reference emissions, with green bars indicating that the results for the trained emissions agree more closely with the measurements by at least a standard deviation (σ). The corresponding error intervals are indicated by darker colours.

The results for the annual mean 10 µm DAOD for the year 2009 from the reference and validation simulations are shown in Fig. 5c and d. Both simulations yield a DAOD distribution similar to that obtained from IASI. Note that the model results are obtained from daily average values and thus consider data for all model time steps of the year, whereas the IASI result is based on the instantaneous values, which are retrieved at most twice a day and include a substantial amount of missing values; therefore, the IASI distribution is more noisy and has a larger year-on-year variability.

The fact that the trained emissions reproduce the high DAOD values that are expected over the dust belt, but low levels elsewhere, is quite remarkable as the model architecture treats all regions worldwide equally and there is no a priori selection of dust source regions involved. Moreover, the spatial pattern of the DAOD resulting from the trained emissions agrees more closely with the IASI observation than the reference result. The correlation coefficient of the modelled annual grid cell DAOD values and the corresponding IASI values (re-sampled to the model grid) is 0.82 for the reference simulation, but 0.88 for the trained emissions.

The annual mean AOD at 550 nm from the reference and validation simulations are compared to the corresponding MODIS observations in Fig. 9. Again the result based on the trained emissions is closer to the observations than the reference results. The correlation coefficient of the modelled annual grid cell AOD values and the corresponding MODIS values (re-sampled to the model grid) is 0.72 for the reference simulation and 0.73 for the trained emissions.

To complement the satellite observations, we considered ground-based observations by Aerosol Robotic Network (AERONET) stations in the seven dust-affected regions defined by Klingmüller et al. (2018) (Fig. 10). We quantified the degree of agreement of AERONET and EMAC daily AOD values using the Pearson correlation coefficient and the skill score S. The AERONET data were interpolated to the 550 nm wavelength of the EMAC output using the Ångström exponent.

Figure 11 reveals that the validation simulation using the trained emissions achieves substantially higher skill scores and correlation coefficients over Africa and the Middle East than the reference simulation. Over regions where the atmospheric aerosol is less dominated by aeolian dust, including Asia, the Americas and Australia, the effect of the new emissions is smaller but generally enhances both the skill scores and the correlation coefficients.

Figure 12 compares the modelled annual mineral dust deposition to observations at different sites worldwide (Huneeus et al., 2011; Checa-Garcia et al., 2021). Both the correlation coefficient r and the skill score S are higher for the validation simulation, whereas the root-mean-square error (RMSE) and the mean average error (MEA) are smaller. All metrics indicate that the data-driven emission scheme used in the validation simulation results in a more realistic global distribution of the mineral dust deposition than in the reference simulation.

Figure 12 Comparison of the mineral dust deposition simulated by EMAC and the observations from the AeroCom dust benchmark data set (Huneeus et al., 2011). Correlation coefficient (r), root-mean-square error (RMSE), mean absolute error (MAE) and skill score (S) are provided for both the reference simulation and the validation simulation using the trained dust emissions.

8 Conclusions

We presented a trainable aeolian dust model that successfully learnt from satellite observations to compute realistic global dust distributions and emissions based on the environmental conditions. On the one hand, the model architecture incorporates few a priori assumptions about the nature of dust sources, and on the other hand we explicitly modelled the better known components of the dust cycle, in particular the atmospheric transport.

We focussed on the emission submodel of the trained model, representing a novel dust emission scheme that was derived independently of existing methods. Using the EMAC model as an example, we demonstrated its superior performance when applied within an atmospheric chemistry–climate model.

Future developments will include model improvements by adjusting hyperparameters and resolution, addressing characteristics of the DAOD data set through preprocessing and model adjustments, considering additional input data (e.g. on biological soil crusts,  Rodriguez Caballero et al., 2022), and utilising related machine learning models (e.g. for soil moisture,  Klingmüller and Lelieveld, 2021).

The analysis of data-driven models will provide new insights into the nature of dust mobilisation processes and their dependency on environmental conditions. This might also allow the refinement of existing parameterisations of mineral dust emissions, e.g. the surface friction velocity or soil moisture terms.

Relating emission sources or surface concentrations to observed column burdens of gases and particles, being influenced by transport and other processes at the surface and in the atmosphere, is an emerging challenge in applications of satellite observations. The concept presented here in the context of aeolian dust can be generalised and applied to other atmospheric constituents for which satellite retrievals are available, for example additional aerosol components, nitrogen oxides or water vapour.

Appendix A

The submodels used in the EMAC simulations are listed in Table A1.

Lauer et al. (2007); Klingmüller et al. (2014)Pozzer et al. (2006)Kaiser et al. (2012)Jöckel et al. (2006)Dietmüller et al. (2016)Jöckel et al. (2006)Jöckel et al. (2006)Kerkweg et al. (2006a)Pringle et al. (2010)Sander et al. (2014)Tost et al. (2007)Sander et al. (2019)Kerkweg et al. (2006b)Kerkweg et al. (2006b)Dietmüller et al. (2016)Tsimpidi et al. (2014)Dietmüller et al. (2016)Tost et al. (2006)Kerkweg et al. (2006a)Kerkweg et al. (2006b)

Table A1 EMAC submodels used in the present study.