Position correction in dust storm forecasting using LOTOS-EUROS v2.1: grid-distorted data assimilation v1.0

When calibrating simulations of dust clouds, both the intensity and the position are important. Intensity errors arise mainly from uncertain emission and sedimentation strengths, while position errors are attributed either to imperfect emission timing or to uncertainties in the transport. Though many studies have been conducted on the calibration or correction of dust simulations, most of these focus on intensity solely and leave the position errors mainly unchanged. In this paper, a grid-distorted data assimilation, which consists of an image-morphing method and an ensemble-based variational assimilation, is designed for realigning a simulated dust plume to correct the position error. This newly developed grid-distorted data assimilation has been applied to a dust storm event in May 2017 over East Asia. Results have been compared for three configurations: a traditional assimilation configuration that focuses solely on intensity correction, a grid-distorted data assimilation that focuses on position correction only and the hybrid assimilation that combines these two. For the evaluated case, the position misfit in the simulations is shown to be dominant in the results. The traditional emission inversion only slightly improves the dust simulation, while the grid-distorted data assimilation effectively improves the dust simulation and forecasting. The hybrid assimilation that corrects both position and intensity of the dust load provides the best initial condition for forecasting of dust concentrations.

1 Introduction

Dust storms are a result of wind erosion liberating particles from exposed dry surfaces (World Meteorology Organization, 2019). They occur commonly in arid or semi-arid regions, e.g. North Africa, the Middle East, Southwest Asia and East Asia (Shao et al., 2013). During dust events, fine dust particles can be lifted several kilometres high into the atmosphere and carried over thousands of kilometres (Zhang et al., 2018). It is estimated that 2000 Mt dust is emitted into the atmosphere annually (Shao et al., 2011). Such a huge amount of atmospheric mineral dust has profound effects on the Earth system, e.g. the cycles of energy, carbon and water (Mahowald et al., 2010). Specifically, dust particles are recognized in fertilizing terrestrial and ocean ecosystem (Shepherd et al., 2016), enhancing precipitation by acting as droplet nuclei (Benedetti et al., 2014) and interacting with atmospheric radiation, and may therefore significantly modify the Earth radiative balance (Balkanski et al., 2007; Wu et al., 2016). Apart from the influence on the environment, dust storms pose a great threat to human health by carrying thousands of tonnes of particulate matter as well as bacteria, viruses and persistent organic pollutants to densely populated regions (World Meteorological Organization, 2017; Basart et al., 2019). Reported illnesses include dust pneumonia, strep throat, cardiovascular disorders and eye sicknesses (Shao and Dong, 2006; Ozer et al., 2007; Benedetti et al., 2014; World Meteorological Organization, 2018). The low visibility caused by dusts can also lead to severe disruptions of air and other traffic. For example, more than 1100 flights were delayed or cancelled in Beijing after it was struck by an extreme dust event in May 2017.

Together with growing interest in dust storms, the understanding of the physical processes associated with dust storms has increased rapidly over the last decades (World Meteorological Organization, 2018). Large efforts have been made to develop dust modelling systems (Marticorena and Bergametti, 1995; Shao et al., 1996; Marticorena et al., 1997; Alfaro et al., 1997; Wang et al., 2000; Liu et al., 2003; Basart et al., 2012), which mathematically simulate the life cycle of dust including emission, transport and deposition. Large-scale global dust transport models, e.g. CAMS-ECMWF (Morcrette et al., 2009), or regional ones, e.g. NASA-GEOS-5 (Colarco et al., 2010) and BSC-DREAM8b (Mona et al., 2014), are essential parts of larger Earth system models. The most important application of these models is to forecast dust concentrations over a few hours to a few days in order to reduce the potential threats to society. Though these systems are usually able to predict the starting and ending of a dust event, large differences are found in emission and deposition burdens and spatial distribution of dust clouds (Huneeus et al., 2011; Koffi et al., 2012). Dust simulations could differ from observations by up to 2 orders of magnitude (Uno et al., 2006; Gong and Zhang, 2008). The modelling skills are limited due to several aspects, e.g. the insufficient knowledge of aerosol size distribution (Mokhtari et al., 2012), mismatch in aerosol removal (Croft et al., 2012) and in particular to the inaccurate quantification of erosive dust emission (Gong and Zhang, 2008; Ginoux et al., 2012; Escribano et al., 2016; Di Tomaso et al., 2017). In addition, the quality of the meteorological data, e.g. wind fields and soil moisture, might strongly impact the prognostic quality of dust emission and transport.

In addition to the efforts of upgrading the physical descriptions in numerical models, data assimilation techniques have been developed to improve simulation of dust loads. Data assimilation aims here to estimate the state of dust concentrations by combining a dynamical model with available observations. An assimilation system could for example adjust model parameters within an allowed range such that a simulation is in better agreement with the observations. Various types of observations have been used to adjust dust simulations, for example particular matter (PM) measurements (Lin et al., 2008a; Wang et al., 2008) and visibility records (Niu et al., 2008; Gong and Zhang, 2008) from ground-based monitoring networks, aerosol optical depth (AOD) from sun photometers in the global Aerosol Robotic Network (AERONET), (Schutgens et al., 2012) and the satellite-retrieved AOD (Khade et al., 2013; Yumimoto et al., 2016; Di Tomaso et al., 2017; Dai et al., 2019). Those studies focused either on updating atmospheric dust concentrations directly or on optimizing emission parameters that lead to better simulations. In both cases, only the intensity of either concentrations or emissions is adjusted, while other input parameters are assumed to be known, and processes of transport and removal are assumed to be certain.

