2.1 The MPTRAC Lagrangian transport model

For the simulation of the dispersion and depletion of SO₂ from the Ambae eruption, we applied the Massive-Parallel Trajectory Calculations (MPTRAC) model (Hoffmann et al., 2016, 2022a). Mainly, trajectories of air parcels are calculated by given horizontal winds and vertical velocities of meteorological input data, with additional stochastic perturbations being added to simulate diffusion and subgrid-scale wind fluctuations. Additionally, several improved or newly developed chemical and physical modules of MPTRAC are applied in the simulations. In particular, the SO₂ mass of the air parcels is decomposed by OH oxidation and wet removal processes instead of simply using an exponential decay with a fixed lifetime as applied in our earlier studies. The new and revised modules are described in the following sections.

The transport simulations with MPTRAC are driven by the European Centre for Medium-Range Weather Forecasts’ (ECMWF’s) fifth-generation reanalysis (ERA5) (Hersbach et al., 2020) with hourly meteorological data at ${\mathrm{0.3}}^{\circ }\times {\mathrm{0.3}}^{\circ }$ horizontal resolution on 137 vertical levels. This is a significant improvement in spatiotemporal resolution compared with the previous generation ERA-Interim reanalysis (Dee et al., 2011), providing only ${\mathrm{0.75}}^{\circ }\times {\mathrm{0.75}}^{\circ }$ horizontal resolution on 60 vertical levels at 6 h time intervals. The higher-resolution ERA5 data are expected to lead to more accurate Lagrangian transport simulations compared to ERA-Interim (Hoffmann et al., 2019).

MPTRAC is a Lagrangian transport model developed with a hybrid MPI–OpenMP–OpenACC parallelization scheme for application on CPU/GPU heterogeneous supercomputers, aiming for good parallel performance and scaling efficiency. The computing tasks in MPTRAC are distributed over the compute nodes and compute cores by MPI and OpenMP and offloaded to GPUs by means of OpenACC for faster and more-energy-efficient computational performance (Liu et al., 2020; Hoffmann et al., 2022a). In this study, each individual simulation is conducted in parallel with 48 OpenMP threads. MPI parallel multi-processing for different model parameter settings is applied in the sensitivity tests, which significantly improves the overall runtime of the simulations. On the state-of-the-art high-performance computing system Jülich Wizard for European Leadership Science (JUWELS) at the Jülich Supercomputing Centre (Jülich Supercomputing Centre, 2019), individual Lagrangian transport simulations with 10⁶ air parcels over a time period of 15 d require about 1.5 h of total runtime on a compute node. GPU acceleration was not considered here as there is no increase in computational speed at this problem size.

2.1.1 Hydroxyl radical oxidation in the gas phase

The hydroxyl radical (OH) is an important oxidant in the atmosphere, which causes rapid decay of many gas-phase species. In MPTRAC, the mass decay of a trace gas of an air parcel is described by an exponential formula,

$$\begin{array}{}\text{(1)}& m(t+\mathrm{\Delta }t)=m\left(t\right)\exp \left(-k\mathrm{\Delta }t\right).\end{array}$$

Assuming that the concentration of OH is in near steady state, the pseudo-first-order reaction coefficient k is calculated as

$$\begin{array}{}\text{(2)}& k=k_f\times \left[\mathrm{OH}\right],\end{array}$$

with an effective second-order coefficient k_f and a prescribed monthly zonal mean OH field. The oxidation of SO₂ with OH is a termolecular reaction,

$$\begin{array}{}\text{(3)}& {\mathrm{SO}}_{\mathrm{2}}+\mathrm{OH}\to [{\mathrm{HOSO}}_{\mathrm{2}}{]}^{*}\stackrel{M}{\to }{\mathrm{HOSO}}_{\mathrm{2}},\end{array}$$

where the excited intermediate [HOSO₂]^* requires an inert molecule M (e.g., N₂ or O₂) to remove the energy and stabilize it into sulfate. In the high-pressure limit, the rate-limiting step is the production of [HOSO₂]^*, while in the low-pressure limit, the reaction rate depends on the abundance of M and the production of HOSO₂. Thus, the effective second-order rate coefficient of the SO₂–OH oxidation process is temperature- and pressure-dependent, which is described here by using a formula given by the NASA Jet Propulsion Laboratory (JPL) data evaluation (Burkholder et al., 2019) as

