New strategies for vertical transport in chemistry-transport models: application to the case of the Mount Etna eruption on March 18, 2012 with CHIMERE v2017r4

Excessive numerical diffusion is one of the major limitations in the representation of long-range transport by chemistry-transport models. In the present study, we focus on excessive diffusion in the vertical direction, which has been shown to be a major issue, and we explore three possible ways to address this problem: increase vertical resolution, use an advection scheme with antidiffusive properties, and represent more accurately the vertical wind. This study is done with the CHIMERE chemistry-transport model, for the March 18, 2012 eruption of Mount Etna, which has released about 3 kt of sul5 phur dioxide in the atmosphere into a plume that has been observed by satellite instruments (IASI and OMI) for several days. The change from the classical Van Leer (1977) scheme to the Després and Lagoutière (1999) antidiffusive scheme in the vertical direction has been shown to bring the largest improvement to model outputs in terms of preserving the thin plume emitted by the volcano. To a lesser extent, improved representation of the vertical wind field has also been shown to reduce plume dispersion. Both these changes help reducing vertical diffusion in the model as much as a brute-force approach (increasing 10 vertical resolution).


Introduction
Among many other uses such as operational forecast of air quality, chemistry-transport models (CTM) have been used successfully in the past to represent long-range transport of polluted plumes from different types of sources (mineral dust, volcanic eruptions, biomass burning etc.). The CHIMERE 1 CTM (Mailler et al., 2017) has previously been used to assess Eyjafjalla- 15 jökull's 2010 eruption possible impact on air quality (Colette et al., 2011). Ash transport from this eruption has been modelled and compared to LiDAR vertical profiles, showing that the CHIMERE model represented correctly the advection of this volcanic plume from its source in Iceland to the LiDAR (Light Detection And Ranging) facility located in Palaiseau (France), thousands of kilometers away. The altitude, location and timing of the plume was represented correctly, but the authors have shown that their simulation presented a strongly overestimated vertical spread of the plume. Similar studies focusing on volcanic plume dispersion from Boichu et al. (2013) have highlighted overestimations of plume diffusion on 2010 Eyjafjallajökull.
Several parameters can influence the evolution of the modelled plume: the emission fluxes and time profile; the injection height and vertical profile; chemical processes involving the considered species; the wind field; the numerical advection schemes and 5 the vertical resolution. Boichu et al. (2015) have focused on volcanic plume dispersion sensitivity to altitude of injection, combining CHIMERE and Atmospheric Sounding Interferometer instrument (IASI) for a Mount Etna -study case of an eruption of moderate intensity in April 2011. The eruption presented an emission profile centered at 7 000 m.a.s.l., with weaker emissions at 4 000 m.a.s.l.. These authors found a strong sensitivity of model outputs to the altitude of injection. In Mailler et al. (2017), advection of the volcanic SO 2 plume emitted by a major eruption of the Puyehue Cordó Caulle volcano (Chile) is simulated 10 with the CHIMERE model and compared to satellite measurements and to analyses provided by Klüser et al. (2013). This work shows that, after about one week of simulation, circumpolar transport of the plume has been represented correctly and the final position of the leading edge of the plume is simulated in a reasonable way, but that plume dilution is excessive compared to the observed shape and concentration of the plume.
Most CTMs have been built as offline models forced by meteorological fields, in particular wind and air density, taken from 15 a forcing meteorological model, typically global forecast data such as outputs from the IFS (Integrated Forecasting System) or GFS (Global Forecast System), data from operational forecast centers or data generated by the modellers themselves from a locally run meteorological model. These meteorological data, after interpolation in time and space onto the CTM grid, are used to drive advection within the CTM. However, grid type, grid structure, transport schemes and time discretization are generally different in the CTM from their formulations in the forcing models, deriving into mass-wind inconsistencies: once 20 interpolated onto the CTM grid, the mass and wind fields do not obey the continuity equation anymore (Jöckel et al., 2001).
While theoretical pathways to mitigation of this problem exist (Jöckel et al., 2001), this problem has historically been solved practically in regional CTMs in a straightforward way as described in Emery et al. (2011): reconstructing the vertical wind from the density field and the horizontal mass flux divergence in order to artificially enforce verification of the continuity equation to the expense of the realism of the vertical mass fluxes, in particular in the free troposphere. This approach is justified by the fact 25 that, in the lowest atmospheric layers, the reconstructed vertical mass flux is not very different from the real mass flux from the forcing model, and that this explicitly resolved vertical transport is usually dominated by mixing inside the Planetary Boundary Layer (PBL). Therefore, it has long been thought that this approach generates little if any problem since the main purpose of regional CTMs is to provide a reliable forecast of the concentration of pollutants within the PBL. This historic focus on the PBL has also led to the habit of chemistry-transport modellers to use very loose vertical resolution in the free troposphere. 30 Emery et al. (2011) describes the side effects on such an approach in two of the most used CTMs (CAMx and CMAQ). They show that oversimplification of transport processes in the free troposphere can not only be detrimental to the representation of transport in the free troposphere itself but also defeats its own purpose: focus on obtaining a correct representation of pollutant concentrations in the PBL. They have shown that oversimplified representation of vertical transport and vertical mass fluxes in the free troposphere spuriously increases vertical transport of stratospheric ozone into the troposphere, resulting into degraded scores for forecast and analysis of ozone concentration in the PBL, particularly over complex and elevated terrain (springtime in the United States, in the case of Emery et al., 2011). To solve this problem, the authors tried different approaches. While trying to reduce spurious vertical velocities by applying mass filters, smoothers/desmoother filters or divergence minimizers to the forcing velocity field either brings little improvement to the issue or introduces numerical artefacts, improvement of the vertical transport scheme and increase in the vertical number of layers in the free troposphere did bring substantial improvement 5 to the issue of excessive transport of stratospheric ozone into the stratosphere.
Apart from this wind-mass inconsistency issue, and more specifically for the representation of polluted plumes that are transported over a long range, Zhuang et al. (2018) have shown that correct representation of long-range transport of polluted plumes in the free troposphere is severely limited by the insufficient vertical resolution. They show, through dimensional and theoretical arguments, that if ∆x is not at least several hundred times ∆z, representation of long-range transport of plumes 10 in the free troposphere is hindered primarily by this coarse vertical resolution, and increasing horizontal resolution does not bring substantial added value in terms of reducing numerical diffusion of the plume. Since the ∆x ∆z in typical chemistry-transport models is around or below 20 (with a horizontal resolution of, e.g., 20 km for continental scale studies and vertical resolution of, e.g., 1 km), these authors claim that no major improvement will be reached in the representation of long-range transport plumes unless vertical resolution is refined drastically compared to current typical configurations. However, they do not examine the 15 use of anti-dissipative transport schemes, which can be a possibility to reduce vertical diffusion without dramatically increasing vertical resolution.
In the present study, three options have been tested in terms of the accuracy of representation of the long-range advection of thin layers. One option is to choose the Després and Lagoutière (1999) anti-diffusive advection scheme for vertical transport, the second option is to use realistic vertical mass fluxes instead of reconstructed vertical mass fluxes, the third option being 20 refinement of vertical resolution. We have chosen the March 18, 2012 Etna volcano eruption to evaluate the impact of these new strategies for vertical transport. This volcano is well monitored, thus allowing to gather detailed model inputs and correlative data. This eruption has been relatively strong, so that the resulting volcanic plume has been distinctly observed and followed by satellite instruments, permitting comparison of modelled to observed plume at different stages of plume evolution over more that two days. Etna volcanic activity is monitored continuously to estimate SO 2 fluxes and plume injection height (Salerno 25 et al., 2009;Mastin et al., 2009;Sellitto et al., 2016;Salerno et al., 2018). Volcanic gas and aerosol are also subject to numerous physical and chemical evolution processes, such as sulphates production or CCN activation, whose have been recently studied (Sellitto et al., 2017;Guermazi et al., 2019;Pianezze et al., 2019) but are not accounted for in this study.
The manuscript is structured in the following way: Material and method (section 2) presents the CHIMERE model configuration for these simulations, including a detailed presentation of the transport formulation and its discretization in the 30 CHIMERE model, adaptation of the Després and Lagoutière (1999) scheme to the CHIMERE framework and presentation of the method for compensation of mass-wind inconsistencies that permits us to use realistic vertical mass fluxes instead of reconstructed mass fluxes. Also in this section, we discuss the satellite data that we used as a comparison point for our model outputs, the different settings of the performed sensitivity tests and the SO 2 emission fluxes that we use. In the results and discussion (section 3), eruption injection altitude impact on plume transport is first investigated and compared to the plume transport constructed from satellite based instruments. In addition, sensitivity to vertical profile of injection has been evaluated.
Then, the dispersion and trajectory of the simulated plumes is discussed and compared to the available data, with a focus on the impact of the various tested parameters on plume dispersion (vertical wind representation, vertical advection scheme and number of vertical levels).
2 Material and methods 5