In our previous studies, ground-based PM₁₀ (total particulate matter with diameter less then 10 µm) measurements (Jin et al., 2018, 2019a) and geostationary satellite AOD (Jin et al., 2019b, 2020) were assimilated with the LOTOS-EUROS simulation model for dust storm forecasts over East Asia. Also these studies solely focused on correcting emission intensities. Data selection (Jin et al., 2019b) and observation bias correction (Jin et al., 2019a) were important aspects here to ensure that the available measurements were used correctly. In addition, an adjoint method was used to identify potential new dust emission sources in case the empirical dust emission and its uncertainty scheme cannot fully resolve the observation (Jin et al., 2020). Severe dust storm events in May 2017 over East Asia were used as test cases, and the assimilation procedure was shown to improve the simulated dust concentrations at the time of observation but also to improve forecasts of dust levels over windows of up to 24 h. During these studies it was noted that although the modelling system in general provided an accurate forecast of the dust plume, a severe position error was present when the plume travelled over a large distance. Specifically, forecasts by the model simulation reported the dust arrival and departure 1 to 10 h prior to reality, as is also illustrated in Sect. 3.

Position errors are a common problem in meteorology, for example in forecasting hurricanes, thunderstorms, precipitation (Ravela et al., 2007; Nehrkorn et al., 2014, 2015) or meteorology-governing events like wildfires (Beezley and Mandel, 2008). In geophysical disciplines, a positional error is often considered together with intensity errors to explain differences between two estimates (Nehrkorn et al., 2015). A misfit in position usually leads to significant degradation of forecasts (Jones and Macpherson, 1997).

When discussing the accuracy of a dust forecast, the shape and position of the plume is a key element as well as the intensity. The position forecast determines which locations will be affected, when the storm will arrive and for how long it will last, while the intensity only describes the actual dust level. A dust forecast with position misfit directly results in incorrect timing profiles of dust loads. The information about dust arrival and departure is sometimes more important than the magnitude of dust load in the early warning system, but until now it has attracted only little attention.

Facing the unresolved positional mismatch, the aforementioned data assimilation focusing solely on intensity correction is less effective, as is illustrated in Sect. 4.1.

Similarly to intensity feature misfits, positional misfits in model simulations can also be adjusted to better resemble observations using data assimilation techniques. Dust simulations suffer from position errors due to for example incorrect emission timing profiles or uncertainties in the transport, both driven by uncertain meteorology fields. To be able to use data assimilation techniques for position correction, it is essential to have a description of these uncertainties. However, position errors are much likely to be non-Gaussian and not easily captured by a static error covariance model (Nehrkorn et al., 2015). For dust simulation, position errors could be caused by uncertainties in the transport, in particular the wind field. These uncertainties accumulate during the time period from emission in remote desert areas to arrival at observation networks in downwind populated areas. Position discrepancies might also arise from incorrect timing profiles of emissions, which is not the case for our test event, as is explained in Sect. 3.1. However, determining the covariance either for transport or for emission timing profile is difficult. Even if there is a complex covariance model that could account for the accumulation of uncertainties along the long track of the plume, a substantial number of observations would then be required to constrain the optimal transport pattern. Data assimilation methods based on static covariance models are therefore often not suitable for dealing with position errors.

Instead, techniques from the field of image processing could be combined with data assimilation to avoid the need for a static covariance that describes the origin of the position error. This has been described as phase-correcting data assimilation in numerical weather prediction (Brewster, 2003), image-morphing ensemble Kalman filter (EnKF) for wildfire models (Beezley and Mandel, 2008), grid distortion data assimilation in oil reservoir modelling (Lawniczak, 2012), and in general as position error correction in variational data assimilation (Nehrkorn et al., 2015). The common approach in all these applications is to reposition the simulation using an image-morphing technique, where the optimal morphing parameters are adjusted to obtain the best fit with the observations using data assimilation techniques. In an application with dust plume simulations, the use of image morphing in the data assimilation avoids the need for developing a complex covariance model to describe uncertain transport or emission timing.

In this study, we propose a grid-distorted data assimilation method to correct position misfits in a simulated dust plume, which is a novel approach in the context of atmospheric dust modelling. The implemented method offers an efficient way to correct for a phase misfit between a dust simulation and available observations without changing the transport scheme and/or the emission timing profile. The grid-distorted data assimilation is then combined with the emission intensity inversion described in Jin et al. (2019b) for a hybrid method. The hybrid method is capable of optimizing the dust plume in case both position and intensity misfits are presented in a dust simulation. Starting from the initial condition using the hybrid assimilation posterior, dust forecasting accuracy (in terms of both arrival and departure and in actual dust load) is further ensured.

The paper is organized as follows. Section 2 introduces the simulation model and observations used to represent the dust intensity. Section 3 shows an example of a dust position error in a dust simulation. The error source is explained and identified to be the uncertainty in long-distance transport process, and it is illustrated that this uncertainty cannot be explained from the known spread in meteorological forecasts. In Sect. 4, the necessity of position error correction is emphasized first, and then the methodology of grid-distorted data assimilation is introduced. A hybrid assimilation method is designed by combining the grid-distorted data assimilation and emission inversion in Sect. 5. The new method is evaluated against assimilation focusing solely on emission intensities or position correction. Section 6 summarizes the conclusion and the added value of using grid-distorted data assimilation to resolve model position error.