$$\begin{array}{}\text{(4)}& k_f(T,[M\left]\right)={\displaystyle \frac{k_{\mathrm{0}}\left(T\right)\left[M\right]}{\mathrm{1}+\frac{k_{\mathrm{0}}\left(T\right)\left[M\right]}{k_{\mathrm{\infty }}\left(T\right)}}}{\mathrm{0.6}}^{{\left\{\mathrm{1}+{\left[{\log }_{\mathrm{10}}\left(\frac{k_{\mathrm{0}}\left(T\right)\left[M\right]}{k_{\mathrm{\infty }}\left(T\right)}\right)\right]}^{\mathrm{2}}\right\}}^{-\mathrm{1}}}.\end{array}$$

The high-pressure limit rate $k_{\mathrm{\infty }}=\mathrm{2.9}\times {\mathrm{10}}^{-\mathrm{31}}\times (\frac{T}{\mathrm{298}}{)}^{-\mathrm{4.1}}{\mathrm{cm}}^{\mathrm{3}}{\mathrm{molecule}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$ and the low-pressure limit rate $k_{\mathrm{0}}=\mathrm{1.7}\times {\mathrm{10}}^{-\mathrm{12}}\times (\frac{T}{\mathrm{298}}{)}^{\mathrm{0.2}}$ ${\mathrm{cm}}^{\mathrm{6}}{\mathrm{molecule}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ were also obtained from the JPL evaluation.

As for the prescribed zonal mean OH fields, monthly mean climatology of OH calculated from simulation of the Chemical Lagrangian Model of the Stratosphere (CLaMS) is used by default (Pommrich et al., 2014). Another OH dataset considered in this study was obtained by integrating Copernicus Atmosphere Monitoring Service (CAMS) reanalysis OH data into monthly latitude- and pressure-dependent fields, keeping consistency with the CLaMS climatology (Inness et al., 2019). The CLaMS data have 18 latitude bins and 34 pressure levels, while CAMS data have 241 latitude bins and 25 pressure levels. CLaMS has a complete chemistry scheme on the stratospheric chemistry, while the CAMS reanalysis uses a chemical mechanism most suitable for the troposphere. A comparison of the CLaMS and CAMS datasets with in situ measurements is presented in the Appendix. Whether the difference between the two OH datasets affects the MPTRAC simulations for the Ambae case study is discussed.

The formation of OH is driven by the photolysis of ozone in the troposphere and H₂O in the stratosphere, which causes a strong correlation with diurnal variability (Minschwaner et al., 2011). In contrast to previous versions of MPTRAC, which did not take diurnal variability into account, the mean OH climatology concentration [OH]₀ is multiplied here by a scaling factor depending on the solar zenith angle θ_SZA as proposed by Minschwaner et al. (2011) to model the diurnal variations,

$$\begin{array}{}\text{(5)}& f\left({\mathit{\theta }}_{\mathrm{SZA}}\right)=\exp [-\mathit{\beta }\sec ({\mathit{\theta }}_{\mathrm{SZA}}\left)\right].\end{array}$$

The term sec (θ_SZA) is the approximate air mass factor which represents the ratio of the optical slant path to the effective vertical path. The parameter β represents the vertical optical depth. Based on Minschwaner et al. (2011), a β value of 0.6 is used for simulations covering the UT/LS region. To maintain the same mean values of the scaled data as the monthly OH field, they are divided by a normalization factor obtained by integrating the correlation factor f over longitude λ. The final OH concentration is calculated as

$$\begin{array}{}\text{(6)}& \left[\mathrm{OH}\right]={\displaystyle \frac{[\mathrm{OH}{]}_{\mathrm{0}}f\left({\mathit{\theta }}_{\mathrm{SZA}}\right)}{{\int }_{-\mathrm{180}}^{\mathrm{180}}f\left({\mathit{\theta }}_{\mathrm{SZA}}\right(\mathit{\lambda }\left)\right)\mathrm{d}\mathit{\lambda }/\mathrm{360}}}.\end{array}$$

Figure 1 shows average vertical profiles of the different OH datasets at tropical latitudes (from 23.5^∘ S to 23.5^∘ N). The CLaMS OH data (Pommrich et al., 2014), calculated using methane and ozone data taken from the HALOE climatology for the year 2005 (Grooß and Russell III, 2005), are compared with the CAMS OH reanalysis data for the years 2005 and 2018, respectively. The differences in the CAMS mean profiles between 2005 and 2018 are found to be negligible. Comparing the CLaMS and CAMS OH data, these two datasets show similar concentrations in the troposphere, while in the stratosphere, CLaMS OH concentrations are much larger than CAMS OH concentrations. To further evaluate the OH fields of the climatologies, we compared them to several NASA in situ OH measurements in the troposphere at similar altitudes, solar zenith angles, and time periods and also to Microwave Limb Sounder (MLS) satellite data in the stratosphere. Both CLaMS and CAMS data showed good agreement with the in situ measurements in the troposphere. At altitudes above 20 km, CLaMS data show much better agreement with MLS observations than the CAMS reanalysis. Further details of the comparison are presented in the Appendix.

Figure 1Vertical profiles of (a) OH concentrations over the tropics retrieved from CLaMS and CAMS climatology data, respectively, on 26 July 2018, 00:00 UTC, and (b) OH oxidation rate coefficients at the same time.

Download

2.1.3 Wet deposition

Irreversible in-cloud and below-cloud wet removal processes of trace gases need to be treated separately. In MPTRAC, the cloud liquid water content (CLWC) and the cloud ice water content (CIWC) of the meteorological input data are used to determine whether an air parcel is located within or below a cloud. The total column cloud water c_l and the cloud depth Z are additional input variables to the wet deposition module. Similar to Pisso et al. (2019), the precipitation rate P (in units of mm h⁻¹) is estimated from c_l (in units of kg m⁻²) by means of a regression analysis using the ERA5 meteorological data for a region covering the Ambae eruption (10^∘ S to 30^∘ N, 160^∘ W to 140^∘ E) over the time period from 25 to 31 July 2018, as shown in Fig. 2,

$$\begin{array}{}\text{(12)}& c_l=\mathrm{0.763}\times P^{\mathrm{0.478}}.\end{array}$$

Figure 2Scatter density plot of total column cloud water versus precipitation rate of ERA5 meteorological data for a region and time period covering the Ambae eruption (see text for details). Color coding indicates the normalized density of the data points. The black line represents the regression of the data as given in Eq. (12). The red line shows the expression given by Pisso et al. (2019).

Download

For the modeling of in-cloud wet deposition, a similar exponential removal as in Eq. (1) is used with a scavenging coefficient Λ (in units of s⁻¹). The in-cloud wet deposition process is implemented with two schemes for a choice. The first scheme is based on the rain-out rate of cloud water, which determines the removal of a soluble trace gas taken up by cloud droplets and removal by precipitation according to the solubility of the trace gas. The in-cloud scavenging coefficient is calculated by multiplying the partition ratio of the species in the aqueous phase versus the gas phase α with the cloud water removal rate (Slinn, 1974; Levine and Schwartz, 1982; Garrett et al., 2006),

$$\begin{array}{}\text{(13)}& \mathrm{\Lambda }=\mathit{\alpha }{\displaystyle \frac{P}{LZ}},\end{array}$$

where L is the volume liquid water content of the cloud (in units of m³ m⁻³), Z is the depth of the cloud layer determined by the pressure of cloud top and cloud bottom as taken from the meteorological input data, and α is defined by

$$\begin{array}{}\text{(14)}& \mathit{\alpha }={\displaystyle \frac{N_l}{N_g}}={\displaystyle \frac{H{P}_{x}L}{P_x/RT}},\end{array}$$

according to the ideal gas law, where H is Henry’s law coefficient, R is the universal gas constant, and P_x is the partial pressure of species x.

Combining Eqs. (13) and (14), a formula similar to the approach of the HYSPLIT model (Draxler and Hess, 1998) is obtained,

$$\begin{array}{}\text{(15)}& \mathrm{\Lambda }=\mathit{\eta }HRTP{Z}^{-\mathrm{1}}.\end{array}$$

The factor η represents a temperature-dependent retention coefficient in cloud versus equilibrium concentration in liquid water. The default value of η is set following Webster and Thomson (2014), assuming the cloud between 238.15 and 273.15 K to be in the mixed phase and a retention ratio in ice clouds to be 0.15:

$$\begin{array}{}\text{(16)}& \mathit{\eta }=\left\{\begin{array}{ll}\mathrm{1},& \text{ if }T\ge \mathrm{273.15}\mathrm{K}\\ \mathrm{0.15},& \text{ if }T\le \mathrm{238.15}\mathrm{K}\\ \mathrm{0.15}+\frac{T-\mathrm{238.15}\mathrm{K}}{\mathrm{273.15}\mathrm{K}-\mathrm{238.15}\mathrm{K}}& \text{ if }\mathrm{238.15}\mathrm{K}<T\\ (\mathrm{1}-\mathrm{0.15}),& <\mathrm{273.15}\mathrm{K}.\end{array}\right.\end{array}$$

SO₂ is a moderate soluble gas with a Henry's law constant of 1.3 M atm⁻¹ at 298 K (Sander, 2015). However, the solubility of SO₂ strongly depends on the pH value because it undergoes dissociation. In the simulations for SO₂, the partition of SO₂ in clouds is represented by an effective Henry's law constant related to the pH value to account for the dissolution and dissociation of SO₂ in cloud water (Berglen et al., 2004),

$$\begin{array}{}\text{(17)}& H_{\mathrm{eff}}\left({\mathrm{SO}}_{\mathrm{2}}\right)=H\left({\mathrm{SO}}_{\mathrm{2}}\right)\times \left[\mathrm{1}+K^{\prime }/\left[H^{+}\right]+K^{\prime }{K}^{\prime \prime }/{\left[{\mathrm{H}}^{+}\right]}^{\mathrm{2}}\right].\end{array}$$

Here K^′ and K^′′ are the first and second dissociation constants of SO₂, using the formulation of Berglen et al. (2004), and H(SO₂) is Henry's law constant (Sander, 2015). With this method, the pH value becomes a control parameter of the model affecting the wet deposition, for which we assumed a default value of 4.5, guided by earlier work (Berglen et al., 2004; Koch et al., 1999).

For comparison, we also implemented another in-cloud wet deposition scheme applying an empirical exponential expression given by

$$\begin{array}{}\text{(18)}& \mathrm{\Lambda }=\mathit{\eta }a{P}^b.\end{array}$$

The choice of the parameters a and b follows the NAME model (Webster and Thomson, 2014).

The below-cloud wet deposition is a washout process through impact or diffusion with raindrops. With respect to SO₂, the washout rate has a typical magnitude of $\sim {\mathrm{10}}^{-\mathrm{5}}$ s⁻¹ (Maul, 1978; Martin, 1984; Elperin et al., 2015). Here, we set the parameters for SO₂ to $a=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}$ and b=0.616, which follows the settings in FLEXPART (Pisso et al., 2019).

2.2 Satellite data products

2.2.1 AIRS sulfur dioxide measurements

The Atmospheric Infrared Sounder (AIRS) (Aumann et al., 2003) is an infrared spectrometer aboard the National Aeronautics and Space Administration’s (NASA's) Aqua satellite. Aqua operates in a sun-synchronous low Earth orbit at an orbit altitude of 705 km, providing nearly continuous measurements since September 2002. The AIRS instrument has across-track scanning capabilities. The swath width is 1780 km, consisting of 90 footprints per scan with a footprint size of 13.5 km ×13.5 km at nadir. The measurements take place at about 01:30 and 13:30 local time for the descending and ascending sections of the orbits, respectively.

To detect the presence of volcanic SO₂ using the AIRS radiance measurements, an SO₂ index defined as the brightness temperature difference between 1407.2 and 1371.5 cm⁻¹ in the 7.3 µm SO₂ waveband is used (Hoffmann et al., 2014). The SO₂ index of Hoffmann et al. (2014) is most sensitive in the column density range of about 10 to 200 DU (Dobson units) at altitudes of 8 to 13 km, which covers explosive volcanic eruptions with SO₂ injections into the UT/LS region. For SO₂ index values larger than 4 K, the SO₂ index is able to clearly detect volcanic plumes.

Figure 3Observed (TROPOMI) and simulated (MPTRAC) SO₂ total mass curve (a) and mass loss (b) of the July 2018 Ambae eruption as a function of time. The blue shading represents the total error of the TROPOMI observations, combining random and systematic error components.

Download

3.1 Baseline simulation

In this section, we present a baseline simulation of the dispersion and depletion of the volcanic SO₂ plume of the Ambae eruption in July 2018. In our initial tests, it was found that, with only OH oxidation and wet deposition being considered, the simulations cannot fully explain the observed fast depletion of SO₂. Therefore, we newly implemented the process of in-cloud oxidation with H₂O₂. As the plume transport in this case occurs in a high-CAPE tropical region, the convection module was also included. A CAPE threshold of 1000 J kg⁻¹ was used to include moderate to strong subgrid-scale convection in the baseline simulation. All the parameter settings in the baseline simulation use the default values as introduced in Sect. 2. In the following sections, we will analyze the sensitivity of parameter choices for the above-mentioned chemistry and physics modules of MPTRAC with respect to the baseline simulation setup.

Figure 4Evolution of the Ambae SO₂ plume in AIRS SO₂ observations as well as column densities from TROPOMI observations and MPTRAC simulations. The time of the satellite observations shown here was restricted to ±3 h around 00:00 UTC. The black triangle shows the location of Ambae Island. MPTRAC simulation results have been sampled on the TROPOMI footprints.

Download

Figure 3a shows the SO₂ total mass curve calculated from the TROPOMI data. The total mass of SO₂ in the UT/LS region shows a strong increase on 26 July 2018, reaching a peak on 28 July, and decreases back to pre-eruption levels on 7 August 2018. According to different satellite instrument observations of the Support to Aviation Control Service (SACS, 2022), the volcanic plume is detected directly over the volcano from 26 to 27 July. The bulletin report of the Global Volcanism Program (GVP) (Krippner and Venzke, 2019) states that the main explosion occurred on 26 July and another two intense episodes, producing volcanic lightning, occurred on 27 July. Based on the different satellite measurements, observational reports as well as empirical testing, we initialized the MPTRAC baseline simulation for the Ambae case study by releasing 10⁶ air parcels with a total mass of 450 kt starting on 26 July, 00:00 UTC, over a time period of 36 h, assuming a Gaussian vertical profile with the maximum centered at 14 km and a full width at half maximum (FWHM) of 4 km. Although there are more advanced methods available to accurately estimate the timing and vertical distribution of volcanic emissions (Hoffmann et al., 2016; Wu et al., 2017; Cai et al., 2022), we did not apply these techniques in this study because these techniques do not fully account for SO₂ chemistry. This study focuses on testing the new and revised chemical and physical modules and on conducting sensitivity tests with respect to the parameter choices. As discussed below, the baseline simulation with constant emission rate and Gaussian vertical profile represents the Ambae case reasonably well so that meaningful testing and evaluation can be conducted.

Figure 5SO₂ mass loss relative to the total mass (in percent) due to OH chemistry (a), wet deposition (b), and H₂O₂ chemistry (c) as well as total mass curves (d) derived from simulations with different release heights centered at 11, 13, 15, 17, and 19 km.

Download

To properly assess the total mass evolution and the SO₂ plume patterns of the baseline simulation, the model output was sampled at the exact time and location as the TROPOMI satellite footprints by calculating the column density of all parcels located within a horizontal search radius of 7 km, as introduced in Cai et al. (2022). Air parcels in data gaps of the satellite observations will be excluded in the evaluation of the simulation results. The model results are multiplied by an altitude-dependent sensitivity profile derived from TROPOMI averaging kernel data. The TROPOMI averaging kernel indicates full sensitivity above 10 km of altitude, and the sensitivity only significantly decreases in the lower troposphere Cai et al. (2022). Since in this case most of the SO₂ injections and the plume height were located above 10 km, the inclusion of the averaging kernel does not significantly change the model results. A lower threshold of 1 DU was applied to both the simulation data and the TROPOMI observations to eliminate the effects of noise. With this approach, the sample output of the model can be quantitatively analyzed and compared with the TROPOMI satellite observations. As shown in Fig. 3a, the MPTRAC baseline simulation yields good agreement with the total mass curve derived from the TROPOMI data. The overall SO₂ lifetime from the simulation is quite similar to the observational data.

Figure 3b shows the relative mass loss relative to the total emissions over time with regard to each module. At the end of the simulation, OH chemistry, H₂O₂ chemistry, and wet deposition lead to a mass loss of 45 %, 14 %, and 21 % of the total mass burden, respectively. The remaining mass loss in Fig. 3a is due to air parcels leaving the study region or having SO₂ mass below the filtering threshold so that they are not accounted for in the mass budget anymore. The loss rates due to in-cloud removal processes, including wet deposition and the aqueous-phase H₂O₂ oxidation, strongly depend on the location of the volcanic plume relative to the cloud fields. The loss due to OH oxidation is a step-shaped curve due to the diurnal variations of the OH concentration, leading to a faster loss rate at daytime and nearly zero loss at nighttime. Simulations using an OH field with diurnal variations will vary by 0.7 % compared to using an OH field without diurnal variations. Differences are relatively low in this case because the observations occur at about 13:00 local time. Overall, the diurnal variations of the OH field have little effect on the SO₂ decay rates over the entire simulation period, because the observational timescale is over several days and the diurnal variations are averaging out.

Figure 6Background images represent the distributions of ERA5 cloud top height on 31 July 2018, 00:00 UTC. Dots represent air parcels released at altitudes of (a) 11 km, (b) 13 km, (c) 15 km, and (d) 17 km. A yellow dot indicates that the air parcel is below the cloud top height, triggering the wet deposition module, while a blue dot represents the opposite.

Download

Figure 7Same as Fig. 6, but background images represent the distribution of ERA5 CAPE values on 31 July 2018, 00:00 UTC. Yellow dots indicate air parcels below the equilibrium level and local CAPE values larger than 1000 J kg⁻¹, triggering the convection module, while blue dots represent the opposite.

Download

Figure 4 shows comparisons of horizontal maps of the SO₂ plume from the baseline simulation with AIRS and TROPOMI measurements, respectively. The bulk of the Ambae plume moved eastwards, and the moving speed is faster at lower altitudes. A small plume at heights of 10 to 12 km moved northwards after 28 July at 160 to 180^∘ W, encountering strong wet deposition. The model results match the satellite observations qualitatively well with similar location and shape of the plume. At the end of the simulation the part of the plume below the cloud top is almost depleted, and the remaining part of the plume above the cloud top somewhat shows deviations from the observations due to the unpolished temporal and vertical variations of the emission estimates. Nevertheless, we consider the baseline simulation in its present form to be suitable for further evaluation and sensitivity tests.

Figure 8SO₂ mass loss relative to total mass (in percent) due to the OH chemistry module (a) and total mass curves (b) for simulations with different OH datasets, including CLaMS data in 2005, CAMS data in 2005 and 2018, and CAMS regional data in 2018 (see text for details).

Download

Figure 9SO₂ mass loss relative to total mass (in percent) due the wet deposition module (a) and total mass curves (b), derived from simulations with different pH values in the wet deposition module and comparing with the wet deposition scheme of the NAME model.

Download

3.2 Sensitivity test on SO₂ release height

In this section, a series of simulations with air parcels released at different altitudes are presented to show the sensitivity of the plume injection height on the evolution of SO₂ total mass burden. The tests use different release heights centered at 11, 13, 15, 17, and 19 km with 2 km wide uniform vertical distribution. The convection module was not activated in this test to avoid vertical mixing of air parcels down to the surface during the simulations. As shown in Fig. 5, the SO₂ decay rate has an obvious dependency on release height. Releasing the volcanic SO₂ at different heights leads to rather different horizontal spread and different mass loss due to wet deposition, depending on whether the SO₂ plume encounters cloud regions or not.

As shown in Fig. 5b, the mass loss due to wet deposition shows particularly high sensitivity to the release height. Wet deposition occurs mainly on 27 to 28 July when releasing air parcels at altitudes below 14 km, while most of the wet deposition mass loss occurs around 1 August when releasing above 14 km. Almost no wet deposition occurs above 18 km because it is above the tropopause and the maximum cloud top height. Aqueous-phase H₂O₂ chemistry is most active below 14 km and has a strong sensitivity to release height, while above 14 km the H₂O₂ chemistry is very weak. The OH gas-phase oxidation does not show a clear correlation with release height, with contributions to mass loss ranging from 30 % to 60 % of the total mass during the simulation. However, note that the OH concentrations and temperature- and pressure-dependent rate coefficients of the OH oxidation are not linearly increasing with altitude (compare Fig. 1) over heights from 10 to 20 km. Above 15 km, the OH oxidation rate decreases with height, while below 15 km, the decay rate is controlled mainly by in-cloud removal and the OH oxidation rate quickly decreases.

Figure 10SO₂ mass loss relative to the total mass (in percent) due the wet deposition module (a) and total mass burden evolution (b), derived from different simulations with different retention ratio settings in the wet deposition module.

Download

Figure 11SO₂ mass loss relative to total mass (in percent) due the H₂O₂ chemistry module (a) and total mass burden evolution (b), derived from simulations with different cloud cover settings in the H₂O₂ chemistry module.

Download

Figures 6 and 7 show the cloud top height and CAPE distributions derived from the ERA5 meteorological data as well as the particles released at different altitudes, clearly indicating which part of the plume is affected by wet removal and subgrid convection, as distinguished by yellow and blue color. The plume released at the location of Ambae Island is transported to a region with high CAPE values and abundant clouds. The wet removal and subgrid-scale convection mainly take effect at altitudes below 15 km. The maximum CAPE values range approximately from 1000 to 1600 J kg⁻¹. In simulations with higher injection altitudes, the air parcels are transported mainly above the cloud top and are barely influenced by wet deposition and convection.

3.3 Sensitivity tests on OH chemistry

The SO₂ decay rate in the OH chemistry module is calculated by Eq. (2), which considers the temperature- and pressure-dependent reaction rate k_f and the OH concentration. The OH concentration is obtained from predefined monthly mean zonal mean data and the solar zenith angle correction to account for day- and nighttime conditions. The OH data used by default in MPTRAC are the monthly mean zonal mean climatology from Pommrich et al. (2014), which was calculated using the CLaMS model chemistry scheme. For another alternative, the CAMS global reanalysis provides 3 h OH data with a resolution of $\mathrm{0.75}{}^{\circ }\times \mathrm{0.75}{}^{\circ }$. To reduce memory needs and maintain consistency with the CLaMS data approach, the CAMS data are also converted into monthly mean zonal mean data. To verify the impact of inter-annual differences, we compared the results between the simulations with CAMS reanalysis data for the years 2005 (matching the CLaMS data) and 2018 (matching the Ambae eruption), respectively. For another test, simulations with regional CAMS OH data extracted from the longitude range from 160^∘ E to 140^∘ W are compared with global data.

As shown in Fig. 8, the inter-annual differences in the CAMS OH data are negligible (∼1 percentage point, pp). The differences between simulations with CAMS and CLaMS global OH data are ∼3 pp. The simulation with regional CAMS OH data for the location of the eruption has ∼4 pp difference compared with the global CAMS data. Overall, the inter-annual and regional differences of the OH datasets have a small influence on the amount of OH oxidation (less than 5 pp in this case), which suggests that using a monthly mean zonal mean climatology is a reasonable approach for modeling the OH loss of SO₂ in the UT/LS region.

Figure 12Vertical distribution of total SO₂ mass of the air parcels on 1 August 2018, 00:00 UTC, from simulations with different CAPE thresholds. The red dashed line indicates the ERA5 thermal tropopause height over the Ambae region (Hoffmann and Spang, 2022).

Download

3.4 Sensitivity tests on wet deposition

In the Ambae case study, the wet removal in clouds plays an important role in depleting SO₂ since the main plume passes through a cloud region. As introduced in Sect. 2.1.3, in-cloud wet deposition is handled as a rain-out process, which depends on the partition ratio of SO₂ in the cloud droplets and the precipitation rate. The partition ratio of SO₂ in liquid phase is defined with an effective Henry's law constant, depending on the assumed pH value. Here, we conducted a series of simulations to test the sensitivity to pH values in the range of 3 to 5, along with a comparison with the NAME wet deposition scheme (Webster and Thomson, 2014). As expected, higher cloud pH values lead to stronger wet deposition. The simulation with wet deposition according to the NAME scheme shows a similar deposition to the baseline simulation with the equilibrium scheme at pH 4 to 4.5. The difference between the two schemes is further reflected in the sensitivity to convection, which will be discussed in Sect. 3.6

The retention ratio represents the solute species in the ice phase versus that in the liquid hydrometeor. The retention coefficient may vary depending on temperature, pH, ventilation rate, accretion rate, and impact velocity (Stuart and Jacobson, 2006). Laboratory studies on the retention ratio of SO₂ suggest a range from 1 % to 60 % (Iribarne et al., 1983; Lamb and Blumenstein, 1987; Iribarne et al., 1990). As introduced in Sect. 2.1.3, a retention coefficient dependent on temperature is applied to model the differences of wet deposition in ice and liquid cloud. Figure 10 shows the SO₂ mass loss curves of simulations given different retention ratios in the range of 0.01 to 1. In this case, most of the wet deposition occurs in ice clouds, leading to a strong dependency on the retention ratio.

3.5 Sensitivity tests on H₂O₂ chemistry

Since the meteorological data cannot fully resolve the convection and detailed cloud structures, a parameter cloud cover was defined in the H₂O₂ chemistry module to represent the probability of an air parcel to be located in a cloud. When a random number in the range of 0 to 1 is larger than a given threshold for the cloud cover, the H₂O₂ decomposition will be skipped. Figure 11 shows the dependency of aqueous-phase H₂O₂ oxidation on cloud cover. Compared to full cloud cover, the amount of H₂O₂ oxidation in a case with only 10 % cloud cover is decreased by 8 pp. However, the SO₂ total mass burden is only slightly changed.

3.6 Impact of the extreme convection parametrization

Due to better spatiotemporal resolution, the ERA5 meteorological data provide more accurate information than ERA-Interim to resolve mesoscale features (Hoffmann et al., 2019). However, convective up- and downdrafts are still underrepresented in the ERA5 data. In MPTRAC, the extreme convection parametrization is applied to represent the effects of unresolved convection in the meteorological input data. In convective columns, the convection module will randomly redistribute the air parcels between the surface and the equilibrium level if CAPE exceeds a given threshold CAPE₀. As shown in Fig. 7, the SO₂ plume released by the Ambae volcanic eruption was transported to a high CAPE region, which makes it a valuable case to study the potential effects of unresolved convection. In this section, the effects of parametrized convection for different thresholds CAPE₀ will be discussed.

Figure 12 shows the vertical distributions of the SO₂ mass of the Ambae plume for different thresholds CAPE₀ from 0 to 2000 J kg⁻¹. The case of CAPE₀ being zero represents the extreme case in which convection occurs everywhere below the equilibrium level where CAPE is present. For CAPE₀ larger than 1000 J kg⁻¹, the parametrized convection is restricted to represent moderate to strong convective events only. There is almost no additional convection for CAPE₀ larger than 2000 J kg⁻¹ compared to the case without parametrized convection. When applying a larger CAPE₀, more air parcels will be transported to altitudes below the equilibrium level (∼14 km).

Figure 13 shows the impact of parametrized convection on mass loss in the SO₂ transport simulations for the Ambae case. The largest sensitivity is found for the H₂O₂ oxidation, while OH oxidation and wet deposition show little difference. This is because some air parcels are transported to heights below 14 km where the H₂O₂ oxidation is most active. The H₂O₂ chemistry is much more sensitive to height than the OH chemistry and wet deposition at altitudes where convection occurs (see Sect. 3.2). Regarding the mass evolution, the effect of the convection module enhances the depletion of SO₂ below the cloud top. At the end of simulation, most of the remaining SO₂ is above cloud top where the convection module takes no effect.

Figure 13Same as Fig. 5 but derived from simulations with different CAPE thresholds in the convection module.

Download