Transport of volcanic aerosol from the Raikoke eruption in 2019 through the Northern Hemisphere

Volcanic injections into the upper troposphere–lower stratosphere (UTLS) affect climate by altering Earth's radiation budget and atmospheric chemistry. However, the pathways by which mid-latitude eruptions spread globally remain poorly understood. We combine nighttime Compact Optical Backscatter Aerosol Detector (COBALD) profiles over Lhasa (China) with ERA5-driven Chemical Lagrangian Model of the Stratosphere (CLaMS) backward trajectories and global three-dimensional sulfur dioxide (SO₂)-based tracer simulations. With this integrated framework, we track the plume of the Raikoke eruption (21–22 June 2019) during its transport through the Northern Hemisphere and its interaction with the mature Asian Summer Monsoon Anticyclone (ASMA). Balloon-borne measurements capture the plume's arrival, vertical spreading, and dilution by air in the ASMA interior. Trajectories reveal two principal pathways from distinct Raikoke plumes: (i) an upper-level branch within the summertime stratospheric easterly flow (∼ 460–490 K) carrying the trailing filament of the vorticized volcanic plume (VVP), and (ii) a lower-level branch within the subtropical westerly jet (∼ 390–430 K) carrying the main plume. Although the ASMA can act as a transport barrier at certain potential-temperature levels, it admits mixing into the anticyclone along jet-aligned filaments and redistributes aerosols internally. We show that the transport of the volcanic plume is sensitive to parameterized mixing by varying the mixing strength in CLaMS. Portable Optical Particle Spectrometer (POPS) profiles over Boulder (USA) confirm the plume's timing and altitude, providing an independent evaluation of its transport outside the ASMA region.

1 Introduction

The impact of volcanic eruptions on climate has been a subject of widespread concern (McCormick et al., 1995; Thompson and Solomon, 2009; Solomon et al., 2011; Bourassa et al., 2012; Masson-Delmotte et al., 2021; Schmidt and Black, 2022). Large eruptions can inject significant amounts of ash and sulfur dioxide (SO₂) into the upper troposphere–lower stratosphere (UTLS). While ash typically settles out within weeks, SO₂ is oxidized into sulfate aerosols. These aerosols can persist in the stratosphere for months to years and be transported globally by atmospheric circulation, leading to significant climate impacts. Sulfate aerosols efficiently scatter incoming solar radiation back to space, producing a global net cooling in the lower troposphere (Minnis et al., 1993; Tabazadeh et al., 2002; Thompson and Solomon, 2009; Solomon et al., 2011). For instance, the 1991 eruption of Mount Pinatubo caused a global surface temperature decrease of approximately 0.5 °C over the subsequent 2 years (McCormick et al., 1995). Reconstructions based on satellite observations for 1979–2015 indicate that clusters of moderate eruptions since 2004 have imposed a persistent negative radiative forcing of about −0.08 W m⁻² (Schmidt et al., 2018). For 2014–2022, stratospheric injections from volcanic eruptions and wildfires produced a global-mean effective radiative forcing of about −0.18 W m⁻² (Yu et al., 2023). In addition, enhanced aerosol surface-area density promotes heterogeneous chemistry on polar stratospheric clouds, thereby accelerating ozone depletion (Hofmann and Solomon, 1989; Portmann et al., 1996; Solomon, 1999; Zuev et al., 2015; Solomon et al., 2016).

Several key factors govern the climate impact of volcanic aerosols from an eruption: (1) eruption magnitude; (2) injection height and self-lofting; (3) eruption latitude; and (4) dynamical evolution. (1) Eruption magnitude can be characterized by the amount of SO₂ and ash released, which largely controls aerosol formation potential and radiative forcing. Notably, substantial perturbations of stratospheric water vapour and aerosol microphysics have also been reported for the 2022 Hunga Tonga eruption despite modest SO₂ injection (Carn et al., 2022; Zhu et al., 2022). (2) Eruption products injected directly into the stratosphere, as well as volcanic plumes in the upper troposphere that self-loft into the lower stratosphere through radiative heating, can persist much longer than plumes confined to the troposphere (Toohey et al., 2025). (3) Aerosols from tropical eruptions are transported most efficiently via the Brewer–Dobson circulation (Brewer, 1949; Dobson, 1956; Butchart, 2014), whereas mid-latitude eruption aerosols can still reach the tropics through Rossby-wave breaking or transport by the Asian Summer Monsoon Anticyclone (ASMA) (Konopka et al., 2009; Kloss et al., 2021; Wu et al., 2023). (4) UTLS jet streams, cyclones, anticyclones, and stratospheric circulation govern dispersion patterns and dilution rates.

In the Northern Hemisphere summer, the ASMA is the dominant circulation system in the UTLS. Deep convection injects pollutants into the UTLS, where the ASMA's strong anticyclonic circulation acts as a dynamical transport barrier, confining these air masses over Asia during ascent into the stratosphere (e.g., Park et al., 2007; Randel et al., 2010; Fadnavis et al., 2014; Santee et al., 2017; Vogel et al., 2019). However, as the ASMA boundary is not a strict transport barrier, air masses can be exported from the monsoon circulation into the extratropical UTLS (Vogel et al., 2016; Yu et al., 2017). This combination of confinement and export is important for interpreting aerosol dispersion during boreal summer, including the dispersion of volcanic aerosol. Previous work on the 2011 Nabro eruption debated whether its plume reached the stratosphere directly or was lofted by monsoon ascent. Bourassa et al. (2012) proposed that the plume remained in the upper troposphere and was subsequently transported into the stratosphere by large-scale ascent in the monsoon, whereas later studies showed evidence of direct stratospheric injection, independent of monsoon-driven lifting (e.g., Fromm et al., 2013; Vernier et al., 2013).