2 Dust model and observations

2.1 Simulation model

In this study, the dust storm is simulated using a regional chemical transport model, LOTOS-EUROS v2.1 (Manders et al., 2017). LOTOS-EUROS has been used for a wide range of applications supporting scientific research and operational air quality forecasts both inside and outside Europe. At present, the operational forecasts over China are released via the MarcoPolo–Panda projects (Timmermans et al., 2017; Brasseur et al., 2019) through http://www.marcopolo-panda.eu/forecast/ (last access: July 2020). Additionally, it is also implemented in the World Meteorological Organization (WMO) Sand and Dust Storm Warning Advisory and Assessment System to provide short-time forecasting of the dust load over the North Africa–Middle East–Europe areas; the online forecast product is delivered via http://sds-was.aemet.es/forecast-products/dust-forecasts/compared-dust-forecasts (last access: July 2020).

To establish a dust simulation over East Asia, the model is configured on a domain from 15 to 50^∘ N and 70 to 140^∘ E, with a resolution of about $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$. Vertically, the model consists of eight layers, with a top at 10 km. The dust simulation is driven by European Centre for Medium-Ranged Weather Forecasts (ECMWF) operational forecasts over 3–12 h, retrieved at a regular longitude–latitude grid resolution of about 7 km. An interface to the ECMWF output set is designed, which not only interpolates the default 3 h ECMWF short-term forecast meteorology to hour values but also averages the forecast to fit the LOTOS-EUROS spatial resolutions (Manders et al., 2017). Physical processes included are wind-blown dust emission, diffusion, advection, dry and wet deposition, and sedimentation.

Figure 1 Distribution of the barren region over East Asia and the China Ministry of Environmental Protection (MEP) observing network.

Figure 2 Original PM₁₀ (a.1), bias-corrected dust observations (a.2), the a priori surface dust concentration (b), maximum over the ensemble simulations driven by ensemble meteorology (c), posterior dust simulation of the emis inversion (d), grid-distorted assim posterior (e) and hybrid assim posterior simulation (f) at 15:00 CST on 5 May. SDC: surface dust concentration. Definitions of emis inversion, grid-distorted assim and hybrid assim can be found in Table 1.

2.2 Observation network

The observations used in this study consist of hourly PM₁₀ concentrations from the China Ministry of Environmental Protection (MEP) air quality monitoring network, which is shown in Fig. 1. By now, the network has over 1700 stations and hence offers an opportunity to track the whole dust plume while it moves through the region.

All these PM₁₀ measurements are actually a sum of dust and airborne particles (black carbon, sulfate, etc). Since the analysed event is an extremely severe case, these PM₁₀ measurements were directly used to quantify the dust load in Jin et al. (2019b, 2020). In this study however, an observational bias correction is performed to make the PM₁₀ measurements fully representative of the dust loads. First, non-dust aerosol levels are calculated using a LOTOS-EUROS simulation following the MarcoPolo–Panda configuration but with the dust tracers disabled. Using these simulations, bias-corrected dust observations were calculated by subtracting the non-dust loads from the original PM₁₀ observations. The original PM₁₀ measurements vs. the pure dust observations can be seen in Figs. 2a.1 and a.2 and 7a.1 and a.2. As dust aerosols are far more dominant during the severe dust storm, the bias-corrected dust observations are actually very close to the original PM₁₀ measurements.

3 Position error

Numerical dust models are expected to provide correct timing profiles and intensity of dust loads. However, a discrepancy between observations and simulations is relatively common in terms of both position and intensity. Unlike the intensity estimation that has been widely investigated already, the position error has received less attention, but it has been the main focus of this study.

3.1 Position error in dust simulation

The test case investigated in this study is a severe dust storm event that occurred over East Asia in May 2017. The detailed calibration of the model simulations on this test case can be found in Jin et al. (2019b, 2020). The dust emission occurred from 2 May in the Mongolia, Gobi and Alxa deserts, of which the location can be seen in Fig. 1. The dust particles lifted up from these regions were then transported in the south-east direction. After 2 to 3 d of transport, the dust plume arrived in central China, where according to the surface observations a positional error was present in the simulations.

The position error in the simulation is illustrated in Fig. 2, which shows the original PM₁₀ measurements, bias-corrected dust observations and the a priori surface dust concentration (SDC) simulation on 5 May at 15:00 (China Standard Time, CST). The measurements of PM₁₀ are strongly elevated when the dust plume passes and could increase to values over 2000 µg m⁻³. Under normal conditions the observations (non-dust aerosols) usually do not exceed values of 200 µg m⁻³, and therefore the location of a dust plume is clearly visible in the bias-corrected dust observations as well as in these original PM₁₀ observations. According to the observations in panel (a), the dust plume forms a band from the west to the east over central China. The corresponding simulation in panel (b) shows a plume with a similar shape but at a location farther to the south-east. This is indicated by the markers that are added to the plumes. For the observations the markers for the left part of the plume are around 35^∘ N, and the right one stays around 37.5^∘ N, while for the simulation they are around 32.5 and 36^∘ N. The dust plume is therefore positioned about 200 km too far to the south; with a wind speed of 40 km h⁻¹ this implies a difference in arrival time of 5 h. The simulated plume, in particular the left part, is also broader than the rather sharp band that is seen in the observations.

To quantify the simulation-minus-observation mismatch, the root mean square error (RMSE) between dust simulation and bias-corrected dust observation has been computed over all stations in central China (marked by the black framework in Fig. 2a). The RMSE of the a priori dust simulation is as high as 388.1 µg m⁻³. This vast mismatch is attributed to the sum of intensity and position error (mainly), as is explained in Sect. 5.2.

