Variability and combination as an ensemble of mineral dust forecasts during the 2021 CADDIWA experiment using the WRF 3.7.1 and CHIMERE v2020r3 models

As operational support to define the Clouds–Atmospheric Dynamics–Dust Interactions in West Africa (CADDIWA) field campaign which took place in the Cape Verde area, the coupled regional model WRF–CHIMERE is deployed in forecast mode during the summer 2021. The simulation domain covers West Africa and the eastern Atlantic and allows the modeling of dust emissions and their transport to the Atlantic. On this route, we find Cape Verde, which was used as a base for measurements during the CADDIWA campaign. Meteorological variables and mineral dust concentrations are forecasted on a horizontal grid with a 30 km resolution and from the surface to 200 hPa. For a given day D, simulations are initialized from D−1 analyses and run for 4 d until D+4, yielding up to six available simulations on a given day. For each day, we thus have six different calculations, with better precision expected the closer we get to the analysis (lead D−1). In this study, a quantification of the forecast variability of wind, temperature, precipitation and mineral dust concentrations according to the modeled lead is presented. It is shown that the forecast quality does not decrease with time, and the high variability observed on some days for some variables (wind, temperature) does not explain the behavior of other dependent and downwind variables (mineral dust concentrations). A new method is also tested to create an ensemble without perturbing input data, but considering six forecast leads available for each date as members of an ensemble forecast. It has been shown that this new forecast based on this ensemble is able to give better results for two AErosol RObotic NETwork (AERONET) stations than the four available for aerosol optical depth observations. This could open the door to further testing with more complex operational systems.

1 Introduction

Over western Africa and during the boreal summer, mesoscale convective systems move from east to west and interact with the African easterly jet (AEJ), African easterly waves (AEWs) and mineral dust plumes (Knippertz and Todd, 2010; Marsham et al., 2011; Cuesta et al., 2020). Leaving the African continent to arrive above the Atlantic Ocean, they can generate tropical storms. The magnitude of interactions between these systems and the mineral dust concentrations via the direct and indirect effects of aerosols on meteorology remain unclear (Lavaysse et al., 2011; Price et al., 2018; Martinez and Chaboureau, 2018). This motivated the deployment of the Clouds–Atmospheric Dynamics–Dust Interactions in West Africa (CADDIWA) field campaign (Flamant et al., 2022). The measurements include long-term surface stations and dedicated airborne measurements. Aircraft were located at Sal (Cape Verde), under the wind flow coming from western Africa. In addition to the local study of storm generation, these aircraft measurements were also designed to help with the validation of spaceborne wind and aerosol products such as those of Aeolus, EarthCare and IASI satellites missions (Clerbaux et al., 2009; Illingworth et al., 2015; Martin et al., 2021). In addition to these measurements, numerical modeling is performed with the coupled regional model WRF–CHIMERE. Simulations are initialized at D−1 in order to provide forecasts from D+0 to D+4.

Figure 1 Model domain in red, with the main studied locations from east to west: Bodélé, Zinder, Banizoumbou, Niamey, Cinzana, Dakar and Cape Verde. The black circle is the possible range of aircraft measurements during the campaign.

Before studying the interaction between aerosols, clouds and radiation using a numerical tool, it is important to assess its accuracy. Forecast is a useful tool for this kind of evaluation. Simulating the same day several times with a different meteorology is a way to quantify the model variability. It is close to an ensemble simulation, even if the number of members is lower (Atger, 1999; Toth et al., 2001; Richardson, 2001). Comparing the 6 d of forecast simulations ranging from D−1 to D+4, but for a given date, enables quantifying the variability of the model.

Another aspect will be analyzed in this study: as several forecast leads correspond to the simulation of the same period but with different initial conditions for the meteorology, we can imagine that the leads are equivalent to ensemble modeling members. Ensemble modeling is widely used in forecasts of meteorology and air quality (Delle Monache et al., 2006; Vautard, 2006; Benedetti et al., 2018). But in general, the ensemble is built using the same model with perturbations. Some other techniques exist such as the “poor man's ensemble” and are widely used in meteorology (Ebert, 2001; Buizza et al., 2003; Bowler et al., 2008). To our knowledge, these approaches are not used for chemistry transport models (CTMs). They consist of using different models but making the forecast for the same period. Also in meteorology, they can be used in operational centers to update the covariance matrixes used for the data assimilation.

In this study, we aim to answer the following question: is it possible to use several forecast leads as an ensemble and improve the quality of the forecast? If the result is positive, it means it is possible to run fewer ensemble simulations each day (hence a faster forecast) while still improving the quality of the forecast. And for institutes which do not have sufficient computing resources to perform classical ensemble simulations, it still allows them to have a probabilistic approach to their forecast based on a single model.

The main goal of this study is thus first to quantify the variability in temperature, wind, precipitation rates, aerosol optical depth (AOD) and the surface concentration of mineral dust as a function of the forecast lead time. The second is to try to establish some correlations between possible differences in forecast results. With this quantification, we can assess the robustness of the forecast and the degree of confidence available that experimenters may have during field airborne campaigns such as CADDIWA. Section 2 presents the modeling system and the studied period. Section 3 presents the results of the comparison between the forecast with different lead times. Section 4 presents a tentative approach of mixing several leads for the same day in order to mimic an ensemble forecast. Results are compared to AErosol RObotic NETwork (AERONET) aerosol optical depth measurements. Section 5 presents the conclusions of the study.

