Water vapour isotope anomalies during an atmospheric river event at Dome C, East Antarctica

From 19 to 23 December 2018, an atmospheric river sourced in the Atlantic hit the French–Italian Concordia station, located at Dome C, East Antarctic Plateau, 3233 m above sea level (a.s.l.). It induced a significant surface warming (+18 °C in 3 d), combined with high specific humidity (3-fold increase in 3 d) and a strong isotopic anomaly in water vapour (+17 ‰ for δ¹⁸O). The isotopic composition of water vapour monitored during the event may be explained by the isotopic signature of long-range water transport, and by local moisture uptake during the event. In this study, we used continuous meteorological and isotopic water vapour observations, together with the atmospheric general circulation model LMDZ6iso, to describe this event and quantify the influence of each of these processes. The presence of mixed-phase clouds during the event induced a significant increase in downward long wave radiation, leading to high surface temperature and resulting in high turbulent mixing in the boundary layer. Although surface fluxes are underestimated in LMDZ6iso, near-surface temperature and specific humidity are well represented. The surface vapour δ¹⁸O is accurately simulated during the event, despite an overestimated amplitude in the diurnal cycle outside of the event. Using the LMDZ6iso simulation, we perform a surface water vapour mass budget by decomposing total specific humidity into contributions from individual processes. Our analysis demonstrates that surface sublimation, which becomes significantly stronger during the event compared to typical diurnal cycles, emerges as the dominant driver of the vapour δ¹⁸O signal at the peak of the event, accounting for approximately 70 % of the total contribution. The second largest contribution comes from moisture input via large-scale advection associated with the atmospheric river, accounting for approximately 30 % of the total. Consequently, our results reveal that the isotopic signal monitored in water vapour during this atmospheric river event reflects both long-range moisture advection and interactions between the boundary layer and the snowpack. Only specific meteorological conditions driven by a pronounced moisture intrusion can explain these strong interactions. Given the marked imprint of air–snow exchanges on the vapour isotopic signal, improving the representation of local processes in climate models could substantially improve the simulation of the isotopic signal over Antarctica and provide valuable insight into moisture uptake processes.

1 Introduction

Water isotopes are largely used as a proxy of past climate variations, in particular at Dome C, East Antarctica, where the longest continuous climatic record has been obtained from water isotope measurements in a 3200 m long ice core (EPICA community members, 2004; Landais et al., 2021). Interpreting water isotopic signal in ice cores is however complex due to the impact of the atmospheric water cycle before snowfall and of post-deposition processes. Studying variations of δ¹⁸O in atmospheric vapour, precipitations and surface snow is thus critical to understand how large-scale moisture intrusion and local processes interact, to finally produce the signal recorded in ice cores. This is why, since 2011, measurements of the isotopic composition of surface snow and precipitation have been carried out at this location (Touzeau et al., 2016; Stenni et al., 2016; Casado et al., 2021; Dreossi et al., 2024; Ollivier et al., 2025b). This motivation also led to the installation of an analyser of water vapour isotopic composition at Dome C for the summer seasons, first in 2014 (Casado et al., 2016), and for all summers since December 2018 (Leroy-Dos Santos et al., 2021; Ollivier et al., 2025a). The measurement of water vapour isotopic composition is a challenge in such low humidity environment because of the strong influence of mixing ratio on δ¹⁸O values (Leroy-Dos Santos et al., 2021).

Between 19 and 23 December 2018, a +17.6 °C temperature anomaly was observed at the Automatic Weather Station (AWS) of French–Italian Concordia station (123.4° E, 75.1° S, 3233 m above sea level (a.s.l.)), located at Dome C, on the East Antarctic Plateau (Fig. 1), relative to the December 2018 clear-sky daily cycle averages. On 21 December 2018, the temperature reached −14.7 °C, the second highest value recorded since the installation of the Concordia AWS in 2005 (Grigioni et al., 2022). During this event, humidity levels increased by a factor of three and water vapour δ¹⁸O increased by 17 ‰ in 3 d. The primary cause of these anomalous values was the intrusion of warm and moist air from lower latitudes in the form of an atmospheric river (AR), a narrow long band of enhanced water vapour fluxes originating from the sub-tropics and mid-latitudes (Zhu and Newell, 1998; Nash et al., 2018). This AR was detected using the polar-specific detection algorithm described in Wille et al. (2021), and characterized using 10 d back-trajectory analyses performed every 3 h at 1000 m above ground level (a.g.l.) at Dome C with the FLEXPART Lagrangian Particle Dispersion Model (Pisso et al., 2019; Leroy-Dos Santos et al., 2020; Wille et al., 2022; Gorodetskaya et al., 2023). The event was first identified over Coats Land, west of Dronning Maud Land, on 17 December 2018, before progressing towards the South Pole and Dome C (Fig. 1).

Figure 1 FLEXPART 10 d backtrajectory analysis showing particles total footprint, i.e. the density partition of particles released from Concordia station (red dot, 123.4° E, 75.1° S, 3233 m a.s.l.) 20 December 2018 18:00 at a 1000 m a.g.l.. The French–Italian Concordia station is located at Dome C on the East Antarctic Plateau. The red circles and the “−xd” labels indicate the mean position of the particles during the xth day preceding 20 December 2018. The light blue pie charts correspond to the Dronning Maud Land and Adelie Land regions. Units represent the density of particles in the atmospheric column launched from Concordia standardized by the total number of particles.

ARs are known to have a major impact on the Antarctic surface climate, as they induce surface warming, either through enhanced sensible heat fluxes associated with warm air advection, or via the presence of mixed phase clouds with high supercooled liquid water and ice water content that increase downward longwave radiation towards the surface (Wille et al., 2019). In addition to triggering intense surface warming, AR also enable anomalously high inland moisture advection, making them a major contributor to extreme precipitation events over the desert polar plateau (Gorodetskaya et al., 2014; Turner et al., 2019; Wille et al., 2021; Adusumilli et al., 2021). The event of December 2018 recorded at Concordia station is in line with those characteristics, as Ricaud et al. (2020) highlighted the presence of a supercooled liquid cloud at 1500 m above Dome C and of precipitating ice crystals and/or blowing snow 20 December 2018. The impact of atmospheric rivers on water vapour isotopic anomalies has been studied for a specific event recorded in Greenland in July 2012 (Bonne et al., 2015), during which the large-scale transport of an air mass from Ivittuut to the NEEM drilling site caused a strong positive anomaly in temperature, specific humidity, and δD. However, large-scale transport only accounts for about half of the observed changes in water vapour isotopic composition. At Concordia Station, a +28 ‰ anomaly in water vapour δ¹⁸O was observed during the AR event of 15–19 March 2022, which originated from the Indian Ocean sector (Wille et al., 2024b). While large-scale dynamics appear to be a major driver of surface anomalies in specific humidity and δ¹⁸O, this study aims to better understand how local boundary-layer processes shape the surface signal during large-scale events.

In this study, we focus on the Dome C AR event of December 2018, with the aim of disentangling the impact of large-scale advection and boundary layer processes on the surface air vapour isotopic signal. In Sect. 2, we present our observational dataset, modelling tools and methodology. In Sect. 3, we describe in detail the December 2018 event, from large scale analysis to local observations at Dome C, and evaluate the LMDZ6iso model during this event. In Sect. 4, we study the influence of large-scale advection and boundary layer processes on the isotopic signal recorded in the water vapour and surface snow. Finally, we discuss our results in Sect. 5.

2 Materials and methods

2.1 Observations

2.1.1 Surface meteorology, longwave and shortwave radiative fluxes at the surface, radiosondes, and column-integrated atmospheric water