The mid-latitude Raikoke volcano erupted on 21–22 June 2019. Raikoke injected ∼ 1.5 Tg of SO₂ into the lower stratosphere (Cai et al., 2022; Vernier et al., 2024). Subsequently, ash-driven radiative heating lofted parts of the volcanic plume, raising the volcanic cloud top above 20 km within 4 d (Muser et al., 2020; Gorkavyi et al., 2021). 3 d after the eruption, compact anticyclonic “vorticized volcanic plumes” (VVPs) detached from the main volcanic cloud, trapping more than half of the SO₂ mass (Khaykin et al., 2022). The primary VVP then rose to ∼ 27 km, spanning ∼ 400 km in width but only ∼ 1.5 km in depth. Meanwhile, the residual main plume was diluted at lower altitudes, resulting in two distinct Raikoke plumes (Khaykin et al., 2022). In August 2019, part of the aerosol plume was transported into and through the ASMA during the mature stage of the boreal-summer anticyclone, providing an ideal case to examine how the ASMA modulates volcanic plume transport.

Figure 1 (a) Map showing the locations of the balloon measurement sites over Lhasa (29.7° N, 91.1° E) and Boulder (39.9° N, 105.2° W), and the Raikoke volcano (48.3° N, 153.3° E). The blue contour shows the ERA5 geopotential height at 150 hPa averaged over August 2019 at 00:00 UTC. The 14.32 km contour is used here as an indicator of the ASMA at UTLS levels. (b) Spatial distribution of SO₂ vertical column density (DU) from TROPOMI, retrieved with an assumed plume height of 15 km, using orbits from 22:46 UTC on 24 June to 03:50 UTC on 25 June 2019.

Beginning in August 2019, frequent balloon-borne measurements were conducted over Lhasa, at UTLS levels located within the ASMA as indicated by the ERA5 150 hPa geopotential height field (Fig. 1a). The profiles captured two pronounced aerosol layers at distinct altitudes that are consistent with the trailing filament of the VVP and the diluted main plume. We combine these observations with ERA5-driven CLaMS backward trajectories and global three-dimensional SO₂-based tracer simulations to investigate the source region, transport pathways, and transit times from Raikoke to the Tibetan Plateau. Here, SO₂-based refers to passive tracers initialized using the horizontal extent of the Raikoke SO₂ cloud derived from the Tropospheric Monitoring Instrument (TROPOMI) aboard Sentinel-5P (Theys et al., 2024). We further use POPS profiles from Lhasa and Boulder to evaluate the simulated plume evolution over the Northern Hemisphere. In particular, Boulder serves as an independent site outside the ASMA region to evaluate plume transport. In situ UTLS observations of volcanic aerosol are rare, particularly for the 2019 Raikoke eruption. Here we combine balloon-borne COBALD and POPS profiles with high-resolution ERA5-driven CLaMS backward trajectories and global three-dimensional SO₂-based tracer simulations to identify transport of the Raikoke plume through the Northern Hemisphere and into the ASMA, with particular attention to transport pathways and time scales. Specifically, we address the following research questions: (1) What are the key transport pathways of the Raikoke plume in the vicinity of the ASMA? (2) Which processes dominate plume dilution? In addition, we quantify the sensitivity of the simulations to key model settings (parameterized mixing intensity, injection height, and injection region).

The remainder of the paper is organized as follows: Sect. 2 describes data and methods; Sect. 3 presents balloon observations from Lhasa and Boulder; Sect. 4 reports the backward-trajectory analysis and global CLaMS transport simulations; Sect. 5 discusses modelling sensitivities, including sensitivities to parameterized mixing, injection-region definition, and injection height; Sect. 6 concludes; and Appendix A presents supplementary material.

Table 1 Summary of analyzed balloon flights at Lhasa and Boulder, including sequential flight ID, date, mid-ascent time (UTC), key instruments used in this study, and the potential temperature (θ) and pressure (p) ranges of the observed volcanic plume layer. The plume layer was identified from COBALD measurements where BSR₄₅₅ > 1.1 and CI > 6, and from POPS measurements where the particle number concentration exceeded 150 cm⁻³. Mid-ascent time denotes the midpoint of the ascent period from launch to balloon burst. The “Key instruments” column lists only the instruments analyzed in this work; additional payload components were flown but are not listed here. A dash (–) denotes that no volcanic plume layer was identified for that flight.

2 Data and Methods

2.1 Balloon-borne Instruments

Volcanic aerosol profiles at Lhasa were obtained during the 2019 Asian Summer Monsoon (ASM) season as part of the Sounding Water Vapor, Ozone, and Particles (SWOP) campaign. The SWOP campaign was led by the Institute of Atmospheric Physics, Chinese Academy of Sciences. Balloons were equipped with an electrochemical concentration cell (ECC) ozonesonde, a cryogenic frostpoint hygrometer (CFH), a compact optical backscatter aerosol detector (COBALD), and an International Met Systems (InterMet) iMet-1-RSB radiosonde (GRUAN Lead Centre, 2025). However, ECC ozone measurements are not analyzed in this study. A total of 9 measurements occurred in July–August (1 in July, 8 in August), followed by monthly measurements through January 2020. On 20 August 2019, the payload also included a Portable Optical Particle Spectrometer (POPS). During July–August 2019, Lhasa lay within the interior of the ASMA at UTLS levels (Fig. 1a). The site was chosen to observe the Raikoke plume in the ASMA interior and to document its arrival and vertical evolution. Additional information on SWOP field activities and on campaigns conducted in other years is available in Bian et al. (2012); Li et al. (2017, 2018, 2020); Ma et al. (2022); Yang et al. (2023).

The POPS data over Boulder were provided by the Baseline Balloon Stratospheric Aerosol Profiles (B²SAP) project; we analyze five profiles – 3 influenced by the Raikoke eruption (2 in August 2019 and 1 in November 2019) and 2 background references on 28 June 2019 and 3 December 2019 (Todt et al., 2023).

Figure 1a shows the geographic locations of Raikoke volcano and the balloon-sounding sites at Lhasa and Boulder. An overview of the analyzed balloon flights at Lhasa and Boulder, including the identified plume-layer ranges, is provided in Table 1.

The Compact Optical Backscatter Aerosol Detector (COBALD) was developed by ETH Zurich (Brabec et al., 2012). Its lightweight and portable design makes it suitable for balloon payloads used in the study of cirrus clouds (Brabec et al., 2012; Cirisan et al., 2014; Reinares Martínez et al., 2021; Yang et al., 2023) and aerosols (Vernier et al., 2015; Brunamonti et al., 2018; Hanumanthu et al., 2020). COBALD uses two high-power LEDs emitting at 455 nm (blue) and 940 nm (near‑infrared). It measures backscattered light from air molecules and particles, including aerosols and cloud particles (e.g. ice crystals). The backscatter ratio (BSR) is defined as