2 The modeling system

2.1 The modeling tools

In this study, we use the WRF–CHIMERE model built with WRF 3.7.1 (Powers et al., 2017) and CHIMERE 2020r3 (Menut et al., 2021). These two models are coupled using the OASIS3-MCT external coupler (Craig et al., 2017). The WRF model is forced with the global-scale forecast fields from the NCEP Global Forecast System (GFS) (Halperin et al., 2020). For this experiment, a specific configuration was designed in order to have a lower numerical cost. Indeed, the goal was to launch the simulation for 6 d, from (D−1), the day before the current date, to (D+4) 4 d in advance. This long forecast was designed to allow the aircraft scientists to have enough time to decide what flight plan to use, depending on the meteorological situation to come. Daily simulations were launched at midnight to benefit from the latest forecast meteorological field and needed to be available to scientists by 08:00 local time in Cape Verde (05:00 UTC), including all post-processed figures. These constraints of real-time forecasts led to a light version of the model wherein only mineral dust is modeled. The model is also used in offline mode, meaning that there are no feedbacks of aerosols on meteorology, in order to ensure the stability of the calculation. Mineral dust is modeled with 10 bins from 0.01 to 40 µm. The dust emissions scheme used is the one of Alfaro and Gomes (2001), modified by Menut et al. (2005). Note that the CHIMERE model is also used daily in forecast mode for air quality with all available chemical processes, being operated by operational centers such as Prevair and Copernicus (Rouïl et al., 2009; Marécal et al., 2015).

The model domain in Fig. 1 is defined with the same horizontal grid for WRF and CHIMERE and covers part of the Atlantic Ocean and West Africa, from −40 to +20^∘ E in longitude and 0 to 33^∘ N in latitude. It is constituted of 200×110 cells with a constant resolution of 30 km. The WRF model has 32 vertical levels from the surface to 50 hPa. CHIMERE has fewer vertical levels with 15 layers from the surface to 300 hPa. This domain was designed to be able to simultaneously (i) model mineral dust emissions in Africa from Dakar to Bodélé, (ii) model transport from Africa to the Atlantic Ocean, and (iii) have the measurement site of Cape Verde not too close to the domain boundaries. In Fig. 1, the locations of Cape Verde, Dakar, Bodélé, Zinder, Niamey, Banizoumbou and Cinzana are reported. Model results were extracted daily at these locations for the CADDIWA scientists based in Cape Verde with the aircraft. The circle indicates the possible range of aircraft measurements during the campaign around the island of Sal.

2.2 The observations

The goal is not to perform a comparison of the model with observations. It will be done only at the end of the study with a comparison of measured versus modeled aerosol optical depth (AOD) and for the 2 m temperature and 10 m wind speed for some locations.