CHIMERE simulations
Simulations have been performed using a development version of the CHIMERE model Mailler et al., 2017) including the new developments presented in section 2.2. The simulations have been performed with no chemistry, and an inert gaseous tracer with the molar mass of SO 2 has been emitted at the location of the Etna volcano, with fluxes and injection heights that are presented below. In particular, oxidation of SO 2 and subsequent formation of sulphate or sulfuric acid 10 have not been represented. Simulations last 120 hours starting on March 18, 00 UTC. The CHIMERE model has been forced using WRFv.3.7.1 (Weather Research and Forecasting Skamarock et al., 2008), with an update of the forcing meteorological variables every 20 minutes using the WRF-CHIMERE online simulation framework (Briant et al., 2017). The WRF model has been run with 33 vertical levels from surface to 55 hPa. Boundary conditions meteorological fields and spectral nudging of the WRF simulation have been taken from the NCEP GFS dataset at 0.25 • resolution (NCEP, 2015). The CHIMERE simulation 15 domain contains 799 × 399 cells at 5 km resolution. The geometry of the domain, which has a Lambert-conformal projection, is shown in Fig. 1. The top of model is placed at 150 hPa, with either 20, 50 or 99 vertical layers to evaluate the impact of vertical resolution on the volcanic plume. Horizontal advection in the CHIMERE model has been represented using the classical Van Leer (1977) second-order slope-limited transport scheme.