To place the event in a broader climatological context, we use hourly temperature and specific humidity from the Concordia AWS, provided by the IPEV/PNRA Project “Routine Meteorological Observation at Station Concordia” since 2005 (Grigioni et al., 2022). However, the AWS measurements are subject to two known limitations: (i) a warm daytime bias during summer under calm wind conditions, due to the lack of sensor ventilation (Genthon et al., 2011); and (ii) a dry bias at night, as standard thermohygrometers are unable to measure supersaturated conditions (Genthon et al., 2017). We remind the reader that, throughout this manuscript, the term night refers to an austral summer night, during which the sun remains above the horizon. For the detailed analysis of the event, we rely on measurements from the 45 m meteorological tower, located about 800 m upwind of the main station buildings, where the known AWS biases have been addressed through the CALVA instrumentation: temperature sensors are mechanically ventilated (Genthon et al., 2011), and updated thermohygrometers are specifically designed to detect supersaturation, even under extreme winter conditions (Genthon et al., 2017; Vignon et al., 2022). This dataset includes temperature, relative humidity, specific humidity, and wind speed measured approximately 3 m above the surface (Genthon et al., 2021a, b, 2022). Relative humidity is computed following the protocol described in Vignon et al. (2022) and Ollivier et al. (2025b). In addition, we use surface water vapour fluxes estimated by Ollivier et al. (2025b), based on the bulk method and derived from the same 3 m meteorological tower measurements, combined with snow surface temperature observations. For the bulk method, the stability functions from Holtslag and De Bruin (1988) are used for stable atmospheric conditions and from Högström (1996) for unstable atmospheric conditions. The roughness lengths for heat and moisture are assumed to be equal to the roughness length for momentum. Further details on the estimation of the surface water vapour flux can be found in Ollivier et al. (2025b). We also use upward and downward shortwave and longwave surface radiation from the Baseline Surface Radiation Network (BSRN) installed at Concordia Station (Lupi et al., 2021; Bai et al., 2022). Finally, we use available radiosondes, launched twice a day from Concordia Station throughout December 2018 (Grigioni et al., 2019). The soundings provide vertical profiles of temperature, relative humidity with respect to liquid water, and wind speed. Relative humidity is recalculated with respect to ice.

2.1.2 Isotopic composition of water vapour and surface snow

An infrared cavity ring-down spectrometer (Picarro L2130-i) was used to continuously measure the isotopic composition of water vapour, pumped from 1 m above the snow surface (Casado et al., 2016), and these data are available in Leroy-Dos Santos et al. (2021). This instrument was installed in 2018 at Concordia station, together with a home-made device to generate water vapour of known isotopic composition at low humidity to calibrate the analyser (Leroy-Dos Santos et al., 2021). Surface snow samples (0–1 cm) were collected every 2–4 h during the day in the clean area, approximately 800 m upwind of the main station buildings. The detailed sampling protocol and site description are provided in Ollivier et al. (2025b).

2.2 Models

2.2.1 Atmospheric river detection

The December 2018 AR was identified using an AR detection algorithm specifically designed for polar regions (Wille et al., 2019, 2021). In previous studies, this detection algorithm has been used to assess AR impacts such as surface melt, ice-shelf instability, and snowfall accumulation, and to characterise their synoptic properties (Wille et al., 2022; Pohl et al., 2021; Wille et al., 2024a, b; Gorodetskaya et al., 2023). Essentially, the AR detection algorithm used in this study searches for either integrated water vapour (IWV) or the meridional component of the integrated vapour transport (vIVT) between 37.5 and 85° S for values above the 98th percentile of all monthly climatological values defined from 1980–2021, using 3-hourly fields from ERA5 reanalysis vertically integrated on all levels. If this value forms a continuous segment that extends at least 20° in the meridional direction, then the segment is labelled as an AR. The December 2018 AR event was identified based on the patterns of IWV and vIVT. For more details on this algorithm, see the Methods in Wille et al. (2021).

2.2.2 The isotope-enabled general circulation model LMDZ6iso

We use the atmospheric general circulation model (AGCM) LMDZ6, described in Hourdin et al. (2020), using the version 20231022.trunk with the NPv6.3 physical package (Hourdin et al., 2023). This configuration is closely aligned with the atmospheric setup of IPSL-CM6A (Boucher et al., 2020), developed for phase 6 of the Coupled Model Intercomparison Project (CMIP6) (Eyring et al., 2016). The surface layer scheme and the turbulent diffusion parameterisation in LMDZ6 have been adapted to reproduce the structure and variability of the boundary layer at Dome C (Vignon et al., 2017b, 2018), including the summertime diurnal cycle and the two distinct stable regimes-weakly stable and very stable-that dominate during winter. We use LMDZ6's standard horizontal Low Resolution (LR) longitude-latitude grid (144 × 142), which corresponds to a 2.0° resolution in longitude and 1.67° in latitude. The vertical grid comprises 79 levels, with the lowest atmospheric level approximately 7 m a.g.l. at Dome C. Our simulation is nudged towards 6-hourly 3D fields of temperature and wind of the ERA5 reanalysis (Hersbach et al., 2020) using a relaxation time scale of 12 h. Nudging is excluded below the sigma-pressure level equivalent to 850 hPa a.s.l. (around 1000 m a.g.l. over the ice sheet), allowing the physics and dynamics of the model to operate freely within the boundary layer. Surface ocean boundary conditions are derived from ERA5 reanalysis monthly mean sea surface temperature and sea-ice concentration fields. A summary of the simulation evaluation is provided in Sect. 2.2.4.

2.2.3 Isotopic processes

The isotopic processes in LMDZ6iso, the isotope-enabled version of LMDZ6, are described in Risi et al. (2010), and a summary is provided here. In this study, we focus on the isotopes H${}_{\mathrm{2}}^{\mathrm{16}}$O, H${}_{\mathrm{2}}^{\mathrm{18}}$O, and HDO in water vapour, and use the standard δ notation to express the relative abundance of heavy stable water isotopes compared to the light isotope:

$$\begin{array}{}\text{(1)}& \mathit{\delta }=\left({\displaystyle \frac{R_{\mathrm{sample}}}{R_{\mathrm{VSMOW}}}}-\mathrm{1}\right)\cdot \mathrm{1000},\end{array}$$

where δ (expressed in ‰) refers to either δ¹⁸O or δD. Here, R_sample is the ratio of heavy to light water molecules in the sample (e.g., H${}_{\mathrm{2}}^{\mathrm{18}}$O $/$ H${}_{\mathrm{2}}^{\mathrm{16}}$O or HDO $/$ H${}_{\mathrm{2}}^{\mathrm{16}}$O), and R_VSMOW is the equivalent ratio in the Vienna Standard Mean Ocean Water (VSMOW), used as reference. For moisture transport, LMDZ6 uses the Van Leer advection scheme (Van Leer, 1977; Risi et al., 2010). Under first-order fractionation behaviour, variations in δD are approximately eight times those in δ¹⁸O (Craig, 1961). Deviations from this relationship are quantified using the deuterium excess (d-excess):

$$\begin{array}{}\text{(2)}& \text{d-excess}=\mathit{\delta }\text{D}-\mathrm{8}\cdot {\mathit{\delta }}^{\mathrm{18}}\text{O}.\end{array}$$

The d-excess is particularly sensitive to kinetic fractionation processes (Dansgaard, 1964), and thus serves as a useful tracer of evaporation conditions at the moisture source (Gat et al., 2011), as well as non-equilibrium processes such as ice condensation under supersaturated conditions or evaporation of precipitating droplets (Jouzel and Merlivat, 1984). Equilibrium fractionation coefficients between water vapour and liquid or ice are derived from Merlivat and Nief (1967) and Majoube (1971). Kinetic fractionation effects are parametrized for sea surface evaporation following Merlivat and Jouzel (1979) and for ice condensation and supersaturation based on Jouzel and Merlivat (1984), where the kinetic fractionation coefficient α_cin is defined as:

$$\begin{array}{}\text{(3)}& {\mathit{\alpha }}_{\mathrm{cin}}={\displaystyle \frac{S_{\mathrm{i}}}{\mathrm{1}+{\mathit{\alpha }}_{\mathrm{eq}}\cdot \frac{D}{D_{\mathrm{iso}}}\cdot (S_{\mathrm{i}}-\mathrm{1})}},\end{array}$$

where S_i is the supersaturation during ice condensation (relative humidity with respect to ice, dimensionless), α_eq is the equilibrium fractionation coefficient, and D and D_iso are the diffusivities of most abundant water and the isotope respectively (in m² s⁻¹). In LMDZ6iso, supersaturation used for isotope processes calculation is parametrized as a simple linear function of temperature for freezing conditions and is used exclusively to compute the fractionation coefficient:

$$\begin{array}{}\text{(4)}& S_{\mathrm{i}}=\mathrm{1}-\mathit{\lambda }\cdot \left(T-T_{\mathrm{freeze}}\right),\end{array}$$

where λ is a tuning parameter ranging from 0 K⁻¹ (no supersaturation with respect to ice) to 0.006 K⁻¹ (high supersaturation with respect to ice), T is the air temperature (K), and T_freeze is the freezing temperature (273.15 K). To reduce the bias in surface snow δ¹⁸O and d-excess, we performed our simulation using a value of λ = 0.004 K⁻¹, as described in Dutrievoz et al. (2025b).