Figure 3 PM₁₀ observations (a) and the a priori dust simulation (b) at 15:00 CST on 4 May. SDC: surface dust concentration.

Figure 4 Illustration of intensity-centred assimilation only (a) versus assimilation after position error correction (b).

3.2 Uncertainty in emission timing profile

One potential origin of the position error is an incorrect emission time profile. That is, changes in the time period over which dust is released from the source regions could to some extent alter the position of the simulated plume.

Actually during the first 48 h after dust emission started, the simulated dust plume was still in northern China and showed in general the same pattern as is visible in the observations. For example the aerosol optical depth (AOD) retrieved from the Himawari-8 geostationary satellite showed that the simulated plumes are correctly positioned in northern China (Jin et al., 2019b). The good phase match in general can also be seen from a snapshot of the ground PM₁₀ observation vs. the simulated surface dust concentration on 4 May at 15:00 CST in Fig. 3. There might already be position misfits in the dust simulation at these snapshots that are not easily detected. The magnitudes of the dust concentration showed discrepancies, but these could be corrected by emission inversion through assimilation of those AOD observations or PM₁₀ measurements. The good match in position between simulated and observed dust plume indicates that the emission timing profile is rather accurate too. When the dust plume is transported farther southward, the simulated plume starts to deviate from the available surface measurements.

3.3 Uncertainty in meteorology

Another possible origin of the position error in the simulations is the uncertainty in the meteorological data. In our study, the simulation model is driven by ECMWF meteorological forecasts. The uncertainty in this input is reflected in the ensemble forecasts that are available too (Palmer, 2019). For the studied period, the ensemble forecast of N_meteo=26 different members is available, where each member is a perturbation of the deterministic forecast. The resolution of meteorological ensemble is about 30 km, which is comparable to the LOTOS-EUROS resolution for these experiments.

To estimate the impact of the meteorological uncertainty, the dust simulations have been repeated N_meteo times using input from the meteorological ensemble. The spread in simulated dust concentrations is computed in terms of the maximum over the ensemble via

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{c}_{\text{max}}(x,y,z,t)=\text{max}\left(c_{\mathrm{1}}(x,y,z,t),\mathrm{\dots },c_{N\text{meteo}}(x,y,z,t)\right).\end{array}$$

In here, c_i represents the dust concentration field that results from a simulation with the ith ensemble member. This measure reflects for each location whether in any of the simulations a severe dust load is present. The ensemble maximum here can be used as a quick criterion: only if the dust plume (in observational view) is covered by the maximum is meteorological uncertainty (represented by ensemble meteorological inputs) likely to resolve the dust plume position error.

A snapshot of the ensemble maximum Eq. (1) at 15:00 CST is shown in Fig. 2c. The map shows a broader plume, which implies that some ensemble members result in a dust plume that is more to the north and others more to the south than the a priori forecast. The extended dust field is however not wide enough to cover the area with increased observation values. The uncertainty approximated using the available meteorological ensemble therefore could not be used to fully account for the position error. The origin might be that the required case is not represented in the ensemble but also because the simulated dust transport in the LOTOS-EUROS model does not take all meteorological details into account or is simply not accurate enough.

To resolve the position error, a complex covariance matrix would then be required to fully account for the accumulation of uncertainties along the long track of the plume. The uncertainty in the interface that interpolates and averages the meteorological forecast to fit our LOTOS-EUROS model resolution should also be taken into account here.

4 Grid-distorted data assimilation

The experiments in the previous section showed that the mismatch between dust plume simulation and observations cannot be easily explained by inaccurate emission timing or uncertainty in the meteorological data available. We therefore propose to use a griddistorteddataassimilation to correct for the position errors without attributing this error to a specific part of the simulation model or its input.

4.1 Necessity of position error correction

Position errors pose a great challenge for data assimilation, where it is often easier to adjust amplitudes rather than a position. This strongly limits the forecast skill, and further improvement requires the correction of position errors.

The difference between assimilation of observations with or without correction of position errors is illustrated in Fig. 4. The panels show a hypothetical dust concentration along a coordinate, which could be either spatial or temporal without loss of generality. The a priori simulation (dashed) differs from the observations (stars) both in amplitude and shape (location and width in space or arrival and duration in time). The underlying simulation model is therefore likely to be imperfect in either emission strength, emission timing or transport or a combination of all of these. The left panel illustrates a typical assimilation of observed concentrations that adjust emission strengths only. In such an assimilation, the a priori concentrations are just scaled towards the observations. The posterior concentrations are therefore closer to the observations but only where the a priori simulations have any concentrations at all. On the left side of the axis the simulated concentrations are therefore strongly reduced to match with the zero observations. However, if initially no dust is present in the simulations, as is the case on the right side of the axis, then the assimilation does not suddenly introduce dust out of nothing.

The right panel illustrates how a position error correction could improve this. Before analysing the observations, the a priori plume is shifted and reshaped to have the best match with the observations, ignoring differences in amplitude. If this repositioned plume is analysed with the available observations, the posterior result is in much better agreement with the observations along the entire axis, also where initially no dust was simulated. The assimilation will still adjust the emission strengths, but these are now not adjusted to correct for transport errors.

Figure 5 Illustration of grid distortion technique: (a) original grid map, (b) original dust concentrations as band, (c) distorted grid map, (d) distorted dust concentrations.