20
Total concentration for all species together will be noted C (number of gas particles per unit volume; molec.m −3 ; corresponding to air density), concentration of a particular species s will be noted C s (number of molecules of species s per unit volume; molec..m −3 ), mixing ratio for species s will be noted α s = Cs C . Continuity equation for the motion of air is as follows: where Φ = Cu is the total flux of molecule number. 25 In the following equations, i, j, k are the indices for the two horizontal directions and the vertical direction respectively,  through lateral cell boundaries, positively oriented towards increasing i and j values respectively, are noted F i+ 1 2 ,j,k and F i,j,k . The mass flux components need to verify the discretized form of Eq. 1: The continuity equation for species s is as follows: or equivalently: In CHIMERE, Eq. 4 is discretized as: whereᾱ are the reconstituted values of mixing ratio α on the facet indicated by the indices (α i,j,k+ 1 2 for the top facet of cell i, j, k etc.). If the mass flux values F are such that Eq. 2 is exactly verified, then Eq. 5 ensures that the mixing ratios α s are not affected by wind-mass discrepancies: in particular, a species with an initially uniform mixing ratio will maintain it after transport is applied. This is why it is so critical for chemistry-transport models to have Eq. 2 verified exactly.
Eq. 5 raise two important questions: 1. How to express the interpolated mixing ratiosᾱ s , which is the task of the transport scheme ?
2. How to enforce exact verification of Eq. 2 to ensure the absence of mass-wind inconsistencies ? 5