$$\begin{array}{}\text{(1)}& \text{BSR}={\displaystyle \frac{{\mathit{\beta }}_{\text{air}}+{\mathit{\beta }}_{\text{particles}}}{{\mathit{\beta }}_{\text{air}}}},\end{array}$$

where β_air is the backscatter coefficient of air molecules and β_particles is the backscatter coefficient of aerosols or cloud particles. The instrument measures only the total backscatter (${\mathit{\beta }}_{\text{total}}={\mathit{\beta }}_{\text{air}}+{\mathit{\beta }}_{\text{particles}}$); β_air is calculated in post-processing from Rayleigh scattering using the measured pressure and temperature profiles, and β_particles is obtained as β_total−β_air. Here, BSR₄₅₅ represents the backscatter ratio at 455 nm. The instrument's field of view is ± 6°, and its detection range spans 0.5–10 m. Because sunlight severely interferes with the measurements, COBALD is only operated at night. COBALD BSR uncertainty is typically characterized by an absolute error interval of about ∼ 5 % for the profile, while the precision is better than ∼ 1 % under UTLS conditions (Vernier et al., 2015; Brunamonti et al., 2021; Reinares Martínez et al., 2021).

The COBALD color index (CI) is defined as

$$\begin{array}{}\text{(2)}& \text{CI}={\displaystyle \frac{{\text{BSR}}_{\mathrm{940}}-\mathrm{1}}{{\text{BSR}}_{\mathrm{455}}-\mathrm{1}}},\end{array}$$

which serves as an empirical proxy for particle size within aerosol layers. Note that when BSR₄₅₅ and BSR₉₄₀ are close to 1, CI can become sensitive to small baseline or calibration offsets.

The Cryogenic Frostpoint Hygrometer (CFH) is a high‐precision instrument for measuring water vapor based on the chilled mirror principle. It controls the mirror temperature to maintain a stable frost- or dew-point. Combined with temperature and pressure data from the iMet radiosonde, the relative humidity over ice (RH_ice) is calculated using the empirical equation from Murphy and Koop (2005):

$$\begin{array}{}\text{(3)}& {\text{RH}}_{\text{ice}}={\displaystyle \frac{e_{\text{ice}}\left(T_{\text{mirror}}\right)}{e_{\text{ice}}\left(T_{\text{environment}}\right)}},\end{array}$$

where e_ice(T) is the saturation vapour pressure over ice at temperature T (in Kelvin), T_mirror is the measured frost-point temperature, and T_environment is the ambient air temperature. CFH uncertainties may be as low as ∼ 2 % in the lower troposphere and ∼ 5 % near the tropical tropopause under good operating conditions (Vömel et al., 2016). For UTLS conditions we use a conservative uncertainty of ± 10 % (Fahey et al., 2014; Poltera et al., 2025).

The Portable Optical Particle Spectrometer (POPS) is a lightweight instrument for measuring aerosol number density and size distribution (Gao et al., 2016). POPS employs a 405 nm laser to detect light scattered by individual particles, with scattering intensity related to particle size. It measures particles from 140 to 2500 nm and reports number concentrations in size bins (Todt et al., 2023). Measurement uncertainties are dominated by particle sizing (including sensitivity to the assumed refractive index) and flow-rate calibration (Gao et al., 2016).

2.2 TROPOMI

TROPOMI, the satellite instrument aboard ESA's sun-synchronous Sentinel-5P platform launched in 2017, is a hyperspectral imaging spectrometer that records back-scattered radiation from the ultraviolet to the shortwave infrared (Veefkind et al., 2012). Total-column SO₂ in Dobson units (1 DU = 2.69 × 10¹⁶ $\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}$) is retrieved using differential optical absorption spectroscopy (DOAS) applied to three wavelength windows (312–326, 325–335, 360–390 nm) (Theys et al., 2017, 2024). The publicly available Level 2 product provides four vertical columns: the surface–to–top-of-atmosphere column and three columns assuming an SO₂ plume centered at 1, 7, or 15 km altitude. Following previous Raikoke studies (Muser et al., 2020; de Leeuw et al., 2021; Cai et al., 2022), we adopt the 15 km retrieval, which most closely matches the eruption's mean injection height; values below the instrument's 0.3 DU detection limit are excluded, following Theys et al. (2024). Four consecutive TROPOMI overpasses of the SO₂ vertical column density at 15 km altitude captured the entire Raikoke SO₂ cloud between 22:46 UTC on 24 June and 03:50 UTC on 25 June 2019 (Fig. 1b).

2.3 CLaMS

The Chemical Lagrangian Model of the Stratosphere (CLaMS) is a chemistry–transport model that calculates the three-dimensional motion of air parcels, which in turn represent the model grid (McKenna et al., 2002; Pommrich et al., 2014). The CLaMS hybrid vertical coordinate (ζ) follows orography near the surface and transitions smoothly into θ once σ = $p/p_{\text{surf}}$ reaches 0.3 (usually about 300 hPa) (Pommrich et al., 2014). The isentropic coordinate in CLaMS is well suited for representing transport and mixing in the stratosphere. All CLaMS applications in this study – both backward trajectories and three-dimensional simulations with SO₂-based tracers – are driven by ERA5 reanalysis: the standard runs use the native hourly fields on the 0.3° × 0.3° (latitude × longitude) grid with 137 vertical levels up to 80 km (Hersbach et al., 2020), whereas the sensitivity experiment in Appendix Fig. A2 repeats the simulation with ERA5 fields at 1° × 1° spatial and 6 h temporal resolution.

Diabatic backward trajectories are initialized every second along the balloon's vertical ascent profile, using the in-situ measurements of temperature (for deriving potential temperature), pressure, time, longitude, and latitude to define the start positions. Vertical velocities were computed from ERA5 wind fields and total diabatic heating rates, including clear-sky and cloud radiation, latent heat release, and turbulent and diffusive transport (Ploeger et al., 2021). This approach captures uplifts by resolved convection as represented in ERA5 (Li et al., 2020; Clemens et al., 2024; Vogel et al., 2024). To identify the volcanic influence, only those backward trajectories that passed through the Raikoke eruption region (Fig. 1b) during 22:46 UTC on 24 June and 03:50 UTC on 25 June 2019 were retained.