In the model, the isotopic composition of surface snow is modelled using a snow bucket representing the average isotopic composition of snowfall since the start of the simulation. The isotopic signature of surface snow can only be modified by fractionation during surface vapour condensation, as sublimation is assumed to occur without fractionation, consistent with common isotope-enabled AGCM assumptions (e.g., Hoffmann et al., 1998). Recent studies have nevertheless shown that sublimation and vapour exchange can significantly influence the isotopic composition of surface snow, but it is not yet taken into account in LMDZ6iso (Casado et al., 2021; Hughes et al., 2021; Harris Stuart et al., 2023; Dietrich et al., 2023).

2.2.4 Simulation evaluation

The LMDZ6iso simulation used in this study is similar to the one described and evaluated in Dutrievoz et al. (2025b), but with a ERA5 nudging relaxation time scale of 12 h instead of 3 h, which does not significantly affect the results over Antarctica. A summary of the simulation evaluation is provided here. Over Antarctica for the period 1979–2024, LMDZ6iso captures well the spatial distribution of mean surface temperature with a mean cold bias of 1.4 K. Simulated snow accumulation exhibits a mean bias of +12.2 kg m⁻² yr⁻¹, representing 8.7 % of the mean observed accumulation, independently of the altitude. When ran with an intermediate supersaturation parameter (λ = 0.004 K⁻¹), the model accurately reproduces the spatial distribution of the mean annual surface snow δ¹⁸O across the continent with no bias when considering the entire ice sheet. However, a positive bias is observed inland, where measured δ¹⁸O values are lower than −45 ‰, with a bias of +2.6 ‰ at Dome C, while a negative bias appears in coastal regions, where δ¹⁸O exceeds −40 ‰. Regarding surface snow d-excess, LMDZ6iso systematically overestimates d-excess compared to observations (+6.6 ‰ at Concordia), except at the ice sheet margins, with a maximum positive bias of 11.9 ‰ inland. At Concordia, LMDZ6iso accurately simulates temperature and specific humidity during clear-sky diurnal cycles in December 2018, despite a slight positive bias (1.5 K and 0.04 g kg⁻¹, respectively). The mean water vapour δ¹⁸O over the same period is −66.9 ‰ in the observations and −62.2 ‰ in LMDZ6iso. However, the modelled amplitude of vapour δ¹⁸O is more than twice as large as observed (4.3 ‰ observed vs. 11.5 ‰ modelled). This discrepancy likely reflects an overly strong enrichment during the day and/or an excessive depletion at night.

2.3 Water vapour isotopic budget using tendencies

2.3.1 Computation of isotopic tendencies

To investigate the processes governing the atmospheric water budget and the isotopic composition of water vapour at the surface, we use the LMDZ6iso water vapour tendencies. These tendencies represent the rate of change of each water vapour isotope (in kg kg⁻¹ s⁻¹) due to specific processes driving the water cycle: large-scale advection, cloud and precipitation condensation and sublimation, shallow convection, deep convection (negligible in Antarctica), and vertical turbulent diffusion (which encompasses surface sublimation and condensation).

The total atmospheric water vapour budget and vapour isotopic budget are computed as the sum of contributions from all processes:

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{d}q}{\mathrm{d}t}}=\sum _{\mathrm{process}}{\left.{\displaystyle \frac{\mathrm{d}q}{\mathrm{d}t}}\right|}_{\mathrm{process}},\text{(6)}& {\displaystyle \frac{\mathrm{d}\mathit{\delta }}{\mathrm{d}t}}=\sum _{\mathrm{process}}{\left.{\displaystyle \frac{\mathrm{d}\mathit{\delta }}{\mathrm{d}t}}\right|}_{\mathrm{process}},\end{array}$$

where $\mathrm{d}q/\mathrm{d}t$ is the change in time of atmospheric specific humidity (in kg kg⁻¹ s⁻¹) and $\mathrm{d}\mathit{\delta }/\mathrm{d}t$ is the change in time of the δ of atmospheric vapour (in ‰ s⁻¹).

The rate of change of δ in atmospheric vapour attributable to each process, $\mathrm{d}\mathit{\delta }/\mathrm{d}t{|}_{\mathrm{process}}$ (in ‰ s⁻¹), is calculated as follows:

$$\begin{array}{}\text{(7)}& {\left.{\displaystyle \frac{\mathrm{d}\mathit{\delta }}{\mathrm{d}t}}\right|}_{\mathrm{process}}={\displaystyle \frac{\mathrm{1}}{q}}\cdot {\left.{\displaystyle \frac{\mathrm{d}q}{\mathrm{d}t}}\right|}_{\mathrm{process}}\cdot \left({\mathit{\delta }}_{\mathrm{process}}-\mathit{\delta }\right),\end{array}$$

where q represents the atmospheric specific humidity before the process (in kg kg⁻¹), δ is the atmospheric vapour δ¹⁸O before the process (in ‰), ${\left.\mathrm{d}q/\mathrm{d}t\right|}_{\mathrm{process}}$ is the specific humidity tendency associated to the process (in kg kg⁻¹ s⁻¹), and δ_process is the isotopic composition of this specific humidity tendency at each time step (in ‰). A detailed derivation of this equation is provided in Dutrievoz et al. (2025b), Appendix A.

2.3.2 Quantification of the origin of water and δ¹⁸O in the near-surface atmospheric layer and in the mixed layer

We now aim to identify the origin of the atmospheric water sources. The associated methodology is illustrated in Fig. 2. This study focuses on the analysis of δ¹⁸O, but a similar analysis using δD would yield equivalent results. Based on the tendencies, we intend to decompose the total specific humidity into a sum of contributions from different processes. We define an initial vapour bucket composed of the humidity present in the atmosphere, q_init and $q_{\mathrm{init}}^{\mathrm{iso}}$, at the initial time step t₀. Then at the next time step t + Δt, positive contributions ($\mathrm{d}q/\mathrm{d}t{|}_{\mathrm{process}}\cdot \mathrm{\Delta }t$ > 0 or $\mathrm{d}{q}^{\mathrm{iso}}/\mathrm{d}t{|}_{\mathrm{process}}\cdot \mathrm{\Delta }t$ > 0, Eq. 9) are added to their respective buckets q|_process, which quantify the amount of water originating from each source process (e.g., surface sublimation, cloud or precipitation sublimation, or water advected by the atmospheric dynamics). It leads to intermediate quantities $q^{+}{|}_{\mathrm{process}}(t+\mathrm{\Delta }t)$, the sum of which is larger than q(t+Δt) (Eq. 10). Negative contributions (either $\mathrm{d}q/\mathrm{d}t{|}_{\mathrm{process}}\cdot \mathrm{\Delta }t$ < 0 or $\mathrm{d}{q}^{\mathrm{iso}}/\mathrm{d}t{|}_{\mathrm{process}}\cdot \mathrm{\Delta }t$ < 0) correspond to mass sinks for all vapour buckets (e.g., surface condensation, cloud condensation, water advected by the atmospheric dynamics or shallow convection), so the sum of negative tendencies is subtracted proportionally to each bucket (as illustrated by the rescaling in Fig. 2), resulting in the sum of buckets being equal to the specific humidity q(t+Δt) (Eq. 11). By applying this scheme iteratively at each time step, the initial bucket is progressively depleted-typically within a few days-and replaced by a sum of buckets corresponding to processes with positive contributions to the specific humidity. In other words, when a process adds water to the atmosphere, the origin of the water can be traced back based on the nature of the contributing process. In contrast, when water is removed from the atmosphere, it originates from the mixed vapour from different origins; therefore, the loss is distributed proportionally among all existing buckets. We chose to first sum the positive contributions and then subtract the negative contributions at each time step. This order of operations is not symmetric: applying the negative tendencies first and then adding the positive ones leads to similar but slightly different results. In the following, $\mathrm{d}q{|}_{\mathrm{process}}\left(t\right)=\mathrm{d}q/\mathrm{d}t{|}_{\mathrm{process}}\cdot \mathrm{\Delta }t$, and $\mathrm{d}{q}^{+}{|}_{\mathrm{process}}\left(t\right)=\mathrm{d}{q}^{+}/\mathrm{d}t{|}_{\mathrm{process}}\cdot \mathrm{\Delta }t$:

$$\begin{array}{}\text{(8)}& {\displaystyle }q\left(t\right)=\sum _{\mathrm{process}}q{|}_{\mathrm{process}}\left(t\right)\text{(9)}& {\displaystyle }\mathrm{d}{q}^{+}{|}_{\mathrm{process}}\left(t\right)=\left\{\begin{array}{ll}\mathrm{0}& \text{if}\mathrm{d}q{|}_{\mathrm{process}}\left(t\right)<\mathrm{0},\\ \mathrm{d}q{|}_{\mathrm{process}}\left(t\right)& \text{else},\end{array}\right.\text{(10)}& {\displaystyle }{q}^{+}(t+\mathrm{\Delta }t)=\sum _{\mathrm{process}}\underset{q^{+}(t+\mathrm{\Delta }t){|}_{\mathrm{process}}}{\underbrace{q{|}_{\mathrm{process}}\left(t\right)+{\left.\mathrm{d}{q}^{+}\right|}_{\mathrm{process}}\left(t\right)}},\text{(11)}& {\displaystyle }q(t+\mathrm{\Delta }t){|}_{\mathrm{process}}=q^{+}(t+\mathrm{\Delta }t){|}_{\mathrm{process}}\cdot {\displaystyle \frac{q(t+\mathrm{\Delta }t)}{q^{+}(t+\mathrm{\Delta }t)}},\end{array}$$

where $\mathrm{d}{q}^{+}{|}_{\mathrm{process}}\left(t\right)$ corresponds to the positive contributions of each process to the total and isotopic specific humidity at time step t (in kg kg⁻¹ s⁻¹), $q^{+}(t+\mathrm{\Delta }t)$ refers to the specific humidity from the previous time step incremented by the sum of all positive contributions between t and t+Δt (in kg kg⁻¹) and q(t+Δt) is the specific humidity at time t+Δt (in kg kg⁻¹).

Figure 2 Schematic representation of the methodology used to determine the origin of atmospheric water sources. Depending on whether their contribution is positive or negative, the processes are classified either as sources ($\mathrm{d}q/\mathrm{d}t$ > 0), shown with solid colours, or as sinks ($\mathrm{d}q/\mathrm{d}t$ < 0), shown with hatched colours. Each bucket is supplied by a specific water source and is represented by a distinct colour: the initial bucket, consisting of the water already present in the atmosphere, is shown in grey; the orange bucket corresponds to water originating from surface sublimation; the blue bucket represents water transported by advection; and the purple bucket accounts for water from cloud and precipitation sublimation. In this example, shallow convection, represented in red, does not correspond to a bucket, as it only removes water from the atmosphere. Paler colours indicate newly added water contributions.

Using Eq. (11) for total water and their equivalents for q^iso, we compute the global δ and its anomaly relative to the clear-sky mean diurnal cycle value, as a function of the isotopic contributions associated with each process:

$$\begin{array}{}\text{(12)}& {\displaystyle }\mathit{\delta }{|}_{\mathrm{process}}=\left({\displaystyle \frac{q^{\mathrm{iso}}{|}_{\mathrm{process}}/q{|}_{\mathrm{process}}}{R_{\mathrm{VSMOW}}}}-\mathrm{1}\right)\cdot \mathrm{1000},\text{(13)}& {\displaystyle }\mathit{\delta }=\sum _{\mathrm{process}}\left[{\displaystyle \frac{q{|}_{\mathrm{process}}}{q}}\cdot \mathit{\delta }{|}_{\mathrm{process}}\right],\text{(14)}& {\displaystyle }\mathit{\delta }-{\mathit{\delta }}_{\mathrm{mean}}=\sum _{\mathrm{process}}\left[{\displaystyle \frac{q{|}_{\mathrm{process}}}{q}}\cdot \left(\mathit{\delta }{|}_{\mathrm{process}}-{\mathit{\delta }}_{\mathrm{mean}}\right)\right],\end{array}$$

where δ is the isotopic composition of atmospheric vapour (in ‰), δ|_process is the isotopic composition of atmospheric vapour associated to each process at each time step (in ‰) and δ_mean is the isotopic clear-sky mean diurnal cycle value (in ‰). The derivation of Eqs. (13) and (14) are detailed in Appendix A.

We also performed a water budget analysis within the mixed layer to ensure consistency with the surface signal. The mixed layer is defined as the surface atmospheric layer (7 m a.g.l. in the model) and the layers mixed by shallow convection. During the night, in the absence of convection and under weak turbulent conditions, a decoupling occurs between the surface and the rest of the atmosphere. During the day, shallow convection mixes the air between the surface and approximately 400–500 m a.g.l. The growth of the mixed layer mixes the air from above into the mixed layer, introducing a new process which we refer to as entrainment. The entrainment term, illustrated in Fig. S1 in the Supplement, includes water originating from various atmospheric processes. Therefore, this additional process must be considered when performing a water and isotope budget in a mixed layer with a varying altitude.

3 Description of the event and model evaluation

3.1 Large scale conditions during the event

Although Dome C is geographically closer to the Indian sector of the Southern Ocean – and most events affecting the station generally originate from this sector (Sodemann and Stohl, 2009; Genthon et al., 2016; Schlosser et al., 2017) – analysis of the synoptic-scale atmospheric circulation associated with this event indicates a much more distant moisture source, originating from the Atlantic sector, as shown in Fig. 1. The AR detection algorithm, presented in Sect. 2.2.1, depicts an initial landfall over western Dronning Maud Land on 17 December, with the AR progressing across the Antarctic Plateau. The moisture intrusion reaches Dome C on 19 December, but is not classified as an AR since the AR detection algorithm domain only extends until 85° S and cannot track moisture transport that passes over the South Pole moving from south to north (Fig. 3). These warm air masses, characterised by high integrated water vapour over Dome C, contain large amounts of cloud ice and liquid water, which contribute to surface warming through enhanced sensible heat fluxes and increased downward longwave radiation (Ricaud et al., 2024a). After 22 December, the moisture quantities decreases on the Antarctic Plateau, but temperatures remain elevated around Dome C due to residual moisture in the region. MODIS satellite imagery confirms this timeline, showing the initial inland cloud penetration on 18 December (Fig. 3b) and then traversing the Antarctic Plateau towards Dome C by the 21th (Fig. 3d).

Figure 3 Atmospheric river on (a, b) 18 December 2018 at 00:00 UTC, before the AR crossed the Antarctic Plateau, (c) 21 December 2018 at 06:00 UTC and (d) 21 December 2018 at 00:00 UTC, after the AR had reached Dome C. (a, c) Shape of the detected AR, with anomalies in integrated water vapour transport (IWV) relative to the December climatology (1980–2021) shown in shading. Red and orange contours correspond to ARs detected using the vIVT and IWV based detection methods, respectively. (b, d) MODIS satellite imagery illustrating the progression of the AR across the Antarctic Plateau. Satellite images were obtained from the NASA MODIS instrument via the NASA Worldview application (https://worldview.earthdata.nasa.gov, last access: 2 February 2026)).

This seemingly unusual path of moisture transport over the Antarctic Plateau and the South Pole results from a pronounced atmospheric ridge that extends across the entire continent. The ridge directs the moisture flow towards the coastline during the initial AR landfall on the 18th and then expands across the Antarctic Plateau, reaching the other side of the continent (Fig. S2a). The latent heat release associated with the poleward moisture transport likely contributes to the development of the ridge over the Antarctic continent as evidenced by the prolonged period of positive potential vorticity anomalies around the Weddell Sea extending back to 17 December (Figs. S2 and S3). In addition to the mid-level 500 hPa geopotential height anomalies encompassing most of the Antarctic Plateau, anomalies are also present in the stratosphere at 10 hPa (Fig. S2b). These height anomalies in the stratosphere, observed for most of December and extending downward to the surface during the AR event (Fig. S4a), may indicate a potential weakening of the Southern polar vortex, as persistent lower-level positive temperature anomalies exceeding 5 °C are observed across the continent for most of December, peaking during the major moisture intrusion event from 17 to 22 December (Fig. S4b).

Now that we have characterized the AR, we examine the spatial anomalies in temperature, specific humidity, relative humidity, δ¹⁸O, and d-excess during the event (12:00 UTC on 20 December 2018), relative to the December average (1980–2021) at the sigma-pressure level corresponding to 850 hPa geopotential height in LMDZ, equivalent to an altitude of approximately 1500 m a.s.l. and about 1000 m above the Antarctic Plateau (Fig. S5). The AR pathway exhibits a pronounced temperature anomaly exceeding +10 °C, extending from the oceanic boundary near Dronning Maud Land (Fig. 1) to the edge of the Antarctic Plateau near Adélie Land (Fig. 4a). The AR is also associated with a significant positive anomaly in specific humidity, although spatially more confined compared to the temperature anomaly (Fig. 4b). East of the AR, over Dronning Maud Land and within the 5 °C anomaly contour, the region shows a strong negative anomaly in relative humidity (up to −60 %), with completely dry air in some locations (below 10 %, Figs. 4c and S6c). In the same region, the AR is characterised by a strong negative δ¹⁸O anomaly of −20 ‰ reaching −80 ‰ (Figs. 4d, S6d), and the d-excess displays a strong positive anomaly of +40 ‰, with absolute values reaching up to +80 ‰ (Fig. 4e, S6e). The AR is associated with a positive δ¹⁸O anomaly of approximately +10 ‰, spatially distributed over the same region as the specific humidity anomaly (Fig. 4b, d).

Figure 4 LMDZ6iso anomalies at the model level equivalent to 850 hPa a.s.l. of (a) temperature (°C), (b) specific humidity (g kg⁻¹), (c) relative humidity (%), (d) δ¹⁸O (‰), and (e) d-excess (‰), calculated on 20 December 2018 at 12:00 UTC relative to the December average (1980–2021). The purple contour represents the 5 °C anomaly boundary, while the yellow contour indicates the 10 °C anomaly boundary. The black dot indicates Concordia station. The altitude above the surface corresponding to the LMDZ6iso model level equivalent to 850 hPa is shown in Fig. S5 (around 1000 m a.g.l. over the ice sheet).

3.2 Anomalies monitored at the surface

The AR event at Dome C is identified from the 45 m meteorological tower observations (Sect. 2.1.1) as a period of anomalously high temperature and specific humidity relative to the mean December diurnal cycle, beginning with an increase on 19 December and ending with a return to near-average values on 23 December.

Over the 13-year AWS record period (2006–2018), the December 2018 event is characterized by the second highest maximum temperature (−14.7 °C) and the absolute record for specific humidity (1.54 g kg⁻¹) throughout the 1-hourly summer (December–January) measurements (Fig. 5). Consequently, this event stands out for its exceptionally high temperature and specific humidity compared to the previous 13 summers.

Figure 5 Hourly temperatures (a) and specific humidity (b) recorded at Dome C during summers (December–January) from December 2006 to January 2024 from the AWS. The year of interest, 2018 (December 2018–January 2019), is coloured in dark blue. Each box is delimited by the 0.25 and 0.75 percentiles and the orange bar represents the median. Black bars indicate the values inside 1.5 times the inter-quartile range. Red circles indicate values outside the 0.999 percentiles. The dotted red line indicate the value reached during the December 2018 AR.

The water vapour flux from the surface to the lower atmosphere in December 2018 is estimated using meteorological measurements at Dome C and the bulk method described in Sect. 2.1.1. The roughness length for momentum (z₀) used in the reference calculation is set to 1 mm, corresponding to the value in LMDZ. Sensitivity tests were also performed using two alternative values, 0.01 and 6.3 mm, as estimated in Vignon et al. (2017a). The resulting water vapour flux is found to be highly sensitive to the choice of z₀, with larger values of z₀ leading to increased fluxes (Fig. S7). The water vapour flux is characterised by pronounced daytime surface sublimation, while nocturnal conditions are associated with condensation. During the event, the model underestimates the first sublimation peak (21 December) by 74 % (0.046 kg m⁻² h⁻¹ observed, 0.012 kg m⁻² h⁻¹ modelled), and the second peak (22 December) by 44 % (0.048 kg m⁻² h⁻¹ observed, 0.027 kg m⁻² h⁻¹ modelled). On the night of 21 December, the model simulates night-time condensation that is nearly two times stronger than observed (−0.006 kg m⁻² h⁻¹ observed, −0.011 kg m⁻² h⁻¹ modelled).

To assess the intensity of the observed anomaly and evaluate the performance of the LMDZ6iso model, we focus on analysing surface meteorological data recorded by the 45 m meteorological tower (presented in Sect. 2.1.1) and δ¹⁸O in the vapour recorded by the Picarro instrument (detailed in Sect. 2.1.2). We compare the AR event (19–23 December) to the mean diurnal cycle calculated throughout the month of December 2018, excluding the event. The entire period is shown in Fig. S8.

On 21 December 2018 at 07:00 UTC, the temperature and humidity measurements are −14.6 °C and 1.45 g kg⁻¹ respectively (Fig. 6), in line with records of the AWS station. It exceeds the December 2018 clear sky daily cycle averages by 17.9 °C and 1.16 g kg⁻¹, respectively (Fig. 6a, d). The event is characterized by two pronounced peaks in incoming longwave radiation (Fig. 6b) associated with the presence of mixed phase clouds with high supercooled liquid water and ice water content (Ricaud et al., 2024b). Additionally, the AR event is marked by a nocturnal saturation of relative humidity with no episodes of supersaturation except on 23 December, alongside a doubling of wind speed (Fig. 6c, e). The water vapour δ¹⁸O exhibits a pronounced enrichment of 17.2 ‰ relative to the clear-sky diurnal cycle average for December 2018, with a maximum value of −49.6 ‰ recorded on 22 December at 01:00 UTC, reaching values comparable to those observed in surface snow (Fig. 6f). During the event, despite significant spatial variability (more the 8 ‰), surface snow exhibits a mean δ¹⁸O value of −53.6 ‰, and becomes enriched by 2.3 ‰ after the event, reaching −51.3 ‰ (Fig. 6f).

Figure 6 Evolution of (a) air temperature (°C) from the first level of the Concordia meteorological tower (2.8 m a.g.l.), (b) incoming longwave flux from the Baseline Surface Radiation Network, (c) relative humidity with respect to ice (%), (d) specific humidity (g kg⁻¹), (e) wind speed (m s⁻¹) from the first level of the Concordia meteorological tower (3.0 m a.g.l.) and (f) δ¹⁸O (‰ from the Concordia Picarro instrument (1.0 m a.g.l.). In panel (f), observed surface snow is shown as black dots, while modelled surface snow is represented by a dotted line. The grey rectangle indicates the period of the AR. The blue line represents the model output (first level: 6.7 m a.g.l.), while the black line corresponds to the observations. The grey and blue dashed lines represent the mean diurnal cycle calculated over the entire month of December 2018, excluding the event (19–23 December), for observations and the model, respectively. The orange shading indicates the difference between the mean clear-sky diurnal cycle and the observations during the event, while the red shading highlights the peaks in downward longwave radiation.

During the event, the model captures the two main peaks in downward longwave radiation, albeit with underestimated peak intensities (22 % lower for the first maximum peak and 34 % for the second maximum peak) (Fig. 6b). Outside the event, the model reproduces the diurnal temperature cycle with an amplitude 11 % larger than observed (9.3 °C observed, 10.2 °C modelled), associated with a mean warm bias of 1.6 °C, and simulates the specific humidity cycle with an amplitude 14 % higher than observed (0.24 g kg⁻¹ observed, 0.27 g kg⁻¹ modelled), with a mean moist bias of 0.05 g kg⁻¹. During the event, the maximum temperature bias is −0.8 °C, while the maximum specific humidity bias is −0.01 g kg⁻¹ (1 % lower than observed) (Fig. 6a, d). Outside the event, the model fails to reproduce the nocturnal period of supersaturation. However, during the event, when no supersaturation is observed, the model accurately reproduces the diurnal cycle of relative humidity (Fig. 6c). The LMDZ6iso model also captures variations in wind speed during the event, with a slight mean bias of 0.8 m s⁻¹ (11 % lower than observed) (Fig. 6e). While the model overestimates the amplitude of the diurnal cycle of vapour δ¹⁸O by more than a factor of two (4.3 ‰ observed, 11.5 ‰ modelled), as detailed in Dutrievoz et al. (2025b) and Ollivier et al. (2025a), it successfully captures vapour δ¹⁸O during the event's peak, with a positive bias of 1.8 ‰ (Fig. 6f). The modelled surface snow δ¹⁸O remains constant at a value of −48 ‰, which is 5 ‰ higher than observed during the event (Fig. 6f).

3.3 Vertical atmospheric profile during the event

We examine the vertical structure of the atmosphere during the AR event and during the unperturbed period outside the event (1–18 and 24–31 December 2018). We recall that the LMDZ6iso model is nudged to ERA5 reanalysis for temperature and wind above the boundary layer, and it reproduces these fields consistently with the reanalysis (not shown). In December 2018, outside of the event, the mean temperature increases from −31.5 °C at the surface to −30.1 °C at 89 m a.g.l., before decreasing to −38.1 °C at 2032 m a.g.l. (Fig. 7a). During the AR (21 December 2018 at 11:00 UTC), the temperature is 12.3 °C higher than under non-AR conditions within the first 2000 m a.g.l. Outside the AR period, specific humidity averages 0.3 g kg⁻¹ throughout the atmospheric column up to 2000 m a.g.l., whereas during the event it is on average more than three times higher and decreases progressively from 1.3 g kg⁻¹ at the surface to 0.7 g kg⁻¹ at 1964 m a.g.l. (Fig. 7b). Relative humidity with respect to ice is close to saturation both before and during the event (Fig. 7c). Wind speed outside the AR period averages 3.1 m s⁻¹ throughout the atmospheric column up to 2000 m a.g.l. During the AR, the wind speed increases sharply from 6.7 m s⁻¹ at the surface to 25.9 m s⁻¹ at 1582 m a.g.l. (Fig. 7d). The model accurately captures variations in temperature, specific humidity, and wind speed both during and outside the AR event, but it underestimates relative humidity over the lowest 2000 m a.g.l. by 9 % during the event and by 16 % outside the event. Vertical profiles above 2000 m a.g.l. are shown in the Supplement (Fig. S9).

Figure 7 Vertical profiles of (a) temperature (°C), (b) specific humidity (g kg⁻¹), (c) relative humidity with respect to ice (%) and (d) wind speed (m s⁻¹) from the surface to 2000 m a.g.l., based on average radiosonde measurements at the Concordia station during December 2018 excluding the event (19–23 December) (twice per day, black solid lines) and during the AR event (21 December 2018 at 11:00 UTC, black dashed lines). Corresponding LMDZ6iso model output at the nearest grid point for the same times are shown in blue lines (December averages) and red lines (AR event). The standard deviation is shown in gray for observations and in blue for LMDZ6iso for 52 radiosondes during December 2018 outside the AR event.

4 Attribution of the isotopes anomalies at the surface

We now use the model to obtain a more detailed understanding of the sources of isotopic variability measured at the surface during the AR event. In the model, variations in specific humidity and water isotopes are governed by the four processes outlined in Sect. 2.3: surfaces fluxes, cloud condensation and precipitation sublimation, shallow convection and advection. Our goal is to identify the contributions of these processes to the observed surface variations in specific humidity and δ¹⁸O between 19 and 23 December 2018. At Concordia, the diurnal cycle is characterised by maximum solar insolation at 03:00 UTC and minimum insolation at 15:00  UTC (local time at Concordia is UTC+8). The warm temperature anomaly characteristic of the AR event at the surface extends up to 2000 m a.g.l. throughout the study period (Fig. 7a). This anomaly overlays the diurnal cycle, which is marked by nocturnal surface cooling (Fig. 8a). For specific humidity, a positive anomaly associated with the AR is observed from 20–21 December within the first 1000 m a.g.l., with a peak intensity on 21 December (Fig. 8b). Relative humidity exhibits a sharp decrease to approximately 20 % from 18 December, evident above 200 m a.g.l. This is followed by a period of saturation starting on 20 December and lasting until early 22 December, before transitioning back to a drier phase (relative humidity lower than 20 %) on 22 and 23 December, above 300 m a.g.l. (Fig. 8c).

Figure 8 Modelled vertical profiles of (a) temperature (°C), (b) specific humidity (g kg⁻¹), and (c) relative humidity with respect to ice (%) at Concordia. The red outline indicates regions where the relative humidity is higher than 98 %. The dotted lines indicate the height of the mixed layer during shallow convection. The dashed vertical lines delineate the boundaries of the atmospheric river. All times are given in UTC.

Based on the model's water vapour tendencies, daytime surface sublimation and the turbulent transport of moisture in the first 200 m a.g.l., extending up to 400 m a.g.l. on 21 and 22 December, result in an enrichment of atmospheric vapour in δ¹⁸O (Fig. 9a, b). During nighttime, the temperature decreases, leading to saturation and condensation near the surface, which causes a reduction in δ¹⁸O in surface vapour. The analysis of specific humidity and δ¹⁸O tendencies in the surface layer (0–6.7 m a.g.l.) confirms the strong influence of surface sublimation and turbulence on the surface signal, both during the typical diurnal cycle and the AR event (Fig. 10). Shallow convection mixes this humidity within the mixed layer (from the surface to 200–400 m a.g.l.) during the day, transporting surface moisture upward. This process decreases δ¹⁸O at the surface while enriching vapour at higher altitudes (Figs. 9c, d and 10). Cloud condensation leads to a depletion of both moisture and δ¹⁸O around 400 m during the peak of the event, as well as during the post-event diurnal cycle (Fig. 9e, f). In contrast, the sublimation of clouds and precipitation enriches the air in both moisture and δ¹⁸O bellow the cloud, around 200 m a.g.l. Finally, advection transports large-scale moisture enriched in δ¹⁸O to the surface starting on 21 December, followed by a decrease in both moisture and δ¹⁸O on 22 December (Figs. 9g, h and 10).

Figure 9 Decomposition of the rate of change in time of the different processes controlling specific humidity (a, c, e, g: in g kg⁻¹ h⁻¹) and vapour δ¹⁸O (b, d, f, h: in ‰ h⁻¹) as a function of time at Concordia station. Panels (a) and (b) correspond to surface sublimation and turbulent mixing, panels (c) and (d) to shallow convection, panels (e) and (f) to cloud condensation and precipitation sublimation and panels (g) and (h) to advection. It should be noted that the sign of changes is relative to water vapour so the surface sublimation is positive (increase in water vapour) and the cloud condensation is negative (decrease in water vapour). The dotted lines indicate the height of the mixed layer during shallow convection. The dashed vertical lines delineate the boundaries of the atmospheric river. All times are given in UTC.

Figure 10 Rate of temporal change in the surface atmospheric layer due to the different processes contributing to the (a) specific humidity tendency (g kg⁻¹ h⁻¹) and (b) vapour δ¹⁸O tendency (‰ h⁻¹) at Concordia station. The black line represents the signal obtained from the first level model output. The individual contributions to the rate of change of humidity and δ¹⁸O are shown in colour: orange for vertical turbulent diffusion and surface fluxes, blue for advection, purple for cloud condensation and precipitation sublimation and red for shallow convection. Processes with positive contributions are shown in solid colours, whereas processes with negative contributions are shown using hatched shading. It should be noted that the sign of changes is relative to water vapour so the surface sublimation is positive (increase in water vapour) and the cloud condensation is negative (decrease in water vapour). The dashed vertical lines delineate the boundaries of the atmospheric river.

We now aim to quantify the relative contributions from large-scale processes (advection) and local processes (surface and atmospheric sublimation) to the specific humidity and δ¹⁸O signals during the AR. For this, we track the sources of atmospheric vapour through processes with positive tendencies. We consider that atmospheric vapour is the mixing of 3 buckets corresponding to 3 different sources: surface sublimation, large-scale advection, and cloud and precipitation sublimation (Fig. 11a, b), following the approach described in Sect. 2.3. As shallow convection always transports water from the surface towards the upper layers (shallow convection tendencies always negative in the surface atmospheric layer, Fig. 10), it does not contribute to the surface specific humidity as a water source. At the peak of the event, the large-scale advection accounts for 30 % of the moisture amount (0.44 g kg⁻¹), while surface sublimation accounts for 69 % (1.01 g kg⁻¹). Cloud and precipitation sublimation plays a minor role, contributing for 1 % of the surface vapour (0.02 g kg⁻¹). To facilitate the interpretation of the total isotopic signal, we analysed the δ¹⁸O anomaly relative to the mean vapour δ¹⁸O during the clear-sky diurnal cycles (−62 ‰). The total δ¹⁸O anomaly (Fig. 11d) is computed as the weighted sum of each process’s δ¹⁸O anomaly value (Fig. 11c), scaled by its moisture contribution (Fig. 11b), according to Eqs. (13) and (14). The δ¹⁸O associated with surface sublimation exhibits a pronounced diurnal cycle, closely matching the total δ¹⁸O, with a maximum of −51 ‰ at the peak of the event (Fig. 11c). The δ¹⁸O related to large-scale advection and cloud and precipitation sublimation shows a significant depletion on 20 December (−94 ‰ and −85 ‰, respectively), followed by a sharp increase, reaching maximum values of −51 ‰ and −40 ‰, respectively, at the event's peak. During the peak of the event, surface sublimation is the primary contributor of surface vapour increase in δ¹⁸O, enriching the mean surface vapour by 7.8 ‰ compared to the mean clear-sky diurnal cycle, which accounts for 68 % of the total increase. Large-scale advection contributes for an additional 3.4 ‰, corresponding to 30 % of the increase. Finally, cloud and precipitation sublimation slightly enriches the vapour by 0.2 ‰, contributing to 2 % of the overall increase. The complete period corresponding to Fig. 11a and d is presented in Fig. S10.

Figure 11 (a) Temporal evolution of the positive contributions of different processes to the surface specific humidity (6.7 m a.g.l.) and (b) their respective percentages. (c) Evolution of the δ¹⁸O anomaly relative to the mean clear-sky diurnal cycle value (δ¹⁸O_mean = −62 ‰), associated with different processes. (d) Positive contributions of different processes to the surface δ¹⁸O anomaly relative to this mean diurnal value. The black line represents the signal obtained from total specific humidity and δ¹⁸O. Individual contributions to specific humidity and δ¹⁸O are shown in colour: orange for surface sublimation and vertical turbulent diffusion, blue for advection, and purple for cloud and precipitation sublimation. Vertical dashed lines mark the AR period (19–23 December), while the shaded area highlights the peak of the event, reached on 21 December at 09:30 UTC.

5 Discussion and conclusion

The December 2018 AR event observed at Concordia, which originated from the Atlantic sector and crossed the Antarctic continent, illustrates the interaction between local processes and large-scale mechanisms that control the isotopic signature of surface water vapour. This is the first atmospheric river event observed at Dome C with available isotopic vapour observations, and the only one recorded over the 2009–2019 period originating from the Atlantic sector (Petteni et al., 2025). The processes governing the surface vapour δ¹⁸O signal are presented in Fig. 12 for clear-sky conditions (left) and for the AR event (right). Under clear-sky conditions, the diurnal cycle of δ¹⁸O is primarily controlled by the surface sublimation and condensation cycle. During the day, surface snow sublimation enriches surface vapour, which is mixed within the boundary layer by shallow convection. At night, cooling of the air leads to condensation both in the atmosphere and at the surface, resulting in a depletion in the vapour δ¹⁸O. Additionally, the strong atmospheric stratification at night leads to a decoupling between surface and higher-altitude atmospheric layers. At the peak of the AR event, surface sublimation is identified as the primary source of both specific humidity and vapour δ¹⁸O, accounting for approximately 70 % of the total contribution. The second major contribution comes from large-scale moisture advection associated with the AR, accounting for approximately 30 %. These results emphasize that the observed isotopic signal cannot be attributed solely to large-scale atmospheric transport but requires a detailed consideration of local boundary layer processes. Thus, this study deepens our understanding of the local and synoptic processes that shape the vapour isotopic signal in polar regions. We note, however, that in the model there is currently no isotopic fractionation during surface sublimation. While numerous field studies have demonstrated that the sublimation flux is fractionating (Steen-Larsen et al., 2013; Casado et al., 2016; Madsen et al., 2019), the introduction of a fractionation scheme in the model could reduce the contribution of sublimation. Nevertheless, the analysis of the total water budget suggests that surface sublimation makes a substantial contribution to the observed specific humidity anomaly.

To extend our surface analysis, we also performed a water and isotope budget analysis across the entire mixed layer (Fig. S11). Integrating processes over a variable mixed layer height introduces a contribution from the entrainment with upper layers, as explained in Sect. 2.3.2. Applying the same process attribution as for the surface, we find that 49 % of the specific humidity content at the AR peak is attributed to surface sublimation and turbulence, while 38 % is due to advection and 13 % to entrainment (Fig. S12a, b). Since the entrainment term represents a combination of water originating from surface sublimation and advection, we observe that at least half of the mixed layer vapour content comes from surface sublimation, aligning with the surface-level analysis. For δ¹⁸O, the mixed layer analysis indicates that 48 % of the contributions originate from advection, 39 % from surface sublimation, and 13 % from entrainment compared to the mean clear-sky diurnal cycles (−62 ‰) (Fig. S12d). By attributing all entrained water to advection, we obtain a contribution to total vapour of 61 % by advection and 39 % by sublimation, which we consider to be the upper bound for the contribution of advection to the mixed layer isotopic signal.

Finally, we perform a process-based analysis for positive and negative tendencies to estimate not only the sources of vapour, but also the processes driving the variability of specific humidity and δ¹⁸O during the event. To this end, we compared the unperturbed diurnal cycle prior to the AR event (1–17 December), during which we computed the mean cumulative temporal change in specific humidity associated with each process (Eq. 5; Fig. S13a) and in δ¹⁸O associated with each process (Eqs. 6 and 7; Fig. S13c). We then examined the anomaly during the AR event (19–24 December) relative to this reference diurnal-cycle period, for both specific humidity and δ¹⁸O (Fig. S13b, d). For specific humidity, the first 3 d of the event are characterized by a progressive increase in surface sublimation, to which advection contributes from the third day onward, explaining the main humidity peak on 21 December (Fig. S13b). During the night on 21 December, the surface condensation is partly compensated by increased advection. On 22 December, the combination of surface sublimation and advection explains the second peak in specific humidity. Throughout the event, shallow convection and cloud condensation act to reduce the amount of water in the near-surface atmosphere. For δ¹⁸O, we find a slight positive anomaly in surface sublimation during the first 2 d of the event compared with the typical diurnal cycle, followed by a negative surface-flux anomaly from the fourth day onwards. The negative contributions to δ¹⁸O from condensation and shallow convection are weaker than during unperturbed diurnal cycles. In addition, advection decreases the vapour δ¹⁸O in the model during the peak of the event. Because the amplitude of vapour δ¹⁸O during surface sublimation and condensation cycles is overestimated under unperturbed conditions, as shown here and previously discussed in Dutrievoz et al. (2025b) and Ollivier et al. (2025a), the interpretation of isotopic anomalies during the event remains limited, underscoring the need to correct this bias in the LMDZ6iso model (Fig. S13d).

Figure 12 Summary of the processes controlling the modelled isotopic composition of surface vapour at Concordia during classical diurnal cycle (left) and during the AR (right). The processes of surface fluxes and turbulent mixing are shown in orange, cloud condensation and precipitation sublimation is shown in purple, shallow convection is shown in red and advection is in blue.

This study demonstrates that the LMDZ6iso model performs well in representing surface meteorological variables during the AR event. With regard to vapour isotopes, the model accurately captures the maximum vapour δ¹⁸O during the event, however the overestimated amplitude of vapour δ¹⁸O outside the event, discussed above, raises doubts as to whether this agreement reflects the correct physical processes. This overestimated amplitude may result either from excessive daytime enrichment during sublimation and/or from excessive nighttime depletion during condensation. The analysis of humidity and isotopic tendencies using the LMDZ6iso model proves particularly effective in disentangling and quantifying the contributions of processes influencing the simulated isotopic signal. This methodology offers a promising perspective for better understanding the interactions between local mechanisms and large-scale processes within the surface boundary layer, as well as for more effectively interpreting the isotopic signature of atmospheric rivers in water vapour.

To this end, a first step is to improve the representation of isotopic processes in climate models. Recent studies have demonstrated that isotopic fractionation occurs during sublimation (Madsen et al., 2019; Wahl et al., 2021; Hughes et al., 2021), challenging the assumption of non-fractionating fluxes. Given that local sublimation fluxes dominate over large-scale advection in contributing to surface vapour isotopic variations, incorporating isotopic fractionation during sublimation, as in Wahl et al. (2022) and Dietrich et al. (2023), is a promising perspective to improve the simulation of water vapour isotopes over Antarctica, as it is expected to reduce the amplitude of isotopic diurnal cycle variations by decreasing vapour enrichment during sublimation. Consequently, this could lower the isotopic contribution of surface sublimation during the AR. Additionally, improving the representation of snow in the LMDZ6iso model could further advance the accuracy of vapour isotopic composition simulations. Currently, the model averages snowfall from the start of the simulation into a single snow bucket, yielding a long-term mean isotopic value that obscures variability associated with recent events. Consequently, fresh precipitation, which directly influences the isotopic composition of vapour through sublimation, is not adequately represented. To address this limitation, future developments will aim at implementing a more detailed snow bucket scheme capable of capturing the isotopic impact of recent precipitation on surface water vapour.

In parallel, a more detailed investigation of snowfall and surface snow evolution during AR events is also needed. Surface snow is expected to reflect both the isotopically enriched snowfall associated with ARs and wind-blown snow, as wind speed exceeded 10 m s⁻¹ during this specific event. At such intensities, drifting snow can be lifted several hundred metres above the surface (Palm et al., 2011), where it may readily sublimate, potentially changing the isotopic composition of atmospheric water vapour (Wahl et al., 2024). Consequently, the sublimation of airborne snowflakes may exhibit a distinct isotopic signature than that of sublimation flux from the surface. As this process is not currently represented in LMDZ6iso, understanding its contribution to the surface isotopic signal would be a valuable next step.

To assess whether our analysis is representative of the interactions between large-scale and local processes, it would be valuable to analyse additional AR events, both at Concordia and other locations, to determine whether the contributions of local and large-scale processes vary depending on the event and location. In terms of detection algorithms, this trans-Antarctic AR highlights the need for an AR detection framework that is not restricted to 85° S. Future developments in polar AR detection will hopefully address this limitation. In our study, the AR crosses the Antarctic Plateau for 3 d, from its arrival in western Dronning Maud Land on 17 December to its passage over Concordia on 19 December. During this period, the AR moves through an extremely dry atmosphere, where its passage induces a strong temperature anomaly, leading to significant surface sublimation. This process may have progressively modified the isotopic composition of the AR. As a result, the AR could gradually reach an isotopic composition similar to that of surface snow, thereby weakening the contribution of advection to the observed δ¹⁸O anomaly at the surface in Concordia. To test this further, similar studies should be conducted on other AR events with different moisture origins and at stations located at varying distances from the coast. The March 2022 atmospheric river would be an excellent case study, as a +28 ‰ anomaly in water vapour δ¹⁸O was measured at Concordia (Wille et al., 2024b), with the moisture originating from the Indian Ocean. This would help determine whether the isotopic composition of an AR equilibrates with that of surface snow as it travels across the Antarctic Plateau. Moreover, an isotopic analysis of these events using a Lagrangian approach, combined with a process-based decomposition similar to that of Dütsch et al. (2018), would provide a complementary perspective on the processes involved during the transport of atmospheric rivers.

Finally, it would be valuable to gain a more detailed understanding of the isotopic composition of water vapour throughout the vertical structure of the atmosphere. In this study, we have highlighted the strong contribution of surface sublimation to the isotopic anomaly at the surface and within the boundary layer. However, evaluating our model at higher altitudes remains crucial. For example, using observational data, it would be necessary to evaluate daytime shallow convection, as this process mixes the boundary layer and transports moisture from the surface to higher altitudes. Recently, Rozmiarek et al. (2025) used fixed-wing uncrewed aircraft to obtain isotopic observations (δ¹⁸O and δD) up to 1500 m a.g.l. over the northeastern Greenland Ice Sheet during summer 2022. This innovative approach revealed an “inverted-C” vertical structure for δ¹⁸O and δD under cloudy conditions and a “C-shaped” structure for d-excess between the surface and 1500 m altitude. This pattern reflects the dual influence of local-scale processes, such as surface sublimation and large-scale advection. Reproducing this methodology under clear-sky conditions and during AR events in Antarctica and Greenland would provide a valuable extension to our surface and mixed layer analysis. Furthermore, isotopic measurements at higher altitudes would help better constrain the microphysics of mixed-phase clouds, which are frequently observed over the Antarctic Plateau (Ricaud et al., 2020, 2024b).

Appendix A: Computation of the isotopic delta for different process

To compute the isotopic delta as a function of the isotopic delta associated with each process, we use the definition of δ = $\left(R/R_{\mathrm{VSMOW}}-\mathrm{1}\right)\cdot \mathrm{1000}$. In the LMDZ6iso model, water variables are defined by the specific humidity of total water q (in kg kg⁻¹) and the specific humidity of each water isotopes q^iso (in kg kg⁻¹). The specific humidity of total water and of each water isotope can be decomposed as the sum of the contributions from each process:

$$\begin{array}{}\text{(A1)}& {\displaystyle }q=\sum _{\mathrm{process}}q{|}_{\mathrm{process}},\text{(A2)}& {\displaystyle }{q}^{\mathrm{iso}}=\sum _{\mathrm{process}}{q}^{\mathrm{iso}}{|}_{\mathrm{process}},\end{array}$$

where q|_process and q^iso|_process represent the total specific humidity and the isotopic specific humidity respectively, associated with each process (in kg kg⁻¹), computed as the time integral of the isotopic tendencies:

$$\begin{array}{}\text{(A3)}& {\displaystyle }q{|}_{\mathrm{process}}=\underset{t_{\mathrm{init}}}{\overset{t_{\mathrm{final}}}{\int }}{\left.{\displaystyle \frac{\mathrm{d}q}{\mathrm{d}t}}\right|}_{\mathrm{process}}\cdot \mathrm{d}t+q_{\mathrm{init}},\text{(A4)}& {\displaystyle }{q}^{\mathrm{iso}}{|}_{\mathrm{process}}=\underset{t_{\mathrm{init}}}{\overset{t_{\mathrm{final}}}{\int }}{\left.{\displaystyle \frac{\mathrm{d}{q}^{\mathrm{iso}}}{\mathrm{d}t}}\right|}_{\mathrm{process}}\cdot \mathrm{d}t+q_{\mathrm{init}}^{\mathrm{iso}}.\end{array}$$

Using these variables, the isotopic ratio R in mol mol⁻¹ can be computed as follows:

$$\begin{array}{}\text{(A5)}& R={\displaystyle \frac{q^{\mathrm{iso}}/M^{\mathrm{iso}}}{q/M^{\mathrm{16}}}}={\displaystyle \frac{q^{\mathrm{iso}}}{q}}\cdot {\displaystyle \frac{M^{\mathrm{16}}}{M^{\mathrm{iso}}}}=\sum _{\mathrm{process}}{\displaystyle \frac{q^{\mathrm{iso}}{|}_{\mathrm{process}}}{q}}\cdot {\displaystyle \frac{M^{\mathrm{16}}}{M^{\mathrm{iso}}}},\end{array}$$

where M^iso is the molar mass of the water isotope (in kg mol⁻¹), and M¹⁶ is the molar mass of H${}_{\mathrm{2}}^{\mathrm{16}}$O (in kg mol⁻¹). In this equation, we approximated the specific humidity of H${}_{\mathrm{2}}^{\mathrm{16}}$O, q¹⁶, by the specific humidity of total water q.

Using this framework, the δ of water vapour becomes:

$$\begin{array}{}\text{(A6)}& {\displaystyle }\mathit{\delta }=\left({\displaystyle \frac{R^{\mathrm{iso}}}{R_{\mathrm{VSMOW}}}}-\mathrm{1}\right)\cdot \mathrm{1000},\text{(A7)}& {\displaystyle }=\left(\sum _{\mathrm{process}}\left[{\displaystyle \frac{q^{\mathrm{iso}}{|}_{\mathrm{process}}}{q}}\cdot {\displaystyle \frac{\mathrm{1}}{R_{\mathrm{VSMOW}}}}\right]-\mathrm{1}\right)\cdot \mathrm{1000},\text{(A8)}& {\displaystyle }=\left(\sum _{\mathrm{process}}\left[{\displaystyle \frac{q^{\mathrm{iso}}{|}_{\mathrm{process}}}{q}}\cdot {\displaystyle \frac{\mathrm{1}}{R_{\mathrm{VSMOW}}}}-{\displaystyle \frac{q{|}_{\mathrm{process}}}{q}}\right]\right)\cdot \mathrm{1000},\text{(A9)}& {\displaystyle }=\sum _{\mathrm{process}}\left[{\displaystyle \frac{q{|}_{\mathrm{process}}}{q}}\left({\displaystyle \frac{q^{\mathrm{iso}}{|}_{\mathrm{process}}}{q{|}_{\mathrm{process}}}}\cdot {\displaystyle \frac{\mathrm{1}}{R_{\mathrm{VSMOW}}}}-\mathrm{1}\right)\cdot \mathrm{1000}\right],\text{(A10)}& {\displaystyle }=\sum _{\mathrm{process}}\left[{\displaystyle \frac{q{|}_{\mathrm{process}}}{q}}\cdot \mathit{\delta }{|}_{\mathrm{process}}\right].\end{array}$$

From Eqs. (A2) and (A10), we compute the global δ anomaly relative to the mean diurnal cycle value, expressed as a function of the isotopic contributions associated with each process (Eq. 14):

$$\begin{array}{}\text{(A11)}& {\displaystyle }\begin{array}{rl}\mathit{\delta }-{\mathit{\delta }}_{\mathrm{mean}}=& \sum _{\mathrm{process}}\left[{\displaystyle \frac{q{|}_{\mathrm{process}}}{q}}\cdot \mathit{\delta }{|}_{\mathrm{process}}\right]\\ & -\sum _{\mathrm{process}}\left[{\displaystyle \frac{q{|}_{\mathrm{process}}}{q}}\cdot {\mathit{\delta }}_{\mathrm{mean}}\right],\end{array}\text{(A12)}& {\displaystyle }=\sum _{\mathrm{process}}\left[{\displaystyle \frac{q{|}_{\mathrm{process}}}{q}}\cdot \left(\mathit{\delta }{|}_{\mathrm{process}}-{\mathit{\delta }}_{\mathrm{mean}}\right)\right].\end{array}$$