Vertical wind strategy
In CHIMERE, as in most other chemistry-transport models (Emery et al., 2011), with the notable exception of chemistrytransport modules that are embedded within a meteorological model and use the same grid and time step as it is most notably the case of WRF-CHEM (Grell et al., 2005), the model does not have access to mass flux components F i± 1 2 ,j± 1 2 ,k± 1 2 and a density field C that verify Eq. 2. This is a big problem when it comes to advecting species mixing ratios since as we have seen 10 above, the fact that the CTM is able to transport species with Eq. 5 while maintaining uniformity of initially uniform mixing ratios critically depends on that property. Usually, chemistry-transport models rely on the typical outputs of meteorological models, namely the instantaneous values of winds at the meteorological model cell boundaries, and instantaneous values of density. From these variables, it is possible to evaluate the mass fluxes through the CTM cell boundaries, obtaining interpolated mass flux values that are close to the "real" mass flux values from the model. We will note these values F i± 1 2 ,i± 1 2 ,i± 1 2 , and 15 C i,j,k . With these interpolated values, and after discretization in time is also performed, Eq. 2 is not verified, and is turned into: where ε i,j,k is a spurious matter creation term due to mass-wind inconsistencies in the interpolated density and mass-flux values from the meteorological model.
As discussed above, the wind in CHIMERE is normally reconstructed from the bottom to the top of the model in order to 20 prevent mass-wind inconsistency issues.
To enforce Eq. 2, reconstructed vertical fluxes F i,j,k+ 1 2 are produced from the following constraints: Eq. 7a gives the boundary condition to vertical mass flux reconstruction (no incoming mass flux of air comes from the 25 ground surface). Eq. 7b ensures that, in each CTM cell and over one CTM time step, Eq. 2 is strictly verified (in the form of 7b) and Eq. 5 can be integrated using the interpolated horizontal fluxes F i± 1 2 ,j± 1 2 ,k and the reconstructed vertical fluxes F i,j,k± 1 2 . This traditional approach will be labelled "NODIV" (for NO DIVergence) hereinafter.
In this study, we introduce and test a new approach that permits to bypass the need for a reconstructed vertical mass flux and work directly with the interpolated vertical mass fluxes F i,j,k+ 1 2 while still maintaining conservation of uniform mixing ratios. To explain this approach, we need to expand Eq. 5 as follows: where δᾱ i,j,k+ 1 2 =ᾱ i,j,k+ 1 2 − α i,j,k , and analogous definitions for the other δᾱ terms. Injecting Eq. 6 into Eq. 8 we obtain: After simplification: is initially uniform, then all the δᾱ terms vanish, and mixing ratio uniformity will be maintained after integrating Eq. 8 if, and only if, the mass-wind inconsistency term ε i,j,k is zero. This is already well-known but with this formulation we can obtain a modified version of Eq. 5 that will enforce mixing 20 ratio preservation even if the mass-wind inconsistency term ε i,j,k is not zero: Eq. 11 will be solved in the simulation labeled WRFW, with mass fluxes directly interpolated from the meteorological model winds in the three directions. It must be noted that mass conservation is not enforced by this equation: the additional term is an artificial mass production/loss term that breaks the conservation of total mass of species s over the entire domain. If we summarize this part, the NODIV simulation (classical approach) enforces tracer mass conservation and tracer mixing ratio conservation. This is obtained to the expense of irrealistic vertical transport, since mass-wind consistency is, in this 5 approach, enforced by artificially modifying the vertical mass fluxes. However, this reconstructed wind is significantly different from WRF input data while reaching the upper troposphere ( Figure S1 in supplements), and this approach induces excessive transport across tropopause. Vertical wind distribution comparisons between WRFW and NODIV strategies ( Figure S1) show that more vertical diffusion is expected in the NODIV strategy in the upper troposphere. The WRFW approach that we propose here, on the other hand, permits mixing ratio conservation and the use of realistic vertical mass fluxes, to the expense of mass 10 conservation. While non-conservation of mass is obviously a significant drawback for a transport strategy, we will quantify this problem of non-conservation of mass, as well as the problems introduced by artificial reconstruction of vertical mass fluxes in the representation of vertical transport in the NODIV approach ( Figure 3).

Vertical advection scheme
After discretizing the advection equation for species s in the form of Eq. 5, the point of the vertical transport scheme is to 15 estimate the reconstructed mixing ratiosᾱ i,j,k+ 1 2 , for k varying between 1 and the number of vertical levels nz. The most simple way of doing so is the Godunov donor-cell scheme, simply evaluatingᾱ s,i,j,k+ 1 2 as: This order-1 scheme is mass-conservative but extremely diffusive. It is therefore important to find more accurate ways to 20 estimateᾱ s,i,j,k+ 1 2 .

The Van Leer (1977) scheme
The second-order slope-limited scheme of Van Leer (1977) brought to our notations yields the following expression ofᾱ s,k+ 1 2 (for F i,j,k+ 1 2 > 0). 25 where ν = F i,j,k+ 1 2 ρ i,j,k V i,j,k is the Courant number for the donor cell i, j, k, V i,j,k being its volume: if ν > 1, then more mass leaves the cell than the mass that was initially present and the Courant-Friedrichs-Lewy condition is violated, yielding numerical instability. If α s,k is a local extremum of mixing ratio ((α s,k − α s,k−1 ) (α s,k+1 − α s,k ) ≤ 0), no interpolation is performed and the scheme falls back to the simple Godunov donor-cell formulation (Eq. 12). This order-2 scheme has been used for decades in chemistry-transport modelling, being a good tradeoff between reasonably weak diffusion, at least compared to more simple schemes such as the Godunov donor-cell scheme, and small computational burden compared to higher-order schemes such as the Piecewise Parabolic Method (Colella and Woodward, 1984).
The Després and Lagoutière (1999) The scheme of Després and Lagoutière (1999) is defined by their equations 2 to 4. If F i,j,k+ 1 2 > 0, these equations brought to our notations, adapted to the flux form of Eq.5 and ignoring the i, j indices to shorten the expression, give: with the same notations as for the Van Leer (1977) scheme (above). As above, if ((α s,k − α s,k−1 ) (α s,k+1 − α s,k ) ≤ 0), no interpolation is performed and the scheme falls back to the simple Godunov donor-cell formulation (Eq. 12). As stated by its 10 authors, this scheme is antidiffusive. Unlike other schemes such as the Van Leer (1977) scheme described above, two unusual choices have been made by the authors in order to minimize diffusion by the advection scheme: -Their scheme is accurate only to the first order -The scheme is linearly unstable, but non-linearly stable (their Theorem 1) The idea of the authors has been to make the interpolated valueᾱ s,k+ 1 2 as close as possible to the downstream value (α s,k+1

15
if the flux is upward). This property is desirable because it is the key property in order to reduce numerical diffusion as much as mathematically possible while still maintaining the scheme stability. The authors present 1d case-studies with their scheme obtaining extremely interesting results: fields that are initially concentrated on one single cell do not occupy more than 3 cells even after a long advection time (their Fig. 2), sharp gradients are very well preserved (their Fig. 1), and, more unexpectedly due to its antidiffusive character, the scheme also performs well in maintaining the shape of concentration fields with an initially 20 smooth concentration gradient. After extensive testing, these authors also suggest (their Conjecture 1) that convergence of the simulated values towards exact values occur even if the time step is reduced before the space step: in simpler terms, this means that the scheme performs very well even at small CFL values, a property that is not shared by most advection schemes.