For the AOD, the AErosol RObotic NETwork global remote sensing network (AERONET, https://aeronet.gsfc.nasa.gov/, last access: 20 July 2023) level 1.5 measurements are used (Holben et al., 2001) (Table 1). The AOD values at a wavelength of λ=675 nm are averaged daily and compared to daily averaged modeled values. For the meteorological variables, the measurements provided by the Weather Information website of the University of Wyoming (UWYO) are used (http://www.weather.uwyo.edu/, last access: 20 July 2023). Data are provided for 2 m temperature and 10 m wind speed. It is noticeable that the data are delivered as integer values, restraining the accuracy of the comparison to the model results.

Table 1 List of the AERONET and meteorological UWYO sites used for the comparisons between measured and modeled AOD, 2 m temperature, and 10 m wind speed. Information includes the longitude λ and latitude ϕ for each site.

2.3 The modeled period and the intensive observation periods

The field campaign was carried out from 8 to 21 September 2021, with airborne measurements around Sal island in Cape Verde (Flamant et al., 2022). In order to have a tested and robust forecast modeling system, the daily forecast started on 10 August and ended on 1 November 2021. Among all observations periods, two events were observed: the tropical perturbation called Pierre-Henri, passing south of Sal on 11 September, and the period from 17 to 24 September with the passage over Sal of the two tropical cyclones called Peter and Rose. In this study, the results will be presented over two periods.

• Section 3 is only for the period 1 to 30 September 2021 for the variability of the forecast during the CADDIWA field experiment.

• Section 4 shows results over the whole modeled period from 10 August to 1 November 2021 for the merging of several forecast leads of modeled fields.

Figure 2 Principle of the modeling system in forecast mode. Each day, the global meteorological fields are downloaded to force the regional WRF model. These regional fields are used to drive the CHIMERE chemistry transport model, mainly for mineral dust emissions, transport and deposition. The procedure is repeated every day.

For the result presentation, there are several possibilities. As presented in Fig. 2, each day, the modeling system runs to simulate 6 d, from (D−1), the day before, to (D+4) 4 d in advance. With all these simulations, results may be discussed in two ways.

1. Comparison of all leads for one date. For example, for 11 September at 16:00 UTC, we can display the result of the simulations performed.

   • 12 September, forecast hour −8, (D−1)

   • 11 September, forecast hour 16, (D+0)

   • 10 September, forecast hour $\mathrm{16}+\mathrm{24}=\mathrm{40}$, (D+1)

   • 9 September, forecast hour $\mathrm{40}+\mathrm{24}=\mathrm{64}$, (D+2)

   • 8 September, forecast hour $\mathrm{64}+\mathrm{24}=\mathrm{88}$, (D+3)

   • 7 September, forecast hour $\mathrm{88}+\mathrm{24}=\mathrm{112}$, (D+4)

   This comparison may be achieved with maps and vertical cross-sections.

2. Comparison between leads during the whole period. It is possible to build time series using (D−1) for all days and (D+0) for all days until (D+4) for all days. In this case, we can calculate statistical scores between the time series as if they were different model realizations.

In the following sections and the Appendix, when the analysis consists of maps or vertical cross-sections, we selected 11 September 2021 to present the results, which is the day when the Pierre-Henri tropical perturbation was diagnosed above Sal in the Cape Verde islands (Flamant et al., 2022).

3 Variability of forecast leads during CADDIWA experiment

Results are presented as statistical scores (defined in Sect. 3.3). They are calculated for data over Bodélé and Cape Verde. The main goal being to compare the simulation leads and evaluate the variability from one lead to the next, there are no measurements in the analysis but only model versus model. For the initialization of the model performed using analyzed meteorological fields (Halperin et al., 2020), the simulation of (D−1) is considered to be the reference.

Figure 3 Time series of surface mineral dust concentrations (µg m⁻³) for each lead and of differences between leads.

3.1 Time series of surface mineral dust concentrations

Time series are presented for two sites, Bodélé and Cape Verde. They are located on a Sahelian isolatitude transect and often used to quantify the amount of mineral dust emitted in the Sahara and after long-range transport of the dust (Marticorena et al., 2010). Figure 3 presents time series for the surface mineral dust concentrations (µg m⁻³). The variability of the forecast for these concentrations should be the result of a mix between the variabilities calculated with the 10 m wind speed (an important parameter for the emissions) and the precipitation (a major sink). In Bodélé, the surface concentrations vary a lot between 0 and 5000 µg m⁻³. It seems huge, but it is classical when just over the main Saharan source. The mass is composed of a large mass distribution, and the majority of big particles are deposited before being transported, close to the source. The variability from one lead to another is important and illustrates the impact of the wind speed variability. The most important differences are between −2000 and +2000 µg m⁻³ in Bodélé. A major forecast underestimation is calculated on 11 September with −2000 µg m⁻³, meaning that an important peak of surface concentrations was modeled for (D−1) and (D+0) but not in advance. In Cape Verde, the surface concentrations are lower, which is logical after long-range transport and because this site is not a mineral dust emission hot spot. The concentrations remain high with peaks around 700 µg m⁻³. It is not the case of 11 September, but peaks are noted on 9, 12 and 16 September mainly. Also for Bodélé, large variability is calculated for 11 September, with forecast differences up to −200 µg m⁻³. As for Bodélé, the model underestimates the concentrations when the forecast is in advance. In addition to these results, time series of 2 m temperature, 10 m wind speed and precipitation are presented in Appendix A1.

Figure 4 Maps of surface mineral dust concentrations (µg m⁻³) and aerosol optical depth for 11 September at 12:00 UTC. Values are displayed (a, c) for the forecast lead (D−1) and (b, d) for the differences between the forecast leads (D−4)–(D−1).

3.2 Maps of mineral dust concentrations and AOD

Maps are presented for surface mineral dust concentrations (µg m⁻³) and aerosol optical depth for 11 September at 12:00 UTC in Fig. 4. Note that complementary maps for wind speed and precipitation are presented in Appendix A2. The simulation shows a large mineral dust plume flowing from Africa to the northern Atlantic Ocean. The site of Cape Verde is under this plume, and the trajectory over the ocean corresponds to the low wind speed values. The differences show the same kind of dipole as diagnosed for the precipitation (Fig. A4), showing that the shift between the forecast leads directly impacts the surface concentrations. With large positive values over land and negative values over the sea, it is noticeable that the more recent forecast (D−1) diagnoses a larger wind speed and then a faster transport: the dust plume is more over land for (D+4) but has already arrived over the sea in (D−1). It means that over Cape Verde, the last forecast diagnosed higher dust concentrations than the previous forecasts. The aerosol optical depth represents the behavior of the mineral dust concentration well, even if it diagnoses the radiative effect of all aerosols in the whole atmospheric column. The shape of the plume is slightly different, and a larger spatial spread of the differences between the forecast leads is seen. The differences remain important in absolute values since they can reach ±0.75 when the maximum AOD is 2. The variability in the forecast is important and shows for this day that the forecast of (D+4) underestimated AOD over Cape Verde compared to (D−1). In addition to these horizontal maps, the same variables are analyzed in Appendix A3 as vertical cross-sections.

3.3 Statistical scores

Usually, the variables O_t and M_t stand for the observed and modeled values, respectively, at time t. In the case of this study, as we want to quantify the variability of the forecast, the variable O_t is the model realization at (D−1) and the variable M_t is the model realization at leads (D+0) to (D+4). The mean value $\stackrel{\mathrm{‾}}{X_N}$ is calculated as

$$\begin{array}{}\text{(1)}& \stackrel{\mathrm{‾}}{X_N}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{t=\mathrm{1}}^N{X}_t,\end{array}$$

with N being the total number of hours of the simulation. To quantify the temporal variability, the Pearson product moment correlation coefficient R is calculated as

$$\begin{array}{}\text{(2)}& R={\displaystyle \frac{\frac{\mathrm{1}}{N}{\sum }_{t=\mathrm{1}}^N(M_t-\stackrel{\mathrm{‾}}{M_t})\times (O_t-\stackrel{\mathrm{‾}}{O_t})}{\sqrt{\frac{\mathrm{1}}{N}{\sum }_{t=\mathrm{1}}^N(M_t-\stackrel{\mathrm{‾}}{M_t}{)}^{\mathrm{2}}\times \frac{\mathrm{1}}{N}{\sum }_{t=\mathrm{1}}^N(O_t-\stackrel{\mathrm{‾}}{O_t}{)}^{\mathrm{2}}}}}.\end{array}$$

The spatial correlation, noted R_s, uses the same formula type except it is calculated from the temporal mean averaged values of observations and the model for each location where observations are available:

$$\begin{array}{}\text{(3)}& R_{\mathrm{s}}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^I(\stackrel{\mathrm{‾}}{M_i}-\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{M}})(\stackrel{\mathrm{‾}}{O_i}-\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{O}})}{\sqrt{{\sum }_{i=\mathrm{1}}^I(\stackrel{\mathrm{‾}}{M_i}-\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{M}}{)}^{\mathrm{2}}{\sum }_{i=\mathrm{1}}^I(\stackrel{\mathrm{‾}}{O_i}-\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{O}}{)}^{\mathrm{2}}}}},\end{array}$$