Table 1 Definition of assimilation experiments.

Figure 6 Diagrams of emis inversion, grid-distorted assim and hybrid assim systems.

4.2 Grid distortion

To align the dust plume with the observations, a grid distortion method as described by Lawniczak (2012) is used. The procedure is illustrated in Fig. 5. In transport models, the flow equations are usually solved on a discrete grid. For the LOTOS-EUROS model used here, the grid is Cartesian (perpendicular in longitude and latitude) and regular in spacing (panel a of the figure). Computed concentrations represent an average over a grid cell, and the simulated plume therefore consists of a set of grid cells with a substantial dust load. Panel (b) shows an example with a dust plume as a band from left to right. The grid distortion smoothly transforms the Cartesian grid into a non-Cartesian grid. That is, the corners of the grid cells are repositioned to a nearby location such that each distorted grid cell remains connected to its original neighbours (panel c). The dust concentration in each grid cell (in µg m⁻³) is kept constant after distortion to ensure a smooth variation in dust intensities over neighbouring cells. The dust plume is deformed together with the grid (panel d).

In mathematical formulation, let (x,y) denote the original Cartesian coordinates. A discrete model grid with regular spacing Δx×Δy is defined on points (x_i,y_j), with i and j the integer indices of the grid points in the x and y direction. The grid distortion is defined as a coordinate transformation that projects an original location (x,y) onto a new location (λ,ψ) with

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }\mathit{\lambda }=\mathrm{\Lambda }(x,y)\text{(3)}& {\displaystyle }& {\displaystyle }\mathit{\psi }=\mathrm{\Psi }(x,y).\end{array}$$

Following Lawniczak (2012), the grid distortion is described using a Poisson equation. The elliptic equation is broadly utilized in mechanical engineering and theoretical physics to describe how an object diffuses in space given a charge (Hazewinkel, 1994). The repositioned grid locations (λ,ψ) are the solutions of two 2D Poisson equations with the charges or distortion functions 𝒫 and 𝒬 on the right-hand side:

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle \frac{{\partial }^{\mathrm{2}}\mathrm{\Lambda }}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathrm{\Lambda }}{\partial y^{\mathrm{2}}}}=\mathcal{P}(x,y)\text{(5)}& {\displaystyle }& {\displaystyle \frac{{\partial }^{\mathrm{2}}\mathrm{\Psi }}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathrm{\Psi }}{\partial y^{\mathrm{2}}}}=\mathcal{Q}(x,y).\end{array}$$

The distortion functions 𝒫 and 𝒬 that drive the grid distortion are initially unknown, and their optimal values are to be calculated as part of the data assimilation procedure described in Sect. 4.3.

The second-order derivatives in Eqs. (4) and (5) are discretized on the grid using finite differences. For Eq. (4), the discretization is

$$\begin{array}{}\text{(6)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{{\mathit{\lambda }}_{i+\mathrm{1},j}-\mathrm{2}{\mathit{\lambda }}_{i,j}+{\mathit{\lambda }}_{i-\mathrm{1},j}}{(\mathrm{\Delta }x{)}^{\mathrm{2}}}}\\ & +{\displaystyle \frac{{\mathit{\lambda }}_{i,j+\mathrm{1}}-\mathrm{2}{\mathit{\lambda }}_{i,j}+{\mathit{\lambda }}_{i,j-\mathrm{1}}}{(\mathrm{\Delta }y{)}^{\mathrm{2}}}}=P_{i,j}\end{array},\end{array}$$

and Eq. (5) gives a similar discretization. When this system is solved for a given right-hand side, the result is a grid of 2D locations $({\mathit{\lambda }}_{i,j},{\mathit{\psi }}_{i,j})$ corresponding to the distorted positions of the original grid points (x_i,y_j). This system can be solved using a numerical method for linear equations. In our experiments, we use the red–black ordering Gauss–Seidel method (Saad, 2003) to solve the discrete system of linear equations.

The distorted dust plume is interpolated back to the Cartesian grids using the nearest searching method (Cayton, 2008) for comparison with observations (that are defined on longitude–latitude coordinates), and to serve as initial fields for the following simulation steps.

4.3 Distortion estimation using 4DEnVar

The grid distortion method provides a new way of repositioning the dust plume without adjusting the long-distance dust transport. We use the ensemble-based variational (4DEnVar) data assimilation (Liu et al., 2008b) algorithm to optimize the grid distortion.

To find the optimal distortion, the initial value and covariance of 𝒫 and 𝒬 need to be defined first. Each element in the two distortion equations is assumed to have a zero mean and a standard deviation, empirically chosen to be 0.015. To enforce a smooth grid distortion, we also prescribe a correlation c between two elements 𝒫(x_i,y_j) and 𝒫(x_k,y_l) (and similar for 𝒬):

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }c=e^{-d(x_i,y_j;x_k,y_l)/\mathcal{L}},\end{array}$$

where d represents the spatial distance in kilometres, and ℒ is an empirical length scale that is set to 1000 km. The parameters used in this study (standard deviation, correlation length scale) were chosen based on experiments for the described dust event; for other events they might need to be revised.

In our 3D model, the grid distortion is applied in the horizontal direction only, changing each layer in the same way. This is mainly to reduce the degrees of freedom in the distortion since no information on the 3D structure of the plume is available from the current observations (surface data and satellite-retrieved column information). It is however also possible to use a 3D distortion with a few degrees of freedom in the vertical (Nehrkorn et al., 2015) for dust events where measurements of the vertical structure are available, e.g. lidar backscatter coefficient (Madonna et al., 2015).