Model parameters tested
The various possible parameter combinations between the vertical flux (NODIV for reconstructed vertical fluxes, WRFW for 25 interpolated vertical mass flux), the vertical transport scheme (VL for Van Leer (1977), DL for Després and Lagoutière (1999)) and the number of vertical levels (20, 50, 99) are summarized in Table 1. Following all the possible combinations between these parameters, 12 simulations have been performed. The time and altitude profiles for injection of SO 2 into the atmosphere (Table 2) have been obtained using SO 2 flux measurement data from the ground-based DOAS FLAME (Differential Optical Absorption Spectroscopy FLux Automatic MEasurements) scanning network (e.g. Salerno et al., 2018). This method accurately measures SO 2 fluxes during passive degassing and effusive and explosive eruptive activity using a plume height inverted by an empirical relationship between plume height and 5 wind speed (Salerno et al., 2009). In explosive paroxysmal events, as in the case in this study, the plume is ejected to higher altitudes and this linear height-wind relationship can not be used and mass flux is retrieved in post-processing using the plume height estimated by visual camera and/or satellite retrieval.
On March 18 2012, between 06 UTC and 15 UTC, a total SO 2 emission of 2.94 kt has been reported by this method. 91 % of this mass has been released within 2 hours around 12 km of altitude. In CHIMERE, emissions have been distributed into a 10 single model cell based on altitude injection. Two alternates cases have been tested, considering -1 km and +1 km for all of the injection heights defined in Table 2.

IASI and OMI instruments
To evaluate the numerical parameters tested in our simulations, satellite based information were used to evaluate the SO 2 plume transport and vertical distribution. SO 2 column observations are provided by the Infrared Atmospheric Sounding Interferometer instrument (IASI) on board the Metop-A European satellite and the Ozone Monitoring Instrument (OMI) on board the Aura NASA satellite (Levelt et al., 2006;McCormick et al., 2013). IASI instrument is operating between 3.7 and 15.5 µm, including 5 SO 2 ν 1 (around 8.7 µm) and ν 3 (7.3 µm) bands (Carboni et al., 2012). IASI scheme (Carboni et al., 2012(Carboni et al., , 2016 provides SO 2 total column content and altitude and it is shown in Figure 2 together with OMI SO 2 columns and CHIMERE outputs. In our analysis we used IASI retrievals from March 18, 17UTC (Fig. S2), March 19, 06UTC (Fig. 2) and March 19, 17UTC (Fig. S3).
OMI data are obtained from the NASA GIOVANNI platform 2 . OMI Data are provided with a 0.25 • × 0.25 • resolution and a daily coverage, at 12AM UTC for 18 th , 19 th and 20 th of March over the studied area. All instruments soundings are resumed 10 in Table 3.

The impact of the injection altitude on plume's transport
Alternative injection height scenarios have been tested, either lifting or lowering the injection height at all times by 1 km, thereby lifting maximum injection heights up to 13 km (res. lowering it down to 11 km) instead of 12 km in the reference 15 simulation. These tests have shown that plume trajectories are strongly sensitive to this parameter (Fig. 1), which is an effect of wind shear. For each satellite sounding, the coordinates of the model column with the strongest vertically integrated SO 2 content (molecules SO2 .cm −2 ) have been selected and considered as SO 2 plume centroids. Doing so with the available 3 IASI  soundings and 3 OMI soundings gives 6 points on the satellite-retrieved Etna plume trajectory, ranging from Sicily to Western Iran. The IASI/OMI centroids and constructed plume trajectory are displayed in black on Figure 4.
Simulations and satellite plumes initial position do not correspond to Etna location because OMI first sounding is at 12 UTC on March 18, 6 hours after the beginning of the eruption. Compared to the trajectory reconstituted from OMI and IASI observations, the plume injected at 11 km seems to be transported too far towards the South, while the plume injected at 13 km 5 appears to be shifted to the North compared to observations. This observation is largely independant of the model parameters (NODIV-VL-99 and WRFW-DL-99 are shown on Fig. 1), suggesting that the configuration with a maximum injection height at 12 km is the best configuration. Therefore, only this choice for injection height will be retained for the rest of the study. In addition, sensitivity to injection vertical profile has been investigated with 3 options: -Injection to a unique altitude -Injection with a full width at half maximum of 100 m (Boichu et al., 2015) -Injection with a full width at half maximum of 300 m These sensitivity tests have shown little differences between the various cases, even in 99 vertical levels resolution, with plumes 5 close trajectories and vertical diffusion. Injection to a unique altitude -consequently in a unique cell -has been conserved and used to perform and evaluate the vertical diffusion strategies. Figure 3 shows the evolution of the total mass of tracer inside the simulation domain as a function of time. Several features from this figures need to be commented.

10
As it could be expected, the simulations with the WRFW wind strategy, due to the additional term in Eq. 11, does not enforce mass conservation. In theory, the additional term, designed to ensure mixing ratio conservation in spite of masswind inconsistencies, can result into either an artificial increase or an artificial decrease in simulated tracer mass. In the three simulations that are shown on Fig. 3, the amount of tracer present in the domain just after the end of the eruption overshoots 13 https://doi.org/10.5194/gmd-2020-62 Preprint. Discussion started: 17 April 2020 c Author(s) 2020. CC BY 4.0 License. the expected mass, by 20% in the simulation in the simulation with 20 vertical levels, 10% in the two simulations with a larger number of levels. In the later stages of plume evolution, when the plume is spread over a larger area the spurious evolutions in tracer mass become weaker, less than 5%.

Horizontal transport evaluation with OMI and IASI instruments
In this section we aim to determine how the various parameters tested (Table 1)  It can be seen that the DL vertical scheme has better agreement with the observations than the VL vertical scheme, with respectively a mean gap of 189 km and 316 km -with NODIV wind strategy option fixed. The WRFW wind strategy also shows better agreement with soundings than the NODIV strategy, with a mean gap of respectively 230 km and 316 km -with 25 VL vertical scheme fixed.
To complete the centroids-gap analysis, agreement between satellite measurements and model simulations in the zonal and meridional displacements of the centroids has been calculated, as expressed in Equation 16 for a given simulation (SIM): where combining IASI and OMI instruments information. a) Transport for all simulations, b) plume transport highlighted according to vertical scheme, c) plume transport highlighted according to vertical wind strategy. The points along the trajectories correspond to those listed in Table 3.
where the i index refers to sounding numbers (Table 3), λ OBS,i and Φ OBS,i refer to the geographic coordinates of the observed centroid for sounding i at time t i (Table 3), λ SIM,i and Φ SIM,i to the coordinates of the centroid simulated, and R is the Earth radius. In this case we focus on the displacement of the SO 2 plume between two successive observation points. The 5 intention of Eq. 16 is to build an index that not only qualifies the distance between the observed and modelled trajectories, but also the realism of the displacements followed between two successive satellite snapshots, giving penalties to simulations which would oscillate erratically around the observed trajectory and bonuses to simulations that would follow a realistic trajectory, but slightly shifted towards either side. Results are displayed on Figure 5b, and again mean value for each time series has been calculated and displayed on Figure 5b last boxes. It can be observed, as in the previous indicator case, the DL vertical scheme 10 and WRFW vertical wind strategy have brought better results than the respectively opposed VL and NODIV parameters.
It also appears clearly for both criteria that the 99 vertical levels option shows the best results to both, centroids position and trajectory comparisons to the satellite than the 50 or 20 vertical levels options. 3.4 Comparison of the simulated plume vertical structure with IASI-retrieved structure IASI observations also provide the estimated altitude of the SO 2 plume for each pixel with a valid SO 2 retrieval (Figure 2e), along with the corresponding uncertainty. For each of the three available IASI soundings (Numbers 2, 3 and 5 in Table 3), we have extracted the plume altitude and its associated uncertainty for the pixel with the highest content of SO 2 . Comparison of this plume altitude with the same calculation made on the plume simulated by CHIMERE is shown on Figure S4 (in supplements).