where I is the number of stations. The root mean square error (RMSE) is expressed as

$$\begin{array}{}\text{(4)}& \mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{T}}\sum _{t=\mathrm{1}}^T{\left(O_{t,i}-M_{t,i}\right)}^{\mathrm{2}}}.\end{array}$$

To quantify the mean differences between several leads, the bias is also quantified as

$$\begin{array}{}\text{(5)}& \mathrm{bias}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{t=\mathrm{1}}^N(M_t-O_t).\end{array}$$

Table 2 Statistical scores for the daily averaged 10 m wind speed (m s⁻¹) and 2 m temperature (^∘C) in Bodélé and Cape Verde.

Results are presented in Table 2 for the surface meteorological variables, 10 m wind speed and the 2 m temperature. Results are presented in Table 3 for the total precipitation and the surface mineral dust concentrations. For each variable, the presented values are the mean value (averaged over the whole month of September 2021), the bias and the correlation. The line (D−1) is always empty since the bias and the correlation compared to (D−1) give the values 0 and 1, respectively.

For the 10 m wind speed, it is noticeable that the mean value does not evolve a lot during the forecast. In Bodélé, the bias is always lower than −0.2 m s⁻¹, with the negative value meaning that the (D−1) analysis simulation was the one with the highest wind speed. In Cape Verde, the mean value is larger, but the bias is lower. With a maximum of −0.11 and negative or positive values, there is not a lot of variability for this site. An important point is the correlation of the leads compared to the (D−1) one: the values are very low for the two sites: between −0.27 and 0.27 for the maximum values. It means that, from one day to the next, the time series vary a lot in frequency. If the mean values are close, the maximum wind speeds are not at the same time.

For the 2 m temperature, we can also observe an important lack of variability for the mean values. The 2 m temperature is higher (≈32 ^∘C in average) in Bodélé than in Cape Verde (≈24 ^∘C in average), clearly showing the difference between desert and maritime-influenced air around an island. The bias is negative over Bodélé, showing that the forecast tends to underestimate the temperature compared to the lead (D−1). The bias is positive over Cape Verde, but the values are so low that it is negligible. Over Bodélé, the correlation remains high (between 0.75 and 0.9), showing that the forecast over the desert is very stable. It is not the case over Cape Verde, with large variability between 0.08 and 0.57. And the decrease in the correlation is not linear with the increasing lead: (D+4) has a correlation of 0.29 when (D+2) has a correlation of 0.08. Having a stable forecast over land does not mean that the forecast is stable over the Atlantic Ocean, with the meteorological systems being completely different.

Table 3 Statistical scores for the daily cumulated total precipitation (kg m⁻² h${}^{-\mathrm{1}}\times \mathrm{100}$) and surface mineral dust concentrations (µg m⁻³) in Bodélé and Cape Verde.