An ensemble of random distortion fields is generated using the assumed prior value (zero) and the assumed covariance. Each member is a vector s collecting all elements of 𝒫 and 𝒬 on the discrete grid:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }[{\mathit{s}}_{\mathrm{1}},\mathrm{\dots },{\mathit{s}}_N].\end{array}$$

In our experiments the ensemble size N was set to 100. For each of these ensemble members, the distorted grid (λ,ψ) is solved from the system of the discrete Poisson equations as described in Sect. 4.2. With this an ensemble of distorted dust maps is formed from the a priori dust field x:

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }\left[\mathit{x}\right({\mathit{s}}_{\mathrm{1}}),\mathrm{\dots },\mathit{x}({\mathit{s}}_N\left)\right],\end{array}$$

where x(s_i) represents the distorted dust field using distortion s_i.

Denote the ensemble perturbation matrix or covariance square root by

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }{\mathbf{S}}^{\prime }={\displaystyle \frac{\mathrm{1}}{\sqrt{N-\mathrm{1}}}}[{\mathit{s}}_{\mathrm{1}}-{\mathit{s}}_b,\mathrm{\dots },{\mathit{s}}_N-{\mathit{s}}_b],\end{array}$$

where s_b is the (zero) prior value. In a 4DEnVar assimilation, the optimal distortion vector s_a is defined to be a weighted sum of the columns of the perturbation matrix S^′ using weights from a control variable vector w:

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }{\mathit{s}}_a={\mathit{s}}_b+{\mathbf{S}}^{\prime }\mathit{w}.\end{array}$$

The optimal control variables are then calculated through minimizing of the cost function:

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }J\left(\mathit{w}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{w}}^T\mathit{w}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathbf{HXS}}_b^{\prime }\mathit{w}+\mathit{d}{)}^T{\mathbf{R}}^{-\mathrm{1}}({\mathbf{HXS}}_b^{\prime }\mathit{w}+\mathit{d}).\end{array}$$

In here, d is referred to as the innovation that describes the difference between observations y and simulations on the distorted grid:

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }\mathit{d}=\mathbf{H}\mathit{x}\left({\mathit{s}}_b\right)-\mathit{y}.\end{array}$$

In here, H is the observation operator that simulates the observed value on the distorted grid, which here simply takes the model simulation from the grid cell holding the observation location. The distortion uncertainty is transferred into the observation space through application of H on the ensemble members:

$$\begin{array}{}\text{(14)}& {\displaystyle }\begin{array}{rl}& {\mathbf{HXS}}_b^{\prime }\approx {\displaystyle \frac{\mathrm{1}}{\sqrt{N-\mathrm{1}}}}\left[\mathbf{H}\mathit{x}\right({\mathit{s}}_{\mathrm{1}})-\mathbf{H}\mathit{x}({\mathit{s}}_b),\mathrm{\dots },\mathbf{H}\mathit{x}({\mathit{s}}_N)\\ & -\mathbf{H}\mathit{x}\left({\mathit{s}}_b\right)].\end{array}\end{array}$$

The observation error covariance matrix R describes the possible differences between simulations and observations due to observation representation errors. R here is defined as a diagonal matrix, in which each representation error is set to an observation-dependent value ranging from 100 to 200 µg m⁻³ following Jin et al. (2018).

To ensure that the position correction is not too much influenced by differences in dust intensity, both the observations y and prior dust simulations x are normalized using their maximum values. Elements in R are also then scaled using the square of the maximum observed value.

The computation of the N=100 grid distortions is the most time-consuming part of the 4DEnVa-based grid-distorted data assimilation method; each of them costs around 2 min in our computing platform (CPU: Intel Xeon(R) E5; programming language: Python 3.7.6). The computation of the ensemble distortions could be re-implemented in a more efficient language but also be easily parallelized; the grid-distorted assimilation method is therefore expected to be computationally efficient enough to allow implementation in an operational forecast.

5 Dust storm data assimilation

The grid-distorted data assimilation was introduced for repositioning the simulated dust clouds. To evaluate the effectiveness, assimilation experiments including grid distortion have been performed and compared with a traditional assimilation configuration focusing on intensities only and a hybrid assimilation that combines these two.

An a priori simulation serves as a reference for all assimilation experiments. The emission inversion assimilation corrects for the dust intensity errors only, while the grid-distorted assimilation only corrects for the position error. The hybrid assimilation combines both in order to correct for the intensity as well as the position error.

5.1 Assimilation methods

Figure 6 shows the schematic overview of the three assimilation methods listed in Table 1. The left panel shows the set-up of the emission inversion, as described in detail in Jin et al. (2019b, 2020). The inversion combines the transport model (LOTOS-EUROS) with a four-dimensional variation (4DVar) data assimilation using a reduced-tangent linearization (Jin et al., 2018). The system assumes that the processes of dust transport and removal are simulated correctly, while only the emission is imperfect. The uncertainty in the emissions was parameterized as a sum of two sources: the uncertainty in the friction velocity threshold and in the erosive wind fields. The dust emission intensity in the source regions is then optimized such that the amplitude of the simulated concentrations is as close to the observations as possible. The optimized emission fields could then be used to drive simulations that have a better forecast skill than simulations with the original emissions.

The grid-distorted assim is designed to adjust the position of the simulated dust plume only. As described in Sect. 4.3, the impact of the actual dust concentrations is avoided by normalizing the dust simulations and observations using their maximum values before calculation of the distortion; afterwards, the distorted dust field is multiplied by the same maximum value again.

The right panel of Fig. 6 shows the set-up of the hybrid assim. Different from the emis inversion and grid-distorted assim, the hybrid assim performs two assimilations sequentially. First the grid-distorted assim is conducted for repositioning the simulated dust plume. Then, the position-corrected dust plume is used as a prior in the second assimilation (similar to an emis inversion) to adjust the emissions to have the best possible match between actual (not normalized) observations and position-corrected simulations. The posterior dust field from the hybrid assim is then used as the initial condition for forecast simulations.

In all assimilation tests, only observations from the snapshot of 5 May at 15:00 CST are used for fair comparison. The repositioned plume is only available for this single moment; measurements at earlier times can therefore not be accurately assimilated in hybrid assim since the corresponding simulation still has a position error then. In the emis inversion, the assimilation window is set from 2 May, 08:00 CST, which fully covered the related dust emission for this event.

5.2 Optimized plume position and dust load

The a priori dust plume described in Sect. 3.1 is assimilated with observations using the emis inversion, the grid-distorted assim or the hybrid assim. The posterior surface concentrations are shown in Fig. 2d–f, respectively. The optimized dust plumes are evaluated by their position and the RMSE metric that was introduced in Sect. 3.1 to quantify the difference with observations. Note that observations that are used to evaluate the posterior performance are the same as those that have been assimilated. When evaluating the method over a longer time period (multiple dust events), validation with independent observations should be considered.

Panel (d) shows the posterior dust plume using the emis inversion. The markers indicate that it has in general the same position as the a priori, and hence the position error has not been corrected yet. In terms of root mean square error (RMSE), the emis inversion posterior simulation is improved, but only slightly; the RMSE is reduced from 388.1 µg m⁻³ for the a priori simulation to 362.9 µg m⁻³. The emis inversion also has little effect on the dust simulation at earlier periods of the dust event, which can be found through a comparison of the a priori and emission-inversion-only simulations on 3 May at 13:00 CST in Fig. S1 in the Supplement. The a priori and emis inversion also present a relatively similar performance in the early stage.

Using the grid-distorted assim, the repositioned dust plume in panel (e) matches well with the ground observations shown in panel (a.2). The marker indicating the left side of the plume is now around 35^∘ N, which is in agreement with the observations; also the markers at the centre and the right side are now better positioned. Only the very left part of the repositioned dust plume (west of 110^∘ E) still shows a discrepancy compared to the PM₁₀ observations. This can be explained from the fact that this part of the dust plume has a relatively low dust load, which makes the corresponding position error less important in the cost function Eq. (12). In addition, a rather large grid distortion is required for this part of the dust plume to match the measurements, which is constrained with the assumed covariance of the distortion function. The RMSE of the posterior simulation is now significantly reduced to 251.1 µg m⁻³. Though the dust plume is now correctly repositioned, the simulated dust concentration does not exactly match the actual measurements. Especially in the plume centre, the posterior simulation shows dust concentrations over 1200 µg m⁻³ that are still similar to the a priori simulation, while the bias-corrected observations indicate that the dust intensity at most stations is lower than 1200 µg m⁻³.

The hybrid assim posterior simulation provides the best performance, as shown in panel (f). Not only is the dust plume realigned with the observations, but also the amplitude of the dust loads agrees better with the actual situation. For instance, the dust concentration in the plume centre is reduced from 1500 to 1200 µg m⁻³, and in the upper left part of the plume the concentration level is lifted from 100 to 200 µg m⁻³. As a result, the RMSE in the hybrid assim is reduced to 223.4 µg m⁻³.

Figure 7 PM₁₀ (a.1), bias-corrected dust observations (a.2), the a priori forecast (b), ensemble maximum (c) dust forecast and dust forecast driven by emis inversion posterior (d), the grid-distorted assim posterior (e), and hybrid assim posterior simulation (f) at 21:00 CST on 5 May. SDC: surface dust concentration.

Figure 8 Time series of PM₁₀ measurements, bias-corrected dust observations, a priori simulation and forecasting driven by the initial state from emis inversion, grid-distorted assim and hybrid assim at Baoji (a), Xinxiang (b), Xuchang (c), Bozhou (d), Zibo (e) and Zhaoyuan (f). The vertical dashed black line indicates the start of the forecast.

Figure 9 RMSE of the a priori dust simulation and forecasting using the initial state from emis inversion, grid-distorted assim and hybrid assim.

5.3 Forecasting of dust plume position

In an operational setting the posterior dust concentrations are used as initial conditions for a forecast. Starting from the analysis results, forecast runs have been performed. A snapshot of the resulting forecast of the surface dust concentrations as well as the PM₁₀ measurements, bias-corrected dust observations and the a priori forecast at 21:00 CST on 5 May are shown in Fig. 7.

The ground observations in panel a.1 and a.2 indicate that the dust plume is now located along 35^∘ N. In both the a priori simulation and in the forecast based on emis inversion, the right, centre and left plume markers are about 100, 300 and 200 km farther south, respectively. However, the forecasts based on grid-distorted assim or hybrid assim assimilation both show plumes with positions in better agreement with the observations. The best results are obtained for the hybrid assim, which shows better agreement for the central and upper right part of the dust field (panel f) compared to the grid-distorted result (panel e).

5.4 Time series at stations

Figure 8 shows times series of dust concentrations at six different observation sites. The locations can be found in Fig. 1 and were selected to illustrate the general results but also challenges to be solved in future. The time series show PM₁₀ observations (red circles), bias-corrected observations representing the dust part (red dot), the a priori forecast (black line) and the forecasts driven by the three assimilation tests starting from 15:00 CST on 5 May.

For all six sites, the a priori dust simulations estimate an arrival time of the dust cloud that is at least 4 h too early. The emis inversion focusing on intensity correction does not improve the forecasting of the arrival time since it only changes the emission strength. Ignoring the intensity of the dust load, the temporal profiles of the dust forecasts driven by the grid-distorted assim after 5 May at 15:00 CST are in good agreement with the temporal profile of the dust observations.

For stations on the upper side of the plume, e.g. Baoji in panel (a), the declining trend predicted by the a priori and emis inversion forecasts is well reproduced by the grid-distorted assim. For sites where the descent pattern was not captured by the a priori simulation, the emis inversion helps little, while grid-distorted assim resolves the decreasing trend, as can be seen in Zibo and Zhaoyuan. For stations downwind of the plume like Xinxiang, Xuchang and Bozhou, the dust concentrations show an up-and-down pattern caused by the arrival and departure of the plume. The a priori and emis inversion forecasts are unable to capture the dust profile. For instance in Bozhou, the a priori simulation indicated that the main dust plume arrived earlier than 00:00 CST on 5 May, and it started to decline from 12:00 CST. However, the real observation showed that the dust storm actually arrived around 12:00 CST, with a steady increase in concentration. Starting from the grid-distorted assim, the forecast shows concentrations with a trend similar to the observations, although the increase starts a few hours too early. The observation-minus-simulation discrepancy is further reduced for most stations using the hybrid assim that combines the grid-distorted assim and emis inversion.

5.5 Evaluation of forecast skills

The forecast skill of the three assimilation algorithms is also evaluated using the RMSE indicator that was also used for the a priori and posterior dust simulations in Sects. 3.1 and 5.2.

During the period from 16:00 CST on 5 May to 07:00 CST on 6 May, the a priori RMSE reached values around 300 µg m⁻³. The assimilation based on emis inversion helped to decrease the RMSE of the forecast simulations with about 20 µg m⁻³. The improvement is limited since position errors are dominant and still present. The grid-distorted assim is efficient in enhancing dust forecast skills in terms of the RMSE, which significantly reduce to less than 200 µg m⁻³. When combined with emis inversion in the hybrid approach, an additional decrease in RMSE of about 20 µg m⁻³ is achieved.

These results show that the grid-distorted assim is capable of correcting the position error in the simulated dust plume effectively; the hybrid assim that combines the grid-distorted assim and emis inversion provides the best initial condition to drive the dust forecast in the short term.

6 Summary and conclusions

Evaluation of dust storm forecasts focuses on two main criteria: the intensity of the dust load and the position of the cloud. Various studies on improving dust forecasts focused mainly on correcting the intensity only. However, positional misfits are unavoidable as a result of inaccurate emission timing profile and/or accumulation of uncertainties in long-distance transport and therefore need to be taken into account too.

An extremely severe dust storm in May 2017 over East Asia was used as the test case in this study. A regional chemical transport model, LOTOS-EUROS, was used to reproduce the dust event. PM₁₀ observations are available from the China Ministry of Environmental Protection (MEP) air quality monitoring network; bias-correction was used to process the original PM₁₀ measurements to accurately represent the dust load. The position misfits are obviously detected in the results, especially when the simulated dust plume is transported thousands of kilometres away to central China.

The positional misfit in dust simulation could be corrected by data assimilation too. Assimilation configurations for this type of application usually require definition of a background error covariance or an ensemble perturbation scheme that could resolve the full observation–simulation positional discrepancy too. This covariance could for example include the meteorological uncertainty, as described by a meteorological ensemble forecast. For the dust storm studied here it was however shown that the spread in the available meteorological ensemble and/or the way in which the simulation model is using it is not sufficient to explain the position error in the simulations. Therefore, a complex covariance model that could account for the accumulation of uncertainties along the long track of the plume would be required while using the traditional assimilation method. Meanwhile, a substantial number of measurements would then be required to constrain the optimal transport pattern too.

Alternatively, an image-morphing method, grid distortion, is adopted to reposition the simulated dust plume in this paper. The method is then combined with 4DEnVar for a griddistorteddataassimilation, which focuses solely on correcting the dust field position to best fit the assimilated observations. Since in reality both position and intensity errors might be present, a hybrid assimilation algorithm is proposed. In this hybrid system, the griddistorteddataassimilation and an intensity-centred emission inversion are performed after each other.

Assimilation tests using either the emission inversion or griddistorteddataassimilation only or using the hybrid assimilation have been conducted on the studied dust event. The posterior dust simulation and the forecast are slightly improved by using emission inversion. This indicates that the traditional intensity-centred data assimilation is of little help in the case that positional errors are present. Only using the griddistorteddataassimilation, strongly improved posterior and forecast simulations are obtained. The best results are obtained when the hybrid assimilation is performed, with both the position and intensity errors corrected.

The grid-distorted assimilation should be seen as an extension to traditional intensity-centred assimilation, not as a replacement. In the presence of a position error, grid-distorted data assimilation is a computationally efficient pre-processing procedure to correct for errors that are not resolved otherwise. The method could be used to further explore 3D dust and aerosol structure by combining the 3D grid distortion and observations with vertical layering information.