5
It can be observed that for the IASI first available sounding, soon after the Etna eruption, concentration maximum altitude for CHIMERE is consistent with IASI altitude and is found within IASI uncertainties (12 100 m ± 900 m): Along with the trajectories shown on Fig. 1, this is an indication that the highest emission altitude at 12 km in Tab. 2 is realistic.
In the IASI dataset, a bimodal altitude distribution is observed, indicating the coexistence of two separated sub-plumes during this eruption: one located to the East at higher altitude, another one located to the West at lower altitude (Figure 2), which 10 has also been observed in AEROIASI-sulphates soundings in (Sellitto et al., 2017;Guermazi et al., 2019). This separation is due to the sharp separation between emissions at very high altitudes ( 12 km) and emissions below 7, 5 km (Tab. 2) and to the fact that at the Etna latitude wind shear is generally strong with steady westerly winds in the higher troposphere and more variable winds in the lower troposphere. Since most of SO 2 is emitted around 12 km (Table 2), most of SO 2 mass is found in the Eastern plume. From each available IASI observation of the plume (soundings 2, 3 and 5 in Tab. 2), a 15 transition longitude that separates the western plume from the eastern plume has been identified. These longitudes are given in Tab. 2 and can be compared to Figs. S2 and S3. The same longitudes have been used to split the CHIMERE simulated plumes between a Eastern and Western plume. Figure 6 compares the altitude distribution between the CHIMERE simulations and IASI retrievals -20 vertical levels simulations have been removed, because of the coarse altitude resolution does not permit a useful representation of the maximum concentration's altitude distribution (see Figure S4, in supplements). IASI data -20 vertical levels are clearly insufficient to reproduce correctly even the main features of the plume advection in this case: no evolution of the plume altitude can be seen at such a coarse vertical resolution, and the plume seems to be strongly leaking through the top of model, which is not the case with 50 or 99 vertical levels 15 -The simulation that diffuses less the plume is the WRFW-DL-99 -Using an antidiffusive advection scheme permits to reduce diffusion almost as strongly as increasing the number of levels: for example, simulation NODIV-DL-50 preserves the maximum concentration of SO 2 in the plume and the thin plume structure as well as simulation NODIV-VL-99, with a calculation cost divided by two due to the reduced number of vertical levels. Therefore, the use of an antidiffusive advection scheme is a very attractive means of reducing numerical 5 diffusion on in the vertical direction without increasing the computational cost of the simulation.
-Using the realistic vertical wind WRFW instead of the vertically reconstructed NODIV wind also reduces numerical diffusion and avoids intermittent leaks of the tracer through the upper model boundary, as can be seen by comparison between the NODIV-VL-20 and WRFW-VL-20 simulations for example.
-Combining an antidiffusive advection scheme, the use of real vertical wind and the largest number of levels system-10 atically permit to reduce plume diffusion. Qualitatively, the impact of the antidiffusive transport scheme on reducing vertical diffusion seems to be more prononced than that of using real vertical winds instead of reconstructed ones: in

Parameters impact on SO 2 dispersion
In the aim to evaluate the SO 2 overall diffusion (vertical and horizontal) following the present Etna eruption, we have compared the minimum volume in which 50 % of SO 2 mass can be found for each time step and for each simulations. Volume evolution have been represented on Figure 8a. 20 levels simulations present the highest volume occupied, and 99 levels simulations the 25 lowest volume occupied. The simulation which restricts most the diffusion is WRFW-DL-99 one, and by opposition, the one which contains less efficiently the plume is NODIV-VL-20.
To assess the impact of the applied vertical scheme or vertical wind strategy, the volume ratio evolution between VL and DL vertical scheme for each vertical levels number has been calculated -with the wind strategy fixed (NODIV). Likewise, the volume ratio evolution between NODIV and WRFW wind strategy for each vertical levels number has been calculated - 30 with the vertical scheme fixed (VL). Volume ratio evolution have been summarized on Figure 8b. It appears that the vertical transport scheme DL reduce strongly the diffusion compared to VL, as the volume ratios V ol N ODIV −V L /V ol N ODIV −DL increase with time for the three level numbers configuration. Results are less clear for the vertical wind strategy, as the volume ratios V ol N ODIV −V L /V ol W RF W −V L slightly increase in most of the cases, with a stronger effect in the 20 vertical levels case. Finally, we observe that the combination of WRFW and DL parametrisation has consequently reduced the atmospheric diffusion such as the occupied volume for WRFW-DL-20 is quite similar to the NODIV-VL-99 case: using a realistic vertical wind field and an antidiffusive scheme is, for this criterion, as efficient as refining the model vertical resolution by a factor of 5 5. By extension, it has been observed that volcanic plume shape has been modified by DL and WRFW parameters, reducing the surface area containing 50 % of SO 2 total mass.

Conclusion
In order to analyse atmospheric diffusion, the Etna eruption of March 18, 2012 has been modelled using the CHIMERE CTM in the aim to propose and test strategies to improve representation of atmospheric vertical diffusion which has previously been 10 described in multiple studies (Colette et al., 2011;Boichu et al., 2015;Mailler et al., 2017, e.g.) as over-diffusing, inducing an excessive spread of the simulated plumes. First, the sensitivity to plume injection height and profile have been evaluated, following the plume trajectory with satellites retrievals. It appeared in these tests that the trajectory is highly sensitive to the injection altitude. The intermediate option (injection at 12 km) has been retained and tests and comparisons have been made with this injection altitude. No significant impact of plume injection profile has been observed and the most simple case of a 15 unique altitude emission has been retained.
In order to reduce the excessive spread of the plume in the vertical direction due to numerical diffusion, three possible approaches have been tested: increasing the number of levels (20, 50 and 99 level simulations have been performed), using the anti-diffusive scheme of Després and Lagoutière (1999) instead of the classical Van Leer (1977) second-order slopelimited scheme, and using realistic vertical winds instead of vertically reconstructed vertical winds, as it is usually done in 20 CTMs (Emery et al., 2011), to the expense of tracer mass conservation. Our results show that, as expected and as already shown in earlier studies, 20 vertical levels are clearly not sufficient to usefully represent any property of vertical transport and dispersion of this plume, and that increasing the number of vertical layers to 50 or to 99 brings significant added value in all respects: horizontal trajectories are improved compared to satellite measurements, vertical diffusion is reduced and maximal concentrations are preserved better. Also very effective is the use of the Després and Lagoutière (1999) anti-diffusive transport 25 scheme instead of the Van Leer (1977) scheme. To our knowledge, this scheme has never been used in chemistry-transport studies, and we show here that this strategy has a very strong potential in preventing simulations to be affected by excessive vertical diffusion without dramatically increasing the number of vertical levels. In our simulations, using the Després and Lagoutière (1999) scheme with 50 levels only has lead to performances that are comparable to the ones obtained with the Van Leer (1977) scheme and 99 levels (compare Fig. 7c to Fig. 7h or Fig. 7f to Fig. 7k). With an equivalent number of vertical 30 levels, maximum concentrations in the plume after slightly more that 48 hours of atmospheric transport are about four times as strong in a simulation with the Després and Lagoutière (1999) scheme than in the same simulation but with the Van Leer (1977) scheme. In addition, in principle, increasing vertical resolution is meaningful only in cases where plume injection altitude is known with a very good accuracy, while reducing vertical diffusion by other ways do not require increased accuracy in the a priori knowledge of injection height and profile. Finally, it has been shown that using realistic vertical winds instead of reconstructed vertical winds also improve the horizontal trajectories of the plume, when compared to satellite observations, and reduce plume diffusion in terms of minimum volume containing at least half of the plume mass. It needs to be recalled here that this strategy does not guarantee mass tracer conservation. Even though its impact in this respect has been shown to 5 be quite minor in our study, except in the simulations with 20 vertical levels where it generated an initial excess in tracer mass of about 15%, this characteristic needs to be kept in mind and accounted for by potential users of this strategy.
These different strategies need to be further studied in different cases to determine whether they can be generalized in CTMs in order to reduce vertical diffusion issues for all pollutants and all kinds of problems, if they are useful only for long-distance transport of inert plumes as we simulated here. For example, how does the Després and Lagoutière (1999) perform in preserving 10 smooth gradients, transitions between the PBL and the free troposphere, or gradients in ozone concentrations at the tropopause ? Regarding the use of realistic vertical wind fluxes interpolated from meteorological outputs, does it help reduce excessive ozone transport through the tropopause as identified by Emery et al. (2011) ? Can this method be generalized to the general chemistry-transport modelling of the troposphere in spite of the mass conservation issues that are intrinsic to this method ? This study is a call to reopen the issue of limiting vertical mass diffusion in eulerian CTMs: complementary to Zhuang et al. (2018) 15 who emphasized on the need for sufficient vertical resolution, which is confirmed by the present study, we propose two new approaches solve this long-standing problem, including anti-diffusive transport schemes and better representation of vertical mass fluxes throughout the troposphere.