The same type of score is displayed in Table 3 for the daily cumulated total precipitation (kg m⁻² h${}^{-\mathrm{1}}\times \mathrm{100}$) and the surface mineral dust concentrations (µg m⁻³) in Bodélé and Cape Verde. For the precipitation, it is remarkable to see that the correlation is always close to zero. It means that, from day to day, the precipitation varies a lot for a specific location. It is logical since precipitation is a threshold process not continuous in space and time, contrary to the temperature or the wind. The bias is important in both Bodélé and Cape Verde: it corresponds to having a forecast with precipitation and the next lead without for the same place and time. For this process, the statistical scores show that the forecast is very variable from one day to another.

For the mineral dust concentrations, the correlation is also low for both sites. Over Bodélé, the values are between −0.22 for (D+3) and +0.17 for (D+0), and over Cape Verde, values are between −0.30 for (D+1) and +0.19 for (D+3). The bias is non-negligible and may reach 10 % of the mean values. As for the other parameters, there is not a regular decrease with an increasing lead: the system is chaotic, and the instability of the forecast for dust concentrations reflects the instability of the mean wind speed over source areas, then emissions, then transport, and then concentrations at remote locations. A result common to all parameters is that the best scores (for bias and correlation) are often obtained for the lead time close to the analysis (D−1).

4 Merging the forecast leads to make an ensemble

The previous results showed that small variations of meteorological variables may change mineral dust concentrations a lot after long-range transport. This quantification was made with only model results in order to quantify the model's variability. However, for some locations, it is possible to compare the AOD to AERONET measurements. During the studied period, four stations are present in the modeled domain and have available data: Zinder (Niger), Banizoumbou (Niger), Cinzana (Mali) and Cape Verde. Note that, unfortunately, there is no measurement for this period at Bodélé. Using these data, it is possible to calculate statistical scores between the modeled forecast and the measurements.

An added value in this study is that it is also possible to add two model realizations. Considering that the various forecast leads are performed each time with a new meteorology and then natural emissions (here mineral dust emissions), we can consider all these leads to be independent simulations. They are thus similar to ensemble forecast members, usually made the same day but with perturbed initial conditions. As presented in Fig. 2, for one date we have six simulations. It is possible to hypothesize that these six forecast leads are equivalent to six ensemble members. To test this approach, we use the time series at the four locations where AERONET measurements are available to create two new sets called ENSmean and ENSmedian.

• ENSmean corresponds to the mean averaged value of the six members.

• ENSmedian corresponds to the median of the members. Having only size members, this value is in fact the mean average of the third and fourth members.

Table 4 For the period of 1 to 31 September 2021, correlation (R), RMSE and bias calculated between the AERONET aerosol optical depth measurements and the modeled results. Results are presented for four sites, Cape Verde, Cinzana, Banizoumbou and Zinder, and for six forecast leads from (D−1) to (D+4). Two additional forecast leads called ENSmean and ENSmedian represent the mean average of the previous six leads and the median, respectively. The best scores for each site and among all leads are bolded.

4.1 Scores during the CADDIWA period

Statistical scores are first calculated for the period of the CADDIWA experiment from 1 to 31 September 2021. Results are presented in Table 4. For each site and each parameter (correlation R, RMSE and bias), the best score is bolded. It shows that model realizations are close to each other but remain different to the measurements. The variability of the forecast is lower than the difference between measurements and models. It means that the model systematically underestimates the AOD whatever the perturbations included in each forecast realization. The bias is negative for all sites over land (Zinder, Banizoumbou and Cinzana) and over sea around the Cape Verde islands. The bias is smaller over the latter site.

For the correlation, the best value is obtained for the ENSmean lead for two sites out of four: Cape Verde and Zinder. For Cinzana and Banizoumbou the best correlation is obtained for (D+3) and (D+2), respectively. It means that (i) the best scores are not for the “analysis” lead (D−1) as could be expected, and (ii) the combination of leads leading to ENS (ensemble) may be the best forecast. For the bias, results are different: the lower biases are not for ENS. In Cape Verde, the lower bias is for (D−1). For the other sites, they are for (D+3) and (D+4) forecasts. Note that even if ENS does not have the best score, it is also not the worst. These scores show that the best forecast lead is not always the closest to the analysis. It also shows that the use of an “ensemble” lead may provide good results.

For the RMSE, some of the best scores are also obtained for the ENSmean configuration, showing that merging the leads may reduce the model error. For the bias, the values remain very close from one lead to another and there is not really a best configuration.

Table 5 For the period of 1 to 31 September 2021, spatial and temporal correlation, RMSE, and bias for each lead and as an average for the four stations for the AOD measured by AERONET and modeled by WRF–CHIMERE. The best scores for each site and among all leads are bolded.

Table 5 summarizes the results presented in Table 4 by recalculating the scores but for the four stations (Cape Verde, Cinzana, Banizoumbou and Zinder) together. As for the previous results, the best scores are bolded as a function of the forecast lead. The spatial correlation R_s is the best for (D+3), and the correlation R is the best for the ENSmean forecast. The RMSE and bias remain the same for all model realizations.

Table 6 For the period of 1 to 31 September 2021, spatial and temporal correlation, RMSE, and bias for each lead and as an average for the four stations for the 2 m temperature measurements provided by UWYO and modeled by WRF–CHIMERE. The best scores for each site and among all leads are bolded.

The same type of score is presented for the 2 m temperature (Table 6) and the 10 m wind speed (Table 7). Scores are calculated using the UWYO meteorological data. They are hourly, but the problem is that they are recorded in integer form, decreasing their accuracy and possibly biassing the calculation of differences between observations and model results. It is interesting to explore the statistical scores of these two parameters since they are good proxies for mineral dust emissions: the 10 m wind speed is directly used for the saltation process via the friction velocity u_*, and the 2 m temperature is used to diagnosed the additional free convection velocity w_* (Menut et al., 2013).

For the 2 m temperature, the bias is positive as confirmed by the time series presented in Fig. 5. This bias varies a lot between leads, and the lower bias is for (D+4). The spatial correlation, R_s, is high but more or less constant between leads with values from 0.89 to 0.92. It means that the differences between sites remain close between the leads. The temporal correlation R ranges from 0.24 to 0.34 for (D−1). The ensemble leads provide correct scores with R=0.32 and 0.33.

For the 10 m wind speed, results are more variable. The spatial correlation ranges from 0.79 to 0.89 for (D+3). There is no regular decrease in the score with the lead. The ensemble is not the best score, but with R_s=0.84, the spatial correlation is better than (D−1) or (D+0). The temporal correlation is not correct and close to zero. As presented in Fig. 5, the modeled wind speed does not follow the day-to-day variations observed with the measurements. But comparing observed and modeled wind speed remains challenging: first for the integers recorded with the stations and second with the differences of representativity with a specific observation site; on the other hand, there is a model cell of a few tens of square kilometers. Finally, the scores for 2 m temperature and 10 m wind speed are not able to completely explain the scores obtained with the ensemble lead.

Table 7 For the period of 1 to 31 September 2021, spatial and temporal correlation, RMSE, and bias for each lead and as an average for the four stations for the 10 m wind speed measurements provided by UWYO and modeled by WRF–CHIMERE. The best scores for each site and among all leads are bolded.

Figure 5 Time series of 2 m temperature (^∘C) and 10 m wind speed (m s⁻¹) measured (UWYO database) and modeled during the forecast and for several leads. The last time series, called ENSmean and ENSmedian, are the mean averaged and the median values of the previous leads from (D−1) to (D+4).

4.2 Scores during the extended forecast period

In order to have more statistically robust results, the complete modeled period is now presented: 15 August to 1 November 2021. This period is around the CADDIWA experiment and corresponds to the period when the forecast system was running, i.e., 2.5 months. Results are presented as time series in Fig. 6 for the daily averaged AOD.

Figure 6 Time series of aerosol optical depth (AOD) measured by AERONET and modeled during the forecast and for several leads. The last time series, called ENSmean and ENSmedian, are the mean averaged and the median values of the previous leads from (D−1) to (D+4).

It is noticeable that the month of September (compared to August and October) is not the month with the highest AOD: values are of the same order of magnitude over the whole period, except in Cape Verde where the largest peaks are observed during September 2021, corresponding to the CADDIWA measurement campaign.

Table 8 For the period 15 August to 1 November 2021, correlation (R), RMSE and bias calculated between the AERONET aerosol optical depth measurements and the modeled results. Results are presented for four sites, Cape Verde, Cinzana, Banizoumbou and Zinder, and for six forecast leads from (D−1) to (D+4). Two additional forecast leads called ENSmean and ENSmedian represent the mean average of the previous six leads and the median, respectively. The best scores for each site and among all leads are bolded.

Statistical scores are presented in Table 8 in the same way as in Table 4 but this time for a longer period. Over this period, the availability of hourly measurements is 62.5 %, 77.8 %, 79.2 % and 77.8 % for Cape Verde, Cinzana, Banizoumbou and Zinder stations, respectively. For this longer period, the best correlations are not for the ensemble leads, except for the RMSE in Cape Verde and Zinder and for the bias in Cape Verde. For the correlation, the best scores are now for the first forecast leads, i.e., (D−1) and (D+0). All in all, the scores are very close from one lead to the next one.

Table 9 Comparison between observations and the model for the AOD. For the period 15 August to 1 November 2021, spatial and temporal correlation, RMSE, and bias for each lead and as an average for the four stations. The best scores for each site and among all leads are bolded.

Table 9 summarizes the results presented in Table 8 as in Table 5. The best spatial correlation is again for (D+3) when the best temporal correlation is obtained for the leads (D−1) and (D+0). For the RMSE, scores are very close between leads, and the best values are for (D−1), (D+0) and (D+1), but also for the ensemble, with both ENSmean and ENSmedian.

5 Conclusions

In this study, the first goal was to examine the variability of the forecast as a function of the lead time and for each forecasted day. This forecast was performed daily for 6 d during the period August to October 2021 and as support for the CADDIWA field campaign. For meteorological variables (2 m temperature, 10 m wind speed, total precipitation rate) and surface concentrations of mineral dust, the day-to-day variability was quantified. The performances of the forecast over two sites were evaluated: Bodélé (desert area and important source of dust) and Cape Verde (where the measurements of DACCIWA were coordinated). It has been shown that the wind speed is highly variable for day-to-day forecast, while the temperature is stable over land but more variable over sea and shores (Cape Verde being a group of little islands). The less stable parameter is the precipitation at one location when the model may forecast an event one day and not at all the day after.

First, one goal of the study was to examine whether large forecast variability at one site (such as Bodélé) may have a visible impact at a downwind remote site (such as Cape Verde). No evidence of a transport of variability (or a transport of stability) was found during the forecast. The large variability of wind speed, precipitation and temperature induces large variability of the surface concentration of mineral dust. Between forecast leads, large differences were found both for the correlation and the bias. Considering the model configuration used for this study, wherein no direct or indirect effects of aerosols on meteorology and only mineral dust as natural emissions were taken into account, this variability could be underestimated. A next study could be to replay this forecast with a model version including all anthropogenic and natural emissions in the CHIMERE model with an exhaustive evaluation with the measurements of the experiment to come.

Second, a new way of combining forecast leads was tested to improve the predictions. Considering that several forecast leads may be considered to be members of an ensemble, they are combined from (D−1) to (D+4) for all coinciding dates by computing the mean and median values. These new “forecast leads” are compared, with all others members, to the aerosol optical depth measurements of AERONET using correlation, RMSE and bias statistics. It is noticeable that the forecast is not impaired when increasing the lead time. But it is also noticeable that out of four sites, the best scores for two sites are with the ensemble for the period of the CADDIWA campaign. It is not the case for an extended analyzed period, highlighting that the scores are close from one lead to another. The ensemble methodology provides the best scores when the AOD values are the most important and the most variable in time. This result opens perspectives for forecasting in general. It would be interesting to test this hypothesis on operational systems: if the combination of the previous forecasts allows improving the initial conditions of a new forecast, it would allow performing fewer ensemble simulations for the same day and thus considerably reducing the computing cost.

Appendix A: Complementary analysis with the meteorological variables

A1 Time series of meteorological variables

In addition to the time series presented in Sect. 3.1, the same results are presented here for 2 m temperature, 10 m wind speed and precipitation rate.

Figure A1 presents time series and differences for 2 m temperature (^∘C) in Bodélé and Cape Verde. The days of 11 and 18 September are noted in the figure with a black vertical line. In Bodélé the temperature is higher than in Cape Verde, with values between 30 and 35 ^∘C. During the month of September the temperature decreases regularly. The days of 11 and 18 September correspond to periods with the highest temperature values. In Cape Verde, there is no similar trend: the daily averaged temperature remains around 25 ^∘C, showing the maritime characteristic of the Cape Verde environment. The differences are low and oscillate between −2 and +2^∘C. The longer the forecast, the greater the variability. In Cape Verde, the variability is lower and between −1 and +1^∘C. As in Bodélé, the largest differences with (D−1) are obtained with (D+4). The forecast of temperature appears to be relatively stable, with the differences logically growing with the increasing leads.

Figure A1 Time series of 2 m temperature (^∘C) for each lead and of differences between leads.

Figure A2 presents results for 10 m wind speed. Values are lower in Bodélé (middle of the desert) than in Cape Verde (a group of islands). In Bodélé, daily averaged values are between 2 and 7 m s⁻¹, which are values lower than the minimum value generally required for dust erosion over barren soils. But hourly values may be larger and the model uses a Weibull distribution to take into account the sub-hour and the sub-grid spatial variability (Menut, 2018). It is noticeable the two days of 11 and 18 September do not correspond to a high wind speed value, and the days before also do not. In Cape Verde, the wind speed values are between 3 and 10 m s⁻¹, with day-to-day variability higher than in Bodélé. Some days shows high values such as 5 and 11 September. It is the signature, close to the surface, of large-scale meteorological motions.

In Fig. A2, differences are also presented. Differences are of the same order of magnitude between the two locations. For (D+0)–(D−1), differences are maximum ± 0.5 m s⁻¹ when higher values are calculated for (D+4)–(D−1) with a maximum around ± 3 m s⁻¹. It is noticeable that the differences increase with the lead: the more distant the forecast, the greater the difference between the leads. For these differences, there is no systematic bias: they can be negative or positive, showing variability not due to large-scale and/or persistent atmospheric systems, but much more regional variability, with a higher temporal frequency. More specifically for 11 September, when the absolute values shows a peak in Cape Verde, the differences show that this peak was predicted late: 4 d before, for the (D+4) forecast, the peak is 6 m s⁻¹, when it is 9 m s⁻¹ for (D+0). The difference is then −3 m s⁻¹, which is one of the most important during the whole modeled period.

Figure A2 Time series of 10 m wind speed (m s⁻¹) for each lead and of differences between leads.

Figure A3 presents the same kind of time series but for total precipitation in kg m⁻² h${}^{-\mathrm{1}}\times \mathrm{100}$. The time series show that only a few periods had precipitation episodes for the two sites of Bodélé and Cape Verde. In Bodélé, the two periods with rain are 6 and 15 September. In Cape Verde, three episodes are modeled: 6 and 14 September and 5 October (the last one outside the current analyzed period). For the first episode in Bodélé, 6 September, it appears only for the forecast lead (D+4). For the other forecasts, closer in time, there is no precipitation. For the second episode, time variability is observed: depending on the lead, the precipitation episodes have similar magnitude but are forecasted on 14, 15 or 16 September. The difference shows that the forecast is mostly overestimated compared to the analysis of (D−1). In Cape Verde, several precipitation episodes also vary in time. If the first episode is overestimated for the (D+4) lead, it is finally underestimated by the other leads, from (D+1) to (D+3). The second episode is forecasted with less variability in time, with all forecasts being for 13 or 14 September only. The magnitudes are close between leads, with only a low underestimation compared to (D−1). Finally, the forecast is less variable in Cape Verde than in Bodélé.

Figure A3 Time series of total precipitation (kg m⁻² h⁻¹ ×100) for each lead and of differences between leads.

Figure A4 Maps of 10 m wind speed (m s⁻¹) and total precipitation (kg m⁻² h⁻¹) for 11 September at 12:00 UTC. Values are displayed (left) for the forecast lead (D−1) and (right) for the differences between the forecast leads (D−4)–(D−1).

Figure A5 Vertical cross-section of mineral dust concentrations (µg m⁻³) and precipitation rate (in kg kg${}^{-\mathrm{1}}\times {\mathrm{10}}^{\mathrm{6}}$) at isolatitude 17^∘ N. Figures represent the same date, 11 September 2021 at 12:00 UTC. (a, c) Absolute values for forecast lead (D−1) and (b, d) differences between forecast leads (D+4)–(D−1). The line in red is the boundary layer height.

A2 Maps of wind speed and precipitation

Results are presented as maps for one date, 11 September at 12:00 UTC. Figure A4 first presents maps for the 10 m wind speed (m s⁻¹) and total precipitation (kg m⁻² h⁻¹). In the left panel, the absolute value of the forecast lead (D−1) is presented. In the right panel, the differences between the leads (D−4) and (D−1) are shown. Note that for the wind speed, the wind vectors are superimposed. For the 10 m wind speed, the values of (D−1) show moderate values (between 0 and 3 m s⁻¹), except over Mauritania with maximum values ≈15 m s⁻¹. The wind speed is larger over the Atlantic Ocean with values ≈8 m s⁻¹ near Cape Verde. The Cape Verde site is between two different air masses: one coming from the south and evolving along the African coast and the second one on the west side of Cape Verde coming from the north. It results in low wind locally in Cape Verde. The map of differences shows the same pattern, meaning that this structure changes during the forecast: the (D+4) forecast shows negative values, meaning that the wind speed is higher for (D−1) than (D+4). It means that the strong gradient, from northeast to southwest and flying over Cape Verde, observed for (D−1) was not present for the forecast 4 d in advance.

For the total precipitation in Fig. A4, the results on the map show very localized events. Over Africa and over Mauritania, the large amount of precipitation (≈1 kg m⁻² h⁻¹) is colocated with the large 10 m wind speed values. Other precipitation events are modeled more in the south over both the Atlantic Ocean and the Gulf of Guinea for a latitude below 10^∘ N. The difference map shows negative and positive values: it is the mark of a change in wind speed and direction, then location of the precipitation. But whatever the location, each time precipitation was forecasted, it occurred each time, even if not strictly at the same place. For this day, no precipitation was forecasted over Cape Verde, and this forecast remained stable during several forecast leads.

A3 Vertical cross-sections of mineral dust concentrations and rain

Figure A5 presents a vertical cross-section of mineral dust concentrations (µg m⁻³) and precipitation rate (in kg kg${}^{-\mathrm{1}}\times {\mathrm{10}}^{\mathrm{6}}$) at isolatitude 17^∘ N for the forecast lead of (D−1) and the difference between (D+4) and (D−1). The same day and hour as for the horizontal maps are selected for these results.

The goal of these figures is to present the vertical extent of the possible differences between the leads. For the mineral dust concentrations, the large surface concentrations extend vertically until 3000 m. And concentrations are non-negligible until 7000 m. At the longitude of Cape Verde, −23^∘ W, dust concentrations are large, but the forecast is very variable. The differences show maximum values between −20 and −10^∘ W. Around −20^∘ W, the vertical structure shows negative values close to the surface but positive values between 1500 and 3000 m, above the boundary layer. It means that the wind direction changed between the forecast leads but also the vertical distribution of the dust plume coming from Africa. It explains the differences for the surface concentrations and should also have an impact on AOD (see Sect.  4).

The vertical profile of rain shows a large event for this day at longitude $\approx -\mathrm{13}$^∘ W. It corresponds to the event seen in Fig. A4 over the southwest of Mauritania. The vertical cross-section of differences shows negative then positive values: because the wind is faster as the forecast is close from the current day, the precipitation is transported faster and then appears as positive for longitude −13^∘ W and negative in longitude −16^∘ W. If the horizontal transport changes with leads, the vertical structure remains the same with a maximum at 6000 m.