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.

2007; Wu et al., 2016). Apart from the influence on the environment, dust storms pose a great threat on the human health by carrying thousands tons 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 5 traffic. For example, more than 1,100 flights were delayed/canceled 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 modeling systems (Marticorena and Bergametti, 1995;Shao et al., 1996;Marticorena et al., 1997;Alfaro et al., 1997;10 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 on society. Though these systems are usually able to predict the starting and 15 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 up to two orders of magnitudes (Uno et al., 2006;Gong and Zhang, 2008). The modeling 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;20 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 25 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., 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), as well as the satellite retrieved AOD (Khade et al., 2013;Yumimoto et al., 2016;Di Tomaso et al., 2017;Dai et al., 2019). Those studies 30 either focused 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 10 (total particulate matter with diameter less then 10 µm) measurements (Jin et al., 2018(Jin et al., , 2019a and geostationary satellite AOD (Jin et al., 2019b(Jin et al., , 2020 were assimilated with the LOTOS-EUROS simulation model for dust storm forecasts over East Asia. Also these studies soleley 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 5 improve the simulated dust concentrations at time of observation, but also to improve forecasts of dust levels over windows of up to 24 hours. During these studies it was noted that although the modeling system in general provided an accurate forecast of the dust plume, a severe position error was present when the plume traveled over a large distance. Specifically, forecasts by the model simulation reported the dust arrival and departure 1 to 10 hours prior to reality, as will also be illustrated in Section 3.
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. 15 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 will be illustrated in Section Similar as 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 25 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 area. Position discrepancies might also arise from incorrect timing profiles of emissions, which is not the case for our test event as will be explained in Section 3.1. However, determining the covariance either for transport or for emission timing 30 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 amount 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 EnKF for wildfire models (Beezley and Mandel, 2008), grid distortion data assimilation on oil reservoir modeling (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 re-position the 5 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, 10 which is a novel approach in the context of atmospheric dust modeling. 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 that both position and intensity misfits are presented in a dust simulation. Starting from the initial condition using the hybrid assimilation 15 posterior, dust forecast 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 Section 4, the necessities of position error correction is emphasized first, and then 20 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 Section 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 25

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 30 2020). Besides, it is also implemented in the WMO Sand and Dust Storm Warning Advisory and Assessment System to provide short-time forecast of the dust load over the North Africa-Middle East-Europe areas, the online forecast product are 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°N to 50°N and 70°E to 140°E, with a resolution about 0.25°× 0.25°. Vertically, the model consists of 8 layers with a top at 10 km. The dust simulation is driven by European Center for Medium-Ranged Weather Forecast (ECMWF) operational forecasts over 3-12 hours, retrieved at a regular longitude/latitude grid resolution of about 7 km. Physical processes included are wind-blown dust emission, diffusion, advection, dry and wet deposition, and sedimentation.

Observation network
The observations used in this study consist of hourly PM 10 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 10 measurements are actually a sum of dust and airborne particles (black carbon, sulphate, etc). Since the 10 analyzed event is an extremely severe case, these PM 10 measurements were directly used to quantify the dust load in Jin et al. (2019bJin et al. ( , 2020. In this study however, an observational bias correction is performed to make the PM 10 measurements fully representative for 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 10 observations. The original PM 10 measurements vs. mation that has been widely investigated already, the position error has received less attention, but it has been the main focus of this study.

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. (2019bJin et al. ( , 2020. The dust emission occurred   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/m3. This vast mismatch is attributed to the sum of intensity and position error (mainly) as will be explained in Section. 5.2.
3.2 Uncertainty in emission timing profile 25 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 hours after dust emission started, the simulated dust plume was still in north China and showed in general the same pattern as 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 north China (Jin et al., 2019b).

30
The good phase match in general can also be seen from a snapshot of the ground PM 10 observation vs. the simulated surface dust concentration at May 04 15:00 (CST) in Fig. 3. There might already be position misfits in the dust simulation at these snapshots, but not easily detected. The magnitudes of the dust concentration showed discrepancies, but these could be corrected

Uncertainty in meteorology
Another possible origin of the position error in the simulations is the uncertainty in the meteorological data. In our study, the 5 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.
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 10 ensemble via: In here, c i represents the dust concentration field that results from a simulation with the i th ensemble member. This measure reflects for each location whether in any of the simulations a severe dust load is present. A snapshot of the ensemble maximum Eq.
(1) at 15:00 is shown in Fig. 2 (c). The map shows a broader plume, which implies that some ensemble members result in a 15 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 in the meteorological data therefore could not be used to fully account for the position error, and the assimilation system has to correct for this in a different way.
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 grid distorted data assimilation to correct for the position errors, without attributing this error to a specific part of the simulation model or its input. 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 concentrations along a coordinate, which could be either spatial or temporal without loss of 10 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 strengths, 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 15 closer to the observations, but only where the a priori simulations has any concentrations at all. At 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 at 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 analyzing the observations, the a 20 priori plume is shifted and reshaped to have the best match with the observations, ignoring differences in amplitude. If this re-positioned plume is analyzed 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.

Grid distortion 25
To align the dust plume with the observations, a grid distortion method as described by Lawniczak (2012)  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 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: https://doi.org/10.5194/gmd-2021-10 Preprint. Discussion started: 9 March 2021 c Author(s) 2021. CC BY 4.0 License.
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 re-positioned grid locations (λ, ψ) are the solutions of two 2D Poisson equations with at the right hand side the charges or distortion functions P and Q: The distortion functions P and Q 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 Section 4.3.
The second-order derivatives in Eq. (4) and (5) are discretized on the grid using finite differences. For Eq. (4), the discretization is: and a similar discretization for Eq. (5). When this system is solved for a given right hand side, the result is a grid of 2D locations (λ i,j , ψ 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. 15 The distorted dust plume is interpolated back to the Cartesian grids using nearest searching method (Cayton, 2008) for comparison with observations (that are defined on longitude/latitude coordinates), and to serve as initial fields for following simulation steps.

Distortion estimation using 4DEnVar
The grid distortion method provides a new way for re-positioning the dust plume without adjusting the long-distance dust 20 transport. We use the ensemble-based variational (4DEnVar) data assimilation  algorithm to optimize the grid distortion.
To find the optimal distortion, the initial value and covariance of P and Q 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 P(x i , y j ) and P(x k , y l ) (and similar for Q): where d represents the spatial distance in km, and L 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 the 3D model, the grid distortion is applied in horizontal direction only, changing each layer in the same way. This is 30 mainly to reduce the degrees of freedom in the distortion, since no information on the 3D structure of the plume is available from the 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).
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 P and Q on the discrete grid: 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 section 4.2. With this an ensemble of distorted dust maps is formed from the a priori dust field x: where x(s i ) represents the distorted dust field using distortion s i .

10
Denote the ensemble perturbation matrix or covariance square root by: 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: The optimal control variables are then calculated through minimizing of the cost function: In here, d is referred to as the innovation that describes the difference between observations y and simulations on the distorted grid: 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: The observation representation error R describes the possible differences between simulations and observations due to model 25 and observation errors, and is here defined 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 scaled using the square of the maximum observed value.

Dust storm data assimilation
The grid distorted data assimilation was introduced for re-positioning the simulated dust clouds. To evaluate the effectiveness, assimilation experiments including grid distortion have been performed and compared with a traditional assimilation focusing on intensities only and a hybrid assimilation that combines these two. An a priori simulation serves as reference for all assimilation experiments. The emission inversion assimilation corrects for the dust intensity errors only, while the grid distorted 5 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. Fig. 6 shows the schematic overview of the three assimilation methods listed in Table. 1. The left panel shows the setup of the emission inversion, as described in detail in (Jin et al., 2019b(Jin et al., , 2020. The inversion combines the transport model 10 (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 emissions intensity in the source regions is then nudged to make simulation and observations in harmony. The optimized emission fields could then be used to drive simulations that have a better forecast 15 skill than simulations with the original emissions.

Assimilation methods
The grid distorted assim is designed to adjust the position of the simulated dust plume only. As described in Section. 4.3, the impact of the actual dust concentrations is avoided by normalizing the dust simulations and observations using their maximum values before calculation the distortion; afterwards, the distorted dust field is multiplied with the same maximum value again.
The right panel of Fig. 6 shows the setup of the hybrid assim. Different form the emis inversion and grid distorted assim, 20 the hybrid assim performs two assimilations sequentially. First the grid distorted assim will be conducted for re-positioning the simulated dust plume. Then, the position-corrected dust plume will be used as 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.

Optimized plume position and dust load
The a priori dust plume described in Section 3.1 is assimilated with observations using the emis inversion, or the grid distorted 5 assim, or the hybrid assim. The posterior surface concentrations are shown in Fig. 2 (d)~(f), respectively. The optimized dust plumes are evaluated by their position, and the RMSE metric that was introduced in Section. 3.1 to quantify the difference with observations. 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 10 emis inversion posterior simulation is improved but slightly; the RMSE is reduced from 388.1 µg/m 3 for the a priori simulation to 362.9 µg/m 3 .
Using the grid distorted assim, the re-positioned 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 center and the right side are now better positioned. Only the very left part of the re-positioned dust plume (west of 110°E) still shows a discrepancy compared to the PM 10 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 5 simulation is now significantly reduced to 251.1 µg/m 3 . Though the dust plume is now correctly re-positioned, the simulated dust concentration does not exactly match the actual measurements. Especially in the plume center, the posterior simulation show dust concentrations over 1200 µg/m 3 that are still similar to the a priori simulation; while the bias-corrected observations indicate that the dust intensity in most stations are lower than 1200 µg/m 3 .
The hybrid assim posterior simulation provides the best performance as shown in panel (f). The dust plume is not only 10 re-aligned 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 center is reduced from 1500 µg/m 3 to 1200 µg/m 3 , and in the upper left part of the plume the concentration level is lifted from 100 to 200 µg/m 3 . As a result, the RMSE in the hybrid assim is reduced to 223.4 µg/m 3 .

Forecast of dust plume position
In an operational setting the posterior dust concentrations are used as initial conditions for a forecast. Starting from the analysis

Evaluation of forecast skills
The forecast skill of the three assimilation algorithms are also evaluated using the RMSE indicator that was also used for the a priori and posterior dust simulations in Section. 3.1 and Section. 5.2.
During the period from 16:00 May 05 to 07:00 May 06, the a priori RMSE reached values around 300 µg/m 3 . The assimilation based on emis inversion helped to decrease the RMSE of the forecast simulations with about 20 µg/m 3 . The improvement 5 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 3 . When combined with emis inversion in the hybride approach, an additional decrease in RMSE of about 20 µg/m 3 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 10 drive the dust forecast in short term.
assimilation analysis forecast observation missing Figure 9. RMSE of the a priori dust simulation, and forecast using initial state from emis inversion, grid distorted assim and hybrid assim.

Summary and conclusion
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, 15 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 10 observations are available from the China Ministry of Environmental Protection(MEP) air quality monitoring network; bias-correction was used to process the original PM 10 measurements to accurately represent the dust load. The position misfits are obviously detected in the results especially 20 when the simulated dust plume is transported thousands kilometers away to center China.
The positional misfit in dust simulation could be corrected by data assimilation too. Traditional assimilation approached require definition of a background error covariance that should account for the observation/simulation positional discrepancy. 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 meteorological conditions is not sufficient to explain the position error in the simulations.

5
In this paper, an imaging morphing method, grid distorted, is adopted to re-position the simulated dust plume. The method is then combined with 4DEnVar for a grid distorted data assimilation, 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 at present, a hybrid assimilation algorithm is proposed. In this hybrid system, the grid distorted data assimilation and a intensity-centered emission inversion are performed after each other. 10 Assimilation tests using either the emission inversion or grid distorted data assimilation 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-centered data assimilation is of little help in case positional errors are present. Only using the grid distorted data assimilation, 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 15 errors are corrected.
The grid distorted assimilation should be seen as an extension to traditional intensity-centered assimilation, not as a replacement. Under the presence of position error, applying grid distortion should be a pre-processing procedure to correct for errors that are not resolved otherwise.
Code and data availability 20 The source code and user guidance of the CTM, Lotos-Euros, can be obtained from https://lotos-euros.tno.nl. The grid distorted data assimilation algorithm is in python environment, and is archived on Zenodo (https://doi.org/10.5281/zenodo.4579960).

Author contribution
JJ and AS conceived the study and designed the grid distorted data assimilation. JJ and AS performed the control and assimilation tests and carried out the data analysis. AS, HL, HXL, BH, XW, and AH provided useful comments on the paper. JJ prepared the manuscript with contributions from AS and all others co-authors.