Global three-dimensional transport simulations use passive SO₂-based tracers initialized from the horizontal extent of the Raikoke SO₂ cloud in the TROPOMI 15 km product (Fig. 1b; four overpasses between 22:46 UTC on 24 June and 03:50 UTC on 25 June 2019; values ≥ 0.3 DU) (Theys et al., 2024). At 00:00 UTC on 25 June 2019, parcels inside this mask are initialized to 1 (0 outside) within 3 potential-temperature ranges (380–400, 400–420, and 420–440 K) and then advected and mixed globally, yielding tracer fractions between 0 and 1.

Table 2 Summary of tracer simulations using different injection-region definitions (satellite SO₂ plume mask and a rectangular domain), including injection layer, mixing interval, and critical Lyapunov exponent (λ_c).

To assess how different mixing intensities influence the reconstruction of volcanic plume transport processes, two simulations were conducted: (i) a control simulation with mixing every 24 h and λ_c=1.5 (a commonly used set-up for CLaMS simulations driven by high-resolution ERA5 reanalysis data; Vogel et al., 2025); (ii) a modified simulation with mixing every 6 h and λ_c=3.5. In the CLaMS mixing scheme, parameterized mixing is triggered when the integral deformation between neighboring air parcels exceeds an empirical critical deformation γ_c=λ_cΔt, where λ_c is the critical Lyapunov exponent and Δt is the advective time step (mixing interval). For a given Δt, a larger λ_c (thus a larger γ_c) requires stronger deformation to trigger mixing and therefore corresponds to less frequent parameterized mixing (and vice versa) (Konopka et al., 2004, 2007). In our setup, the modified simulation corresponds to enhanced parameterized mixing compared to the control simulation. Throughout most of the paper we show results from the modified simulation, as these agree better with the observations. Sensitivity to parameterized mixing intensity and comparisons with the control simulation are discussed in Sect. 5.1. The mixing configurations and the additional sensitivity runs (rectangular mask and coarser ERA5 input) are summarized in Table 2.

Figure 2 COBALD and CFH profiles above Lhasa from 30 July 2019 to 4 January 2020. COBALD backscatter ratio at 455 nm (BSR₄₅₅, blue) and COBALD color index (CI, black) are shown together with CFH relative humidity over ice (RH_ice, magenta). Pressure is given on the left axis and the corresponding potential temperature on the right axis. Orange shading marks the Raikoke influenced aerosol layers used to initialize the backward trajectories. Light blue shading indicates lower altitude aerosol enhancements attributed to ATAL that are not used for the plume initialization. Gray shading indicates cirrus conditions.

3 Measurement Results

The COBALD–CFH tandem provides a practical method for identifying cirrus clouds and aerosol layers because collocated backscatter and frost-point humidity measurements allow a straightforward separation of ice-cloud signatures from aerosol enhancements (Brabec et al., 2012; Cirisan et al., 2014; Reinares Martínez et al., 2021; Yang et al., 2023). Cirrus clouds are identified using the criteria BSR₄₅₅>1.2, RH_ice > 70 % and CI > 7 (Vernier et al., 2015; Brunamonti et al., 2018; Hanumanthu et al., 2020; Yang et al., 2023). These clouds are shown as gray-shaded regions in Fig. 2. On 30 July 2019, RH_ice reached ∼ 90 % at 423 K (70 hPa), but no corresponding enhancement in BSR₄₅₅ was observed; thus, this layer is not classified as cirrus by our joint criterion and likely represents a moist, ice-free layer. When cirrus and aerosols occur within the same profile, aerosol and cloud signals can generally be separated if they are vertically distinct (e.g., using BSR and CI), whereas if a cirrus cloud overlaps the aerosol layer at the same altitude, the aerosol contribution cannot be reliably isolated because cirrus BSR₄₅₅ typically dominates. Thus, the aerosol cannot be reliably detected or quantified under such conditions. Using these criteria, we exclude cirrus-contaminated layers. Remaining enhancements with BSR₄₅₅>1.1 are treated as aerosol layers and then classified using a CI threshold of 6: layers attributed to the Raikoke plume (CI > 6) are highlighted in orange in Fig. 2, while layers more consistent with ATAL aerosol (CI < 6) are highlighted in light blue. We note that CI is used here as an additional indicator within aerosol layers and should be interpreted with caution, and that the ATAL identification in our dataset is limited. In very clean air where BSR values at both wavelengths approach 1, CI can become unstable; we therefore apply CI only within aerosol layers with BSR₄₅₅>1.1 after excluding cirrus. On 3 August 2019, 42 d after the Raikoke eruption, two distinct BSR₄₅₅ peaks were observed: the upper peak at 459–472 K (59–56 hPa) reached the maximum value recorded in Lhasa (∼ 1.8), and the lower peak at 410–422 K (78–71 hPa) had a somewhat weaker magnitude (∼ 1.4). On 6 August 2019, another sharp peak appeared at 404–416 K (82–75 hPa). After 8 August 2019, the peaks became less sharp and more vertically extended, which persisted until 24 November 2019. By 4 January 2020, the aerosol backscatter enhancement had further weakened, and no distinct peak was observed in the profile.

Figure 3 (a) COBALD measurements in Lhasa. The black line represents the median BSR₄₅₅ for the 2019 profiles influenced by the volcanic event, and the gray lines show the corresponding individual profiles. The red line shows the median BSR₄₅₅ during the ATAL in 2013. The blue line corresponds to the BSR₄₅₅ on 30 July 2019, which was not affected by the Raikoke eruption. (b) POPS measurements in Boulder and Lhasa.

In general, ATAL-related enhancements in COBALD BSR₄₅₅ profiles are largely confined to 360–400 K (core near 370–390 K), with occasional extensions up to 420–440 K depending on region and year (Vernier et al., 2015, 2018; Appel et al., 2022). The ATAL profile from 2013 shown in Fig. 3a is taken from the COBALD measurements over Lhasa reported by Vernier et al. (2015). The 2019 median peak occurs near 417 K, about 33 K above the 2013 ATAL peak (384 K), corresponding to roughly 1.7 km in altitude in the UTLS (based on the 30 July 2019 background sounding). The vertical extent of ATAL enhancements varies across regions and years; therefore, we emphasize the peak magnitude as the more robust difference: the 2019 Raikoke-related median BSR₄₅₅ reaches ∼ 1.25, exceeding the 2013 ATAL peak (∼ 1.10; Fig. 3a).

All POPS particle number densities are converted to standard temperature and pressure (STP; 1013 hPa, 273.15 K) to remove the pressure-driven decline with altitude (Fig. 3b). At Boulder, no clear volcanic aerosol signal was present on 28 June 2019, and by 3 December 2019 the profile shows no obvious volcanic aerosol peak (Fig. 3b). The particle number density at Lhasa on 20 August 2019 was lower than that at Boulder on 7 and 27 August 2019. This is consistent with the global CLaMS tracer simulations (Sect. 4.2; cf. Figs. 7 and 8), which show that the main Raikoke plume is transported along the subtropical westerly jet in the mid-latitude UTLS, while transport into the ASMA interior is more limited due to the anticyclone acting as a partial horizontal barrier. Section 4 presents backward-trajectory calculations and three-dimensional tracer simulations with CLaMS to verify the origin and transport pathways of the plume observed over Lhasa and Boulder.

4 Results of CLaMS transport simulations

4.1 Backward-Trajectory Calculations

To verify that the enhanced aerosol layers observed over Lhasa are linked to the Raikoke eruption, we calculated ERA5-driven diabatic backward trajectories with CLaMS (Fig. 4). Trajectories were released every second along each balloon ascent within the selected plume-layer potential-temperature intervals (orange shading in Fig. 2), using the iMet-derived time and location (longitude, latitude, pressure, and temperature) as initial conditions. All trajectories were calculated backward to a single common reference time (21 June 2019, 18:00 UTC), which serves as the trajectory endpoint for all cases.

Figure 4 Backward trajectories initialized from the Lhasa balloon observations and traced back to the Raikoke eruption reference time (21 June 2019, 18:00 UTC). Only trajectories initialized within the enhanced BSR₄₅₅ potential-temperature range (orange shading in Fig. 2) are shown. Colors indicate trajectory longitude. The x-axis shows date (dd.mm.yy) along the backward integration (right: balloon initialization time; left: eruption reference time).

Figure 5 Backward trajectories for the Lhasa flights (3, 6, 8, 10, and 12 August 2019). Only trajectories that intersect the TROPOMI SO₂ mask between 24 June 2019, 22:46 UTC and 25 June 2019, 03:50 UTC are shown. Trajectories are colored by potential temperature, and the SO₂ mask is shaded in orange. Panel titles give the start day, start θ range, and the fraction of trajectories that satisfy the SO₂ mask criterion. Locations of Lhasa and Raikoke are marked.

Pronounced BSR₄₅₅ peaks were observed over Lhasa on 3 and 6 August 2019 (Fig. 2). For the upper layer on 3 August 2019, backward trajectories indicate that the air parcels were transported approximately along isentropic surfaces (Fig. 4). By contrast, the lower layer on 3 August 2019 includes trajectories lofted from ∼ 370 to ∼ 410 K. In the upper layer on 3 August 2019, the transport pathways differed significantly from other layers, as evidenced by changes in longitude (see Fig. 5 for trajectory details). This suggests that the enhanced BSR₄₅₅ resulted from volcanic aerosols originating from two distinct branches of the volcanic plume at different altitudes, each following a different transport pathway. From 8–20 August 2019, mixing from lower potential-temperature levels began to occur. This coincided with a decrease in BSR₄₅₅ values in the in-situ measurements, as air with relatively low aerosol concentrations dilutes the enhanced signals. Over the following three months, as the ASMA weakened seasonally, air from lower potential-temperature levels increasingly influenced the Lhasa profiles. During this period, peak BSR₄₅₅ within the identified plume layer (orange shading in Fig. 2) decreased from values up to 1.88 in early August 2019 (95th percentile up to 1.77) to ≤ 1.25 from 30 September to 24 November 2019 (95th percentile ≤ 1.24), indicating progressive dilution of the volcanic aerosol layer by relatively aerosol-poor air from the lower troposphere.

To illustrate air masses influenced by the Raikoke eruption, we filtered backward trajectories to include only those passing through the eruption region (Fig. 1b) between 22:46 UTC on 24 June 2019 and 03:50 UTC on 25 June 2019. We then examined several representative pathways to show how the volcanic plume entered the ASMA interior (Fig. 5). In Fig. 5, the reported fraction is defined as (number of trajectories that meet the eruption-region criterion)/(total number of backward trajectories initialized for the corresponding flight date and θ layer shown in Fig. 4). For the 6 cases shown in Fig. 5, the total numbers of initialized trajectories are 79 and 50 (3 August 2019, two θ layers), 72 (6 August 2019), 190 (8 August 2019), 165 (10 August 2019), and 105 (12 August 2019). Because this filtering criterion is highly selective, only a small fraction of trajectories remains, and the resulting fractions can be regarded as conservative estimates. A sensitivity test shows that the fraction of backward trajectories reaching the eruption region is ∼ 3 %–10 % when using the original TROPOMI SO₂ footprint and satellite-pass time window, but increases to ∼ 12 %–28 % when applying a broader rectangular domain (42–73° N, 137–215° E) and extending the time window to 7 d. The retained trajectories still follow the same main transport pathways shown here, providing confidence in the identified transport patterns.

For the calculation on 3 August 2019, the trajectories are divided into two distinct branches. The trajectories in the lower branch originate from 410–422 K and correspond to the main volcanic aerosol plume carried by the subtropical westerly jet. The upper branch originates from 459–472 K. Satellite observations from TROPOMI on Sentinel-5P indicate that the VVP core was entrained into the summertime easterlies around 20–25 July (Gorkavyi et al., 2021; Khaykin et al., 2022). The backward trajectories initialized in the 3 August 2019 plume layer (459–472 K) in Fig. 5 are consistent with the late-July phase of this satellite-tracked pathway. The potential temperature of the VVP during its transit through the ASMA inferred from satellite detections also closely matches the altitudes of enhanced BSR₄₅₅ (Gorkavyi et al., 2021; Khaykin et al., 2022). Meanwhile, lidar measurements in Wuhan (central China) on 30 July reveal that the VVP's core backscatter coefficient is much greater than that of the main aerosol plume (Jing et al., 2023; He et al., 2024), implying that our measurements in this branch captured only the trailing filament of the VVP. Although both sets of air parcels belong to the same balloon profile, they trace back to two distinct plume components: the lower-θ main plume in the subtropical westerly jet (∼ 410–422 K) and the higher-θ VVP filament in the summertime easterlies (∼ 459–472 K). This indicates different source altitudes within the Raikoke plume and distinct dynamical transport pathways.

From 6 to 20 August 2019, the overall transport pattern remained similar to that on 6 August 2019 – corresponding to the main volcanic aerosol plume primarily driven by the subtropical westerly jet. After entering the ASMA, some air parcels deviate from the jet pathway and circulate within the anticyclonic region. This increased spreading of trajectories suggests dispersion and lateral mixing with surrounding air, which may contribute to the observed decrease in BSR₄₅₅.

Figure 6 Same as Fig. 5, but for the Boulder flights (7 and 27 August 2019). For 27 August 2019, the fraction is computed for 385–429 K, while the plotted trajectories are restricted to 405–410 K to better illustrate the transport pathways. Locations of Boulder, Lhasa, and Raikoke are marked.

Backward trajectories for the Boulder flights on 7 and 27 August 2019 are shown in Fig. 6. On 7 August 2019, an enhanced aerosol layer is observed between 373 and 442 K (particle number concentration > 150 cm⁻³). The corresponding backward trajectories indicate transport to Boulder primarily along the subtropical westerly jet. On 27 August 2019, two volcanic layers are present. The upper layer (480–505 K) is associated with the VVP filament transported in summertime easterlies. The lower layer occurs at 385–429 K. To better illustrate the transport pathway, we plot only trajectories initialized between 405 and 410 K. These trajectories indicate that the plume had already traveled once around the globe before reaching Boulder.

Figure 7 SO₂-based tracer fractions within the 400–420 K layer from the modified simulation (simulation details in Table 2). The purple triangle indicates Lhasa. The ASMA edge is shown by the black contour at the 378 × 10³ m² s⁻² Montgomery streamfunction. The white contour marks the 11 PVU potential vorticity isoline.

4.2 Global three-dimensional CLaMS simulations with SO₂-based tracers

In global three-dimensional CLaMS simulations, we use passive tracers initialized from the observed SO₂ plume mask (Fig. 1b). We refer to these as SO₂-based tracers. At 00:00 UTC on 25 June 2019, the tracer is set to 1 inside the mask (0 outside) within 3 potential-temperature ranges (380–400, 400–420, and 420–440 K) and then advected and mixed throughout the global atmosphere in CLaMS. We present results for the 400–420 K release, while the 380–400 and 420–440 K releases are used as sensitivity tests (Sect. 5.2; Figs. 11–12). Figure 7 shows global maps of fractions of the SO₂-based tracer falling within the 400–420 K layer for each measurement day from 30 July to 20 August 2019, using the modified simulation. Here, the tracer fraction (0–1) represents the fractional contribution of air originating from the volcanic plume (identified by the TROPOMI SO₂ mask shown in Fig. 1b) after advection and mixing.

Figure 8 Same as Fig. 7, but for SO₂-based tracer fractions within the 400–420 K layer over the North American continent. The black triangle indicates the location of Boulder.

The ASMA is indicated by the Montgomery streamfunction (Santee et al., 2017) at the 378 × 10³ m² s⁻² contour (layer mean over 400–420 K) as an estimate of its edge. The 11 PVU potential vorticity contour (layer mean over 400–420 K) illustrates the isentropic transport, including ASMA outflow (filaments) and in-mixing of PV-rich air from high latitudes. The strong transport barrier of the ASMA observed within ∼ 370–390 K is much weaker above 400 K (Ploeger et al., 2015). The results shown in Fig. 7 illustrate the overall intrusion of air masses into the ASMA and agree well with the Lhasa COBALD BSR₄₅₅ observations. From 30 July to 1 August 2019, part of the SO₂-based tracer entered the ASMA from its eastern side without reaching Lhasa. This is consistent with the absence of clear aerosol signals in the lower layers of the balloon-borne measurement data. On 3 and 6 August 2019, air masses impacted by the SO₂-based tracer arrived above Lhasa. Without a significant effect of upwelling from below at that time, the SO₂-based tracer signals remained sharp and concentrated. From 8 to 20 August, most of the ASMA contained mixed air masses. However, SO₂-based tracer fractions inside the ASMA remained lower than outside, reflecting the anticyclonic transport barrier. During August, only a minor fraction of the SO₂-based tracer is located inside the ASMA, while the majority remains outside. By later months, progressive mixing reduces the contrast between ASMA and surrounding regions in the tracer distribution (not shown here).

The corresponding SO₂-based tracer fractions within the 400–420 K layer over the North American continent are shown in Fig. 8. Elevated tracer fractions are simulated over Boulder on both 7 and 27 August 2019. By 8 November 2019, the overall tracer fractions are markedly lower, reflecting progressive dilution by mixing during transport.

5 Modelling sensitivities

5.1 Sensitivity to parameterized mixing and injection-region definition

To assess the sensitivity of plume dispersion to model settings, such as parameterized mixing intensity, tracer release regions, and release altitudes, we performed a series of sensitivity experiments. In previous work with CLaMS driven by 1° × 1° resolution and 6 h temporal resolution reanalysis data (e.g., ERA5, ERA-Interim, as applied in Kloss et al., 2021), the intensity of parameterized small-scale mixing was controlled by choosing a critical Lyapunov exponent λ_c=1.5 and a mixing interval of 24 h (e.g., Pommrich et al., 2014). Here, we performed simulations with CLaMS driven by high-resolution ERA5 data (0.3° × 0.3° resolution, 1 h) and found best agreement with observations by enhancing mixing intensity in the simulation slightly by choosing λ_c=3.5 and a mixing interval of 6 h. Taken together, the modified setup corresponds to enhanced parameterized mixing compared to the control setup.

To evaluate the sensitivity of our results with respect to parameterized mixing intensity, Figure A1 in the Appendix shows the SO₂-based tracer fractions from the control simulation (λ_c=1.5; mixing interval Δt = 24 h) and compares them with those from the modified simulation (λ_c=3.5; Δt = 6 h) in Fig. 7. The control simulation in Figure A1 shows a similar overall intrusion pattern but lacks the distinct fractions over Lhasa on 3 August 2019. The tracer distribution in the control simulation is more fragmented, presenting small, isolated pockets compared to the more homogeneous structures in the modified simulation. In a further sensitivity test, we repeat the control simulation with ERA5 fields at coarser spatial and temporal resolution (details in Table 2). This setup reduces computing time and is often used for global CLaMS simulations that cover multiple years. However, changing the ERA5 resolution has only a minor effect and does not affect our conclusions. This is evident from the comparison of Fig. A1 (0.3° × 0.3°, 1 h) and Fig. A2 (1° × 1°, 6 h).

Figure 9 Comparison of tracer-fraction profiles with the BSR₄₅₅ profile over Lhasa. SO₂-based tracer fractions from the control simulation (orange) and the modified simulation (blue) are shown together with rectangle-based tracer fractions from the modified simulation (red); see Table 2 for simulation details. Annotation colors match the corresponding simulation profiles. Values give the Pearson correlation r and the absolute peak-height offset |Δθ| (K) between the model tracer-fraction profile and the COBALD BSR₄₅₅ profile within the analyzed θ range (375–450 K; 375–475 K on 30 September, 28 October, and 24 November 2019).

The different Raikoke tracers from the global three-dimensional CLaMS simulations are interpolated in time and space to the balloon-borne vertical profiles, providing a tracer-fraction profile at each measurement time. By comparing these tracer-fraction profiles from the model simulations with different mixing intensities with the observations (Fig. 9), we assess how the chosen parameterized mixing intensities affect the accuracy of the simulated plume dispersion. To quantify model–measurement agreement for each profile, we compute (i) the Pearson correlation coefficient r between the area-normalized tracer-fraction profile and the COBALD aerosol backscatter ratio at 455 nm (ABSR₄₅₅), defined as BSR₄₅₅−1 (e.g., Cirisan et al., 2014; Reinares Martínez et al., 2021), over the main plume layer (375–450 K; 375–475 K on 30 September, 28 October, and 24 November 2019) and (ii) the absolute peak-height difference |Δθ| between the two profiles. Note that r reflects only linear association and should be interpreted as a relative indicator rather than a definitive measure of agreement. The values of r and |Δθ| are annotated in Fig. 9. In Fig. 9, the SO₂-based tracers are also released at 400–420 K. This setup primarily samples the diluted main plume rather than the higher-altitude trailing filament of the vorticized volcanic plume (VVP). Accordingly, the higher-altitude peak on 3 August 2019 likely originates from the VVP filament near 460–490 K and is not resolved by the 400–420 K release. A sensitivity simulation targeting the upper-level VVP filament shows enhanced tracer fractions near the observation region (not shown here); however, with the currently used CLaMS resolution, this thin and sharp aerosol peak at a single profile location is difficult to reproduce. Consistent with Figs. 7 and A1, the modified simulation produces tracer peaks (blue lines) that closely match the balloon measurements – except on 6 and 8 August 2019, when the peak altitudes deviate slightly. For most dates, the modified simulation shows consistently higher r and smaller |Δθ|. Although the control simulation (orange lines) occasionally aligns with observed peak altitudes, its overall agreement is weaker and it shows large peak-altitude biases on 12 August 2019 and 24 November 2019 (|Δθ_peak| ≈ 25 and ≈ 36 K, respectively; Fig. 9). Overall, the modified simulation captures the vertical structure of the BSR₄₅₅ profile well, despite minor differences in peak altitude for certain dates.

Furthermore, we assess the sensitivity of the model simulations to the plume initialization by carrying out another simulation using the modified configuration, in which Raikoke plume tracers are initialized within a rectangular latitude–longitude mask (from 49 to 62° N and 163° E to 170° W), following Kloss et al. (2021) (see Table 2 for simulation details). Tracers were released on multiple isentropic levels within this rectangular domain. Here, we show the rectangle-based tracer fractions, released at potential‐temperature levels of 400–420 K (red lines in Fig. 9). The peak heights of the rectangle-based tracer fractions (red lines) and the SO₂-based tracer fractions (blue lines) show no significant difference. Because the release domain captures the plume core, even a simple rectangular mask reproduces the observed peak altitude. Hence, for the Raikoke simulations the choice of injection region is less critical than the choice of model mixing intensity.

Figure 10 Same as Fig. 9 but for Boulder: comparison of tracer-fraction profiles with the POPS aerosol number-concentration profile at STP. Annotations report r and |Δθ| (K) relative to POPS within the analyzed θ range (375–450 K for 7 and 27 August 2019, and 375–475 K for 8 November 2019).

Over Boulder, we likewise compare the tracer-fraction profiles on the measurement dates with the POPS aerosol number concentration profiles at STP (Fig. 10). For the computation of r, we evaluate within 375–450 K for 7 and 27 August 2019, and within 375–475 K for 8 November 2019. Overall, the SO₂-based tracers in the modified simulation (blue) reproduce the observed peak structure more closely, with higher r and smaller |Δθ| than the control and rectangle-based simulations.

Figure 11 SO₂-based tracer‐fraction profiles from the modified simulation at different isentropic injection levels (simulation details in Table 2), compared with COBALD BSR₄₅₅ profiles. Blue lines: injections at 400–420 K; brown lines: 380–400 K; purple lines: 420–440 K. Annotation colors match the corresponding injection-layer curves; values are r and |Δθ| (K) between the model tracer-fraction profile and the COBALD BSR₄₅₅ profile within the analyzed θ range (375–450 K; 375–475 K on 30 September, 28 October, and 24 November 2019).

5.2 Sensitivity to injection height

The Raikoke plume was injected primarily within 8–16 km (Kloss et al., 2021; Cai et al., 2022; Vernier et al., 2024). The main injection occurred between 18:00 UTC on 21 June and 03:00 UTC on 22 June 2019 (with SO₂ release possibly extending to ∼ 06:00 UTC) (Cai et al., 2022; Vernier et al., 2024). Radiative heating of ash and nascent sulfate then lofted parts of the plume, raising volcanic cloud tops above 20 km within 4 d of the eruption (Muser et al., 2020; Gorkavyi et al., 2021). These observations indicate that Raikoke's plume experienced an intense updraft during the eruption stage. Aerosol‑radiative lofting is not included in the CLaMS simulations. However, comparison of the vertical profiles of the SO₂-based tracers with in-situ measurements allows estimation of the lofting height. SO₂-based tracers were injected at 20 K intervals between 360 and 500 K to cover the entire possible extent of the Raikoke plume. To assess which initialization range best reproduces the observed plume, we compare tracer‑fraction profiles with the BSR₄₅₅ profiles in Fig. 11.

Figure 12 Same as Fig. 11 but for Boulder: SO₂-based tracer-fraction profiles from the modified simulation at different isentropic injection levels compared with the POPS aerosol number-concentration profile at STP. Blue: 400–420 K; brown: 380–400 K; purple: 420–440 K. Annotations report r and |Δθ| (K) relative to POPS within the analyzed θ range (375–450 K for 7 and 27 August 2019, and 375–475 K for 8 November 2019).

Overall, tracers initialized within 400–420 K (blue lines; ∼ 15–16.5 km at the Raikoke site) best match the enhanced BSR₄₅₅ profiles, compared with 380–400 K (brown) and 420–440 K (purple). For most dates, the 400–420 K initialization attains the highest r and the smallest |Δθ| among the three layers (see annotations in Fig. 11). For instance, on 3, 6, and 12 August 2019, tracers initialized in the 400–420 K layer contribute most strongly to the peak of the BSR₄₅₅ profile, while those from the 380–400 and 420–440 K layers also produce peaks at the same altitude but with smaller amplitudes. On 8 August, 30 September, and 28 October 2019, only the 400–420 K tracer reproduces the BSR₄₅₅ peak at the correct altitude, whereas the peaks from the 380–400 and 420–440 K tracers appear below or above the observed height. For Boulder (Fig. 12), the 400–420 K initialization likewise shows the best agreement with the observations. On 7 August 2019 the 380–400 K layer, and on 27 August 2019 the 420–440 K layer, produce peaks that disagree strongly with the observations (r<0.1).

Cai et al. (2022) estimate a peak injection altitude near 11 km, whereas the tracer level that best matches the observed profiles corresponds to ∼ 15–16.5 km (400–420 K). We hypothesize that this ∼ 4–5 km offset could be consistent with radiative self-heating lofting of the SO₂-rich cloud, but it may also reflect uncertainties in injection-height estimates across techniques and model-related discrepancies (e.g., Kloss et al., 2021; Cai et al., 2022; Vernier et al., 2024).

6 Conclusions

Motivated by the limited availability of in-situ UTLS measurements of volcanic plume and the need to better quantify ASMA-related transport pathways, we combined in-situ balloon observations over Lhasa and Boulder with Lagrangian transport simulations using the CLaMS model, driven by high-resolution ERA5 data, to investigate the transport pathways and mixing processes of the Raikoke plume in the upper troposphere–lower stratosphere (UTLS) during the mature phase of the Asian Summer Monsoon Anticyclone (ASMA). Our main conclusions are:

• Analysis of backward trajectories. 42 d after the Raikoke eruption, the enhanced BSR₄₅₅ was first detected over Lhasa. Two separate transport pathways carried the Raikoke plume to the Tibetan Plateau: (i) a high-altitude pathway along summertime easterly winds (∼ 460–490 K) transporting the trailing filament of the vorticized volcanic plume (VVP) and (ii) a lower-level route via the subtropical westerly jet (∼ 390–430 K) transporting the main volcanic plume. These branches converge over Lhasa on 3 August 2019 but at different altitudes, explaining the dual peaks observed in balloon-borne backscatter profiles.

• Dilution via mixing with relatively aerosol-poor air. During transport through the lower stratosphere, the plume is progressively diluted by background air with low aerosol concentrations. After entering the ASMA, dilution may be driven by upwelling from lower potential temperature levels and by dispersion and lateral mixing with surrounding air. These processes are consistent with the progressively smoothed peaks of enhanced BSR₄₅₅ signals observed along the balloon profiles.

• Validation of CLaMS SO₂-based tracers. The modified CLaMS simulations using SO₂-based tracers successfully capture transport into the ASMA and reconstruct the vertical structure observed over Lhasa, as well as over the North American continent in Boulder. Compared with the control simulation (mixing every 24 h), the modified simulation (mixing every 6 h) produces more coherent tracer distributions that closely match the enhanced BSR₄₅₅ peaks.

In summary, by combining observations and modelling, we provide further insights into how the mid-latitude Raikoke plume is transported through the Northern Hemisphere UTLS over long distances. We show evidence that the Asian summer monsoon plays a key role in this transport. In particular, parts of the plume are entrained into the monsoon anticyclone and are subsequently diluted by horizontal transport and vertical upwelling within the anticyclone. Our findings improve the understanding of how aerosols from a mid-latitude volcanic injection are transported and dispersed in the lower stratosphere, including in the vicinity of the Asian summer monsoon. This improved understanding of plume transport and dilution is relevant for estimating regional radiative forcing and assessing potential impacts on stratospheric ozone chemistry.

We identified two distinct isentropic transport pathways and assessed plume dilution via mixing with air masses outside the plume. Our simulations show that atmospheric mixing processes, as parameterized in CLaMS in relation to flow deformation, are critical for dispersing the volcanic plume during stratospheric transport. Accordingly, changes in model resolution require adjustment of the mixing parameterization. Our results suggest that a small fraction of the Raikoke plume becomes entrained into the ASMA. Within the ASMA, confinement and summertime diabatic uplift can potentially transport these plume fractions to higher altitudes (e.g., Vogel et al., 2019).

Appendix A

Figure A1 Same as Fig. 7, but for SO₂-based tracer fractions from the control simulation (simulation details in Table 2), driven by ERA5 at 0.3° × 0.3° spatial and 1 h temporal resolution.

Figure A2 Same as Fig. A1, but driven by ERA5 at 1° × 1° spatial and 6 h temporal resolution.