To What Extent Could Water Isotopic Measurements Help Us Understand Model Biases in the Water Cycle over Western Siberia

We evaluate the isotopic composition of water vapor and precipitation simulated by the LMDZ (Laboratoire de Météorologie Dynamique-Zoom) GCM (General Circulation Model) over Siberia using several data sets: TES (Tro-pospheric Emission Spectrometer) and GOSAT (Greenhouse gases Observing SATellite) satellite observations of tropo-spheric water vapor, GNIP (Global Network for Isotopes in Precipitation) and SNIP (Siberian Network for Isotopes in Precipitation) precipitation networks, and daily, in situ measurements of water vapor and precipitation at the Kourovka site in Western Siberia. LMDZ captures the spatial, seasonal and daily variations reasonably well, but it underestimates humidity (q) in summer and overestimates δD in the vapor and precipitation in all seasons. The performance of LMDZ is put in the context of other isotopic models from the SWING2 (Stable Water Intercomparison Group phase 2) models. There is significant spread among models in the simulation of δD, and of the δD-q relationship. This confirms that δD brings additional information compared to q only. We specifically investigate the added value of water isotopic measurements to interpret the warm and dry bias featured by most GCMs over mid and high latitude continents in summer. The analysis of the slopes in δD-q diagrams and of processes controlling δD and q variations suggests that the cause of the dry bias could be either a problem in the large-scale advection transporting too much dry and warm air from the south, or too strong boundary-layer mixing. However , δD-q diagrams using the available data do not tell the full story. Additional measurements would be needed, or a more sophisticated theoretical framework would need to be developed.


Introduction
General circulation models (GCMs) have persistent and systematic biases in the representation of modern climate and its associated water cycle (Meehl et al., 2007).One example of these biases is the systematic warm and dry bias over mid latitude continental regions in summer, which is evident especially in the Great Plains of the United States and in Central and Eastern Europe (Kittel et al., 1997;Cattiaux et al., 2013).This bias may reflect a misrepresentation of the water cycle, land-surface feedbacks (Moberg and Jones, 2004;Bellprat et al., 2013), or cloud feedbacks (Boé, 2013;Cattiaux et al., 2013).This misrepresentation casts some doubts on the climate projections simulated for these regions.In particular, there is evidence for a link between the dry and warm bias in mid-latitude continents in summer for present-day and future projections of evapotranspiration (Boé and Terray, 2008) and Published by Copernicus Publications on behalf of the European Geosciences Union.
temperature extremes (Boberg and Christensen, 2012;Christensen and Boberg, 2012).Credible climate projections thus require a correct simulation of the present-day climate, water cycle, and associated processes.
Several approaches have been proposed to investigate the persistent warm and dry bias over mid-latitude continents in summer.Sensitivity tests on the model physics have shown that biases can be reduced by improving parameterizations or tuning parameters related to the representation of soil processes (Polcher et al., 1996;Cheruy et al., 2013).Experiments performed running GCMs in forecast mode with a realistic initialization of meteorological variables have also suggested that the bias is associated with land-surface interactions (Klein et al., 2006).Here, we investigate to what extent comparing the simulated water stable isotopic composition of water vapor to measurements can help us diagnose the source of model biases over continental regions in summer.
The water molecule has several isotopologues.The most common isotopologue is H 16  2 O (hereafter called H 2 O), but heavier isotopologues are also found: HD 16 O (hereafter called HDO, with D standing for deuterium) and H 18  2 O.The water vapor isotopic composition (e.g., the concentration in HDO and H 18  2 O with respect to H 2 O) is affected during phase changes and conserved during transport.Therefore, the water vapor composition records the various evaporation and condensation processes undergone within an air parcel throughout its history.For example, oceanic evaporation and continental evapotranspiration leave different imprints on the water vapor composition (Salati et al., 1979;Gat and Matsui, 1991;Risi et al., 2013;Simonin et al., 2013).Water vapor and subsequent precipitation is further affected by large-scale circulation (Frankenberg et al., 2009;Galewsky and Hurley, 2010) and convective and cloud processes (Bony et al., 2008;Risi et al., 2008;Yoshimura et al., 2010).In turn, water vapor isotopic measurements can help diagnose biases in the way models simulate the water cycle.For example, isotopic measurements have been used to understand systematic biases in the precipitation frequency in the tropics (Lee et al., 2009) and in subtropical tropospheric humidity (Risi et al., 2012a;Risi et al., 2012b).If biases in evapotranspiration, in largescale moisture transport or in cloud processes contribute to the dry and warm model biases over continental regions in summer, then water isotopic measurements may help diagnose them.
Before we can operationally use water isotopic measurements to diagnose causes of model biases, the first step is to understand what controls the water isotopic composition.In this paper, we focus on this first step, and give insight on how water isotopic measurements could be used for model evaluation.With this aim, we exploit a new set of continuous isotopic measurements in low-level atmospheric water vapor and precipitation at Kourovka (70 km northwest from Ekaterinburg) near the western boundary of Western Siberia.This data set is complemented by satellite measurements (Worden et al., 2007;Frankenberg et al., 2013) and compared with existing isotopic GCMs.In particular, using the LMDZ (Laboratoire de Météorologie Dynamique-Zoom) GCM (Risi et al., 2010c), processes controlling water isotopic composition are investigated and sources of model-data mismatches are discussed.
In Sect.2, data sets and models are described.In Sect.3, the LMDZ GCM is evaluated against satellite and precipitation network data and its performance is put in a broader context through comparison with other GCMs.The ability of LMDZ to represent the water cycle is evaluated in Sect. 4 using continuous measurements of water isotopic composition in Kourovka.In Sect.5, physical processes controlling water vapor and precipitation composition are investigated.We discuss the results and conclude in Sect.6.

Isotopic definitions
Isotope ratios are commonly reported relative to a standard in δ-notation: where the standard used is Vienna Standard Mean Ocean Water (VSMOW), δ represents either δD or δ 18 O expressed in ‰, and R is the ratio of HDO or of H 18 2 O to H 2 O. Deuterium excess is defined with respect to the global meteoric water line (GMWL, δD = 8 • δ 18 O + 10 ‰) (Dansgaard, 1964): d excess = δD − 8 • δ 18 O.Values that fall on the GMWL have a d excess of 10 ‰ by definition.Since equilibrium Rayleigh condensation processes roughly follow the GMWL slope of 8, deviations in d excess can provide information about the environmental conditions during nonequilibrium processes.

LMDZ model and simulations
The dynamical equations of the GCM LMDZ are discretized in a latitude-longitude grid, with a resolution of 2.5 • ×3.75 • and 39 vertical levels.The time step for the resolution of the dynamical equations is 1 min.We use the fourth version of LMDZ (called LMDZ4) which was used for CMIP3 (Coupled Model Intercomparison Project Phase 3, Meehl et al. (2007) feeding the fourth IPCC (International Panel on Climate Change) report.The physical package is described in detail by Hourdin et al. (2006) and called every 30 min.Each grid cell is divided into four sub-surfaces: ocean, land, ice sheet and sea ice.
To ensure a realistic large-scale circulation and daily variability, horizontal winds at each vertical level are nudged by ECMWF reanalyses (Uppala et al., 2005).The bare soil fraction as a function of leaf area index is calculated using d' Orgeval (2006) rather than Ducoudré et al. (1993), which leads to an increase in the bare soil fraction.
yes Table 2. List of the SWING2 models used in this study, their respective expended names and references for isotopic implementation and simulations.Whether the winds are nudged by reanalyses or not ("free") is also indicated.Isotopologues of water (H 16 2 O, H 18 2 O and HDO) are transported and passively mixed by the large-scale advection and various air mass fluxes.Fractionation is considered during condensation into droplets and ice crystals, ocean evaporation and evaporation of rain drops (Bony et al., 2008).The isotope-enabled version of the GCM LMDZ is described in detail by Risi et al. (2010c).

Representation of the land surface
The land surface in LMDZ is represented as a simple bucket model.Land surface evaporation is calculated as a single flux, represents all components of evapotranspiration.No distinction is made between transpiration, bare soil evaporation, or evaporation of water intercepted by the canopy.For water isotopes, we assume that transpiration is the dominant component of evapotranspiration (e.g., Williams et al., 2004;Jasechko et al., 2013).This approximation is especially reasonable in Siberia (Iijima et al., 2014).No fractionation is associated with transpiration (Washburn and Smith, 1934;Barnes and Allison, 1988).Thus we neglect fractionation during evapotranspiration, as in most GCMs (e.g., Hoffmann et al., 1998).
To quantify the impact of neglecting fractionating evaporation, we use a few additional simulations in which LMDZ is coupled with a more sophisticated, state-of-the-art landsurface scheme called ORCHIDEE (Organising Carbon and Hydrology in Dynamic Ecosystems, Ducoudré et al., 1993;Krinner et al., 2005) enabled with water isotopes (Risi, 2009;Risi et al., 2013).We do not use these LMDZ-ORCHIDEE simulations in the core of the paper because these simulations are available for the year 2006 only (after sufficient spin-up in perpetual 2006 conditions).Using LMDZ-ORCHIDEE, we compare simulations in which fractionation during bare soil evaporation is activated or disabled.
The representation of land surface is much simpler than in current coupled models used for CMIP3 (Meehl et al., 2007) or CMIP5 (Taylor et al., 2012).Therefore, some of the conclusions reached in this paper regarding the role of land surface processes might be model-dependent and specific to GCMs with very simple land surface schemes.To check to what extent our results are sensitive to the representation of the land surface, we perform different sensitivity tests on land surface parameters, as listed in Table 1.These sensitivity tests were described in Risi et al. (2013) and were shown to result in large differences in evapotranspiration and in δD variability (Risi et al., 2013).

SWING2 models and simulations
To put LMDZ results into a broader context, we compare with other isotopic GCM simulations archived in the SWING2 (Stable Water INtercomparison Group phase 2) database (http://people.su.se/~cstur/SWING2/).We use nine simulations from seven models (listed in Table 2).Some simulations are nudged by reanalyses and some others are not.More details about these SWING2 simulations can be found in Risi et al. (2012b).We use all available years (between 10 and 30) to calculate climatological averages.The choice of the time period used for the climatological average has a very small impact since all models have at least 20 years available over the past few decades.

Water vapor measurements at Kourovka
The measurements of atmospheric surface water vapor isotopic composition are performed by the Picarro isotopic analyzer L2130-i based on wavelength-scanned cavity ring down spectroscopy (WS-CRDS).This instrument was installed in the Kourovka astronomical observatory (57.037The instrumental uncertainty is conservatively estimated to be 1.4 ‰ for δD and 0.23 ‰ for δ 18 O, respectively (Steen-Larsen et al., 2013).However, at low humidity levels (<1500 ppmv) the instrument reveals a strong increase in uncertainty estimated to be 5.6 ‰ for δD and 0.92 ‰ for δ 18 O at 1000-1500 ppmv and 11.2 ‰ for δD and 1.84 ‰ for δ 18 O at 500-1000 ppmv (Bastrikov et al., 2014).The measurements below 500 ppmv are considered to be beyond the instrument measuring capabilities and are not used in the analysis.
From June 2012 until September 2012, the Picarro calibration system was not working properly due to a leakage.For this reason, the data obtained during this period are associated with a reduced accuracy.We thus the use the δD measurements only, for which this higher uncertainty is acceptable.After the replacement of the faulty elements on 18 September 2012, the instrument has revealed good stability and reliability.A detailed overview of the WS-CRDS measurement system setup, calibration, and maintenance can be found in Bastrikov et al. (2014).In summary, every 6 hours of ambient air measurements are followed by a two-standards calibration lasting 30 min for each reference water standard using Picarro Standards Delivery Module.The liquid standards are vaporized at 140 • C using Picarro Vaporizer Module A0211, then mixed with dried room-air dessicated with Drierite (W.A. Hammond Drierite Company, Ltd., USA) and measured by the analyzer.The water standards are calibrated on the VSMOW-SLAP (Vienna Standard Mean Ocean Water-Standard Light Antarctic Precipitation) scale by accurate laboratory mass spectrometer measurements at LSCE (Laboratoire des Sciences du Climat et de l'Environnement).

Other measurements at Kourovka
In addition to the water vapor measurements, precipitation samples have been collected at Kourovka at the daily time scale since the end of October 2012.Measurements are performed on the Picarro isotopic analyzer L2130-i installed at the Ural Federal University (Ekaterinburg), which is equipped with a liquid water analysis system.Instrumental precision is 1.0 ‰ for δD and 0.1 ‰ for δ 18 O.We use the data from the end of October 2012 up to the end of December 2012.
Basic meteorological measurements are performed on an automatic meteorological station.The model is MetPak-II from Gill Instruments.It is installed 1 m below the Picarro input.We use the temperature measurements to better interpret our q measurements.

Precipitation networks
To put the Kourovka measurements into a broader regional context and to evaluate the capacity of LMDZ to capture the spatial and seasonal patterns of precipitation isotopic composition, we combine two networks: GNIP (Global Network for Isotopes in Precipitation, Rozanski et al., 1993) and SNIP (Siberian Network for Isotopes in Precipitation, Kurita et al., 2004).Both networks sample monthly precipitation.Both δ 18 O and δD are measured in the rain samples, allowing us to calculate d excess.However, in the case of SNIP, the d excess needs to be considered with caution due to possible rain evaporation effects after sampling (Kurita et al., 2004).

Satellite data sets
To evaluate the spatial and seasonal patterns of vapor isotopic composition, we use two satellite data sets measuring tropospheric δD.At present, satellite measurements cannot provide δ 18 O with sufficient precision for d excess calculations to be useful for scientific applications.GOSAT (Greenhouse gases Observing SATellite) measurements enable the retrieval of the total-column water vapor content in both H 2 O and HDO (Frankenberg et al., 2013).Since most of the totalcolumn vapor is in the lower troposphere, column-integrated δD is mainly sensitive to the δD of the boundary layer (BL).
TES (Tropospheric Emission Spectrometer) measurements enable the retrieval of some information on the vertical distribution of q and δD in the troposphere (Worden et al., 2006(Worden et al., , 2007)).Recent processing of the data leads to enhanced sensitivity in northern latitudes and to a higher degree of freedom of the signal (Worden et al., 2012).For clear sky scenes these data can distinguish air-parcels at 900 hPa from those at 400 hPa.A correction is applied on observed δD following the calibration study of Worden et al. (2011).More details on the GOSAT and TES data used and the quality selection criteria can be found in Risi et al. (2013).
In order to rigorously compare LMDZ with GOSAT or TES, we take into account spatiotemporal sampling and we apply averaging kernels to the model outputs to account for the vertical resolution and use of a priori constraints of the satellite retrievals.More details on the model-data comparison methodology are in Risi et al. (2012aRisi et al. ( , 2013)).We focus only on comparing spatial and temporal variations in order to minimize uncertainties related to biases in the satellite data (Worden et al., 2011;Frankenberg et al., 2013).

Theoretical curves in δD-q diagrams
To a first order, in mid and high latitudes, the water vapor isotopic ratio of HDO to H 2 O follows Rayleigh distillation (Dansgaard, 1964): where R v and q are the final isotopic ratio and specific humidity, R v0 and q 0 are the initial isotopic ratio and specific humidity, and α is the fractionation coefficient.
Since R v remains close to unity, δD can be approximated by: Therefore, to first order, δD is tightly related to q.When investigating what we can learn from water isotopic measurements, we need to investigate the added value compared to q, for example by analyzing the relationship of δD as a function of q (Worden et al., 2007;Schneider et al., 2010;Galewsky and Hurley, 2010;Noone et al., 2011;Noone, 2012).For even better consistency with a Rayleigh distillation, here we analyze the relationship of δD as a function of ln(q), in which the Rayleigh distillation is approximately a straight line.Figure 1 shows how different processes may affect the relationship of δD as a function of ln(q).An air mass that cools and loses water vapor by condensation follows the Rayleigh distillation line (red).Mixing between air masses that have undergone different degrees of distillation follows a hyperbolic shape with more enriched values for a given q (green) (Galewsky and Hurley, 2010).Recycling of precipitation, either through local transpiration of soil water or through re-evaporation of falling raindrops follows mixing lines with more depleted values for a given q (blue and magenta).Such plots can be instructive qualitatively.They may, however, be difficult to apply quantitatively to the data due to various factors affecting these curves.For example, the Rayleigh distillation line is sensitive to δD in the initial vapor, which may depend on the existence of convective activity (Jouzel and Koster, 1996;Lawrence et al., 2004).The mixing lines are very sensitive to the q-δD properties of the components Theoretical lines in the δD-q plot.The initial vapor (red circle) is assumed to be formed over ocean, with a surface temperature of 25 • C and surface relative humidity of 70% according to the Merlivat and Jouzel (1979) closure assumption.This vapor follows a Rayleigh distillation following equation 1 (red line).The green lines show mixing lines between the initial vapor and a vapors distilled down to different temperatures: 0 • C (thick line), -5 • C (dashed line) and 5 • C (dotted line).These lines are similar to those in Worden et al. (2007).The blue and magenta lines show mixing lines between vapor distilled down to 10 • C and the re-evaporation of rain.We assume that the rain is in equilibrium with the overlying vapor.The blue line represents the evaporation of rain without any fractionation, as is the case for transpiration of soil water.The magenta lines represents the evaporation of rain with fractionation, following ?and Bony et al. (2008): re-evaporation of rain in an equal amount of water vapor with air humidity of 70% (solid), in an equal amount of water vapor with 90% humidity (dotted) or in a twice larger amount of water vapor with 70% humidity (dashed).Note that these re-evaporation lines are different from those in , because they neglected the evolution of the raindrop composition during re-evaporation, whereas we do a more precise calculation, and because they plotted a combination of distillation and re-evaporation, whereas we plotted the effect of re-evaporation only.
figure Figure 1.Theoretical lines in the δD-q plot.The initial vapor (red circle) is assumed to be formed over ocean, with a surface temperature of 25 • C and surface relative humidity of 70 % according to the Merlivat and Jouzel (1979) closure assumption.This vapor follows a Rayleigh distillation following equation 1 (red line).The green lines show mixing lines between the initial vapor and a vapors distilled down to different temperatures: 0 • C (thick line), −5 • C (dashed line) and 5 • C (dotted line).These lines are similar to those in Worden et al. (2007).The blue and magenta lines show mixing lines between vapor distilled down to 10 • C and the re-evaporation of rain.We assume that the rain is in equilibrium with the overlying vapor.The blue line represents the evaporation of rain without any fractionation, as is the case for transpiration of soil water.The magenta lines represents the evaporation of rain with fractionation, following Stewart (1975) and Bony et al. (2008): re-evaporation of rain in an equal amount of water vapor with air humidity of 70 % (solid), in an equal amount of water vapor with 90 % humidity (dotted) or in a twice larger amount of water vapor with 70 % humidity (dashed).Note that these re-evaporation lines are different from those in Worden et al. (2007), because they neglected the evolution of the raindrop composition during re-evaporation, whereas we do a more precise calculation, and because they plotted a combination of distillation and re-evaporation, whereas we plotted the effect of re-evaporation only.
that are being mixed (compare the different green lines in Fig. 1).Re-evaporation lines are sensitive to evaporation conditions (compare the different magenta lines in Fig. 1).Note that in the following, we discuss the relationship between q and water vapor δD in the near-surface air, or in the total column, i.e., between the precipitable water (Q) and the column-integrated water vapor δD.

Consequences for the interpretation of model-data differences in δD
Model-data differences in δD ( δD) can be interpreted in the light of δD-q diagrams: if δD is misrepresented, it could be either because q is misrepresented, or because the δD-q relationship is misrepresented.For example, if δD is controlled by Rayleigh distillation and LMDZ misrepresents the intensity of this distillation, then LMDZ will overestimate both q and δD.On days when LMDZ overestimates q the most, LMDZ will overestimate δD the most.In addition, the slope of δD vs. ln(q) would be the same as that of δD vs. ln(q): indeed, based on Eq. ( 1) and applying the   approximation of Eq. ( 2), δD relates to q following: Similarly, if δD is controlled by mixing and if LMDZ misrepresents the proportion of this mixing, it can be shown that LMDZ will overestimate δD the most on days when it overestimates q the most and that the δD-ln(q) slope would also be similar to the δD-ln(q) slope.
In contrast, if LMDZ misrepresents the initial δD of water vapor, there will be a systematic offset between the data and the model δD, and no particular relationship is expected between δD and ln(q).

Model evaluation of spatial and seasonal variations
To evaluate the spatial and seasonal variations over Siberia simulated by LMDZ, we use satellite and precipitation network data.

Spatial variations
Figure 2 shows the maps of Q, and of annual-mean vertically integrated water vapor δD.Both GOSAT and TES data exhibit the temperature effect, with δD decreasing with latitude, and the continental effect, with δD decreasing along trajectories from west to east (Fig. 2a, c).To first order, LMDZ captures these patterns qualitatively well (Fig. 2b, d, f).´ µ δD  The maps for precipitation δD are consistent with those for vapor δD (Fig. 3a-b).The scarcity of precipitation d excess data makes it difficult to extract any clear spatial signal, although some poleward decreasing d excess trend can be noticed (Fig. 3c).LMDZ captures this pattern qualitatively well (Fig. 3d).For a more quantitative evaluation, Figs. 4 and 5 show north-south and west-east transects around Kourovka and are described below.

Latitudinal gradients in Q and in δD
Observed Q decreases with latitude, as expected from the decrease in temperature.This is very well simulated by LMDZ (Fig. 4a-b).
For satellite δD data, we subtract domain average values to focus on spatiotemporal variations independently of possible systematic biases in the data.The poleward depletion associated with the temperature effect can clearly be seen in    GOSAT, TES, and precipitation network data sets (Fig. 4ce).Compared to both satellite data sets, LMDZ underestimates the latitudinal gradient of δD by about 30 % (Fig. 4cd).This bias was already noticed by Risi et al. (2012a).We investigate the latitudinal gradients at different altitudes using the TES data, which has some vertical information (not shown).LMDZ underestimates the latitudinal gradient from the surface to 800 hPa, and slightly overestimates it from 800 to 550 hPa.800 hPa corresponds approximately to the top of the BL.Above 550 hPa, the TES sensitivity to δD becomes too low to draw conclusions.Therefore, the underestimation of the latitudinal gradient originates from the BL.For example, it could be due to a wrong representation of the latitudinal variations in evapotranspiration, with too much evapotranspiration in the northern part of Siberia, or too little evapotranspiration in the southern part.A seasonal analysis of the latitudinal gradients shows that the latitudinal gradient is underestimated only in summer.This could be consistent with the hypothesis that evapotranspiration is misrepresented, since evapotranspiration occurs mainly in summer.However, sensitivity tests using LMDZ-ORCHIDEE suggest that purely atmospheric processes are responsible for this bias (Sect. 3.3.4).
For the precipitation, the data is noisier due to the scarcity of observations.We cannot detect any obvious bias in the simulation of the latitudinal gradient of δD compared to the data (Fig. 4e).

East-west gradients in Q and δD
Q decreases westward over Siberia, reflecting the progressive rainout along air masses trajectories towards the east.LMDZ overestimates slightly this decrease compared to both GOSAT and TES (Fig. 5a-b).
GOSAT, TES and the precipitation networks all feature an eastward decrease of δD (Fig. 5c-e).This decrease has been attributed to the so-called "continental effect" (Rozanski et al., 1993).This effect is due to the fact that as air masses move eastward, they lose heavy isotopes through precipitation.A fraction of these heavy isotopes are returned to the atmosphere through evapotranspiration (Gat, 2000), but the remaining is lost through river runoff.This is why the amount of continental recycling (i.e., the fraction of the precipitation water which is returned to the atmosphere through evapotranspiration) is known to modulate the continental effect (i.e., the inland depletion of water vapor and precipitation) (Salati et al., 1979;Kurita et al., 2004).
The eastward depletion associated with the continental effect is very well captured by LMDZ compared to GOSAT, TES, and precipitation network data sets (Fig. 5a-c).If the continental effect is modulated by continental recycling, then this suggest that LMDZ represents the east-west gradient in evapotranspiration satisfactorily.

Spatial variations in d excess
For d excess, the data looks much noisier than LMDZ.This could be due to the large uncertainty in the d excess

δD
´ µ ØÓØ Ð ÓÐÙÑÒ Fig. 6.Seasonal cycles of observed (black) and simulated (green) precipitable water observed by GOSAT (a) and by TES (b), δD in the total column water vapor as observed by GOSAT (c) and by TES (d), δD (e) and d-excess (f) in the surface precipitation.We average over the region of Kourovka (50 For satellite datasets, the annual-mean δD is subtracted to focus on the seasonal variations.
For GNIP/SNIP datasets, the spatial standard deviation is also plotted as an envelope.LMDZ outputs are collocated with each dataset and averaging kernels are applied to compare with satellite datasets.For satellite data sets, the annual-mean δD is subtracted to focus on the seasonal variations.
For GNIP/SNIP data sets, the spatial standard deviation is also plotted as an envelope.LMDZ outputs are collocated with each data set and averaging kernels are applied to compare with satellite data sets.
measurement.The extent of post-sampling evaporation effects are difficult to quantify, but they could reach several ‰ (Kurita et al., 2004), which is of the same order of magnitude as the north-south d-excess gradient simulated by LMDZ.
The apparent data noise could also be due to the potentially large spatial heterogeneity of d excess at the scale of a few kilometers: for example, the local surface type could affect d excess (Welp et al., 2012).LMDZ cannot capture this heterogeneity.This could also explain why LMDZ looks smoother than the data.
In spite of the noisiness in the data, a decreasing trend with latitude can be observed.This could be associated with the Rayleigh distillation, which first decreases d excess until about −20 • C and increases it below (Masson-Delmotte et al., 2008).In Siberia, only the decreasing trend can be seen because temperature are infrequently below −20 • C. LMDZ captures the d excess decrease with latitude, with a decrease from 14 ‰ at 35 • N to 5 ‰ at 70 • N (Fig. 4d).
No clear eastward trend can be noticed in observed precipitation d excess.The large noise and lack of clear signal in the data makes it difficult to assess whether the model is consistent or not with the data.

Seasonal variations
The observed seasonal cycle of Q follows that of temperature, with moister air in summer (Fig. 6a-b).LMDZ captures this seasonal cycle qualitatively, but is too dry in summer (consistent with the common dry bias previously mentioned, Kittel et al., 1997;Cattiaux et al., 2013) and is too moist in autumn and winter.The LMDZ seasonal cycle is delayed by a few weeks.The correct simulation of Q noticed in the previous section in annual mean hides seasonal discrepancies.
As is common in most of the Northern Hemisphere, observed δD is more enriched in summer and more depleted in winter (Fig. 6c-e), consistent with the temperature effect.LMDZ captures this seasonality qualitatively well, though δD variations are underestimated compared to both satellite data sets (Fig. 6c-d).The LMDZ seasonal cycle is delayed by a few weeks, consistent with the delay noticed for Q.These problems do not appear in the precipitation networks (Fig. 6e).Both TES and GOSAT have fewer usable data points in the winter.It is possible that there are too few samples to make a robust comparison.
The precipitation networks show large spatial variability for d excess (Fig. 6f).However the comparison suggests that LMDZ misrepresents the observed d excess minimum in summer.

Possible causes of model-data differences
In the previous section, we identified three main modeldata differences: (1) LMDZ underestimates Q in summer, (2) LMDZ underestimates the seasonal variation in δD and (3) LMDZ underestimates the latitudinal gradient in δD.For δD, a common bias of GOSAT and TES towards overestimating δD variations is unlikely given that they used www.atmos-chem-phys.net/14/9807/2014/Atmos.Chem.Phys., 14, 9807-9830, 2014 Table 3. Relationships between precipitable water (Q) and total column water vapor δD as observed by GOSAT or TES and as simulated by LMDZ.We show the relationships for the latitudinal transects as plotted in Fig. 4 and for seasonal cycles as plotted in Fig. 6.For each relationship, we give successively the Q range (minimum and maximum values, in kg m −2 ), the δD range (maximum minus minimum values, in ‰, we do not give absolute values since they might be subject to calibration issues) and the slope of the δD-ln(Q) relationship (in ‰).The correlation coefficients for the δD-ln(Q) relationship are all between 0.82 and 0.84.When we compare LMDZ to GOSAT or TES, we collocate model outputs with each data set and apply the adequate averaging kernels.independent wavelengths and retrieval algorithms.This section explores possible causes for these model biases.

Spatial pattern of the summer bias
The latitudinal gradient in δD is underestimated the most in summer.Therefore, the underestimation of both the δD seasonal variations and the annual-mean δD latitudinal gradient could be due to the underestimation of the δD latitudinal gradient in summer.Figure 7 shows the model-data difference for Q and δD compared to GOSAT and TES.The underestimation of the δD latitudinal gradient is apparent (Fig. 7c-d).This is associated with an underestimation of the Q latitudinal gradient (Fig. 7a-b).The lack of calibration for δD data prevents us to decide whether it is the northern part that is too enriched or it is the southern part that is too depleted, but both GOSAT and TES data for Q shows that the northern part is too moist and the southern part is too dry in LMDZ.As air masses move towards the southeast, their humidity decrease too much in LMDZ compared to the data along air mass trajectories.This suggests that in LMDZ, there is an excess of dehydrating processes, or a deficit of moistening processes, in Western Siberia.The location of Kourovka is ideal to in-vestigate these processes since it is where the air masses start to become too dry along their trajectories.

δD-Q diagrams
As explained in Sect.2.4, the underestimation of δD variations may be associated either with an underestimation of q variations or with an underestimation of the slope of δD-Q slope.Table 3 attempts to separate these two effects.
For the latitudinal gradient, the simulated Q range is comparable with observations, but the δD-Q slope is weaker in LMDZ (Fig. 4a-b, Table 3).Therefore, the δD discrepancies are not due to temperature or Q variations, but rather to the type of hydrological processes (condensation, mixing or re-evaporation).For example, there could be too much diffusion in the advection scheme (Risi et al., 2012a;Risi et al., 2012b) or an over-estimated continental recycling (Risi et al., 2013).Sensitivity tests with ORCHIDEE will however suggest purely atmospheric processes are responsible for the latitudinal gradient mismatch (Sect. 3.3.4).
For the seasonal cycle, LMDZ underestimates Q variations, consistent with the dry bias in summer (Figs.6a-d, 3).In addition, the δD-Q slope is weaker in LMDZ.Therefore,   the underestimation of the δD seasonal cycle is due to both an underestimation of the Q seasonal cycle and of the δD-Q slope.

Impact of the representation of fractionating evapotranspiration
The fact that we neglect the fractionation during the bare soil evaporation component of evapotranspiration may contribute to underestimate the latitudinal gradient of δD near the surface.To quantify this, we compare two LMDZ simulations coupled with ORCHIDEE in which the fractionation during bare soil evaporation is activated (OR) or disabled (ORnofrac).When there is isotopic fractionation during bare soil evaporation, less heavy isotopes are recycled back into the atmosphere.Therefore, we expect a stronger continental effect (Salati et al., 1979), i.e., a stronger east-west gradient in δD.The OR simulation indeed shows a slightly larger east-west gradient in δD than the ORnofrac simulation: −5.4‰ / 10 • compared −5.1 ‰ / 10 • (Table 4).This difference is too small to consider fractionating evaporation as a key improvement.The north-south gradient in δD is not affected (Table 4).This suggests that the north-south gradient is con-trolled by purely atmospheric processes, such as progressive condensation along air mass trajectories.
Regarding d excess, when there is fractionation during bare soil evaporation, the evaporation of HDO is favored relatively to the evaporation of H 18 2 O.This increases d excess in the recycled water vapor (Gat and Matsui, 1991;Risi et al., 2013).The OR simulation indeed shows a significantly larger east-west gradient in δD than the ORnofrac simulation: 0.7‰ / 10 • compared 0.4‰ / 10 • (Table 4).However, the noisiness in the GNIP/SNIP data does not allow one to judge whether considering fractionating evaporation is an improvement.
To summarize, neglecting isotopic fractionation during bare soil evaporation does not appear to be a major caveat of our study.

Impact of the representation of the land surface
To check to what extent our results are sensitive to the representation of the land surface, we compare different sensitivity tests using LMDZ-ORCHIDEE sensitivity tests.We find that the longitudinal q gradient is relatively sensitive to the representation of the land surface (Table 4).In contrast, the latitudinal q gradient is much more robust among tests.This suggests that land surface processes are important for   the Siberian climate, although the latitudinal gradient is controlled mainly by purely atmospheric processes.The latitudinal gradient in δD is a remarkably robust feature of LMDZ and LMDZ-ORCHIDEE simulations.This supports once again that the latitudinal gradient is controlled mainly by purely atmospheric processes, and cannot be used to evaluate the simulation of continental recycling.In contrast, the east-west gradient of water vapor δD is significantly sensitive to the representation of the land surface, with values varying from −5.4 ‰ / 10 • to −7.4 ‰ / 10 • .
The east-west gradient in d excess are very sensitive to the representation of the land surface, with values varying from 0.4 to 1.4 ‰ / 10 • depending on the LMDZ-ORCHIDEE simulations.This suggests that d excess measurements along an east-west transect could be very useful to evaluate the representation of land surface processes.However, the noise in the GNIP/SNIP d excess data is too large compared to the east-west variations.
We note that LMDZ and LMDZ-ORCHIDEE simulations provide significantly different results.We thus reiterate our warning that some of the conclusions reached in this paper might be specific to GCMs with very simple land surface schemes.However, the good match between our LMDZ simulation and TES and GOSAT data suggests that the continental recycling in this simulation is reasonably well represented.

Comparison with other models
As for LMDZ, the European Centre Hamburg Model (ECHAM) model was already compared to the same satellite data sets and showed a good capacity for simulating the observed spatial variations (Butzin et al., 2014).To put the LMDZ results into a broader context, we compare them with those of other model simulations from the SWING2 archive (Fig. 8).Since daily outputs are not available, we cannot collocate the model outputs with the satellite data sets and we cannot convolve them with the appropriate averaging kernels.Therefore, in this section, we just compare the climatological averages of the different models with each other, without any comparison with TES or with GOSAT.SWING2 simulations exhibit a large spread in their representation of the latitudinal gradient in the mid troposphere (Fig. 8a), with gradients ranging from 11 ‰ / 10 • to 25 ‰ / 10 • .Compared to the other models, LMDZ does not simulate the weakest gradient, suggesting that other models might share the same tendency to underestimate this gradient compared to the data.
What might cause this spread?First, some of this spread is related to the simulated q, which reflects the basic climatology.For example, the steepest latitudinal gradient in δD is featured by the free (i.e., not nudged) simulation of LMDZ, which is associated with the largest range of q along the latitudinal gradient (Fig. 8b).The weakest latitudinal gradient in δD is featured by the CAM2 model, which is associated with the smallest range of q along the latitudinal gradient, while both models have a similar δD-q slope (Fig. 8b).Second, Table 5. Daily correlation coefficient (r) and its significance (p value) between observations at Kourovka and LMDZ simulation results for humidity (q), water vapor δD, δ 18 O, d excess and differences of these variables between precipitation and water vapor.During the December-January-February (DJF) season, data is available only for December.If the p value is lower than 5 %, we assume that the correlation coefficient is statistically significant.N/S means "not significant".

Phase period Vapor Vapor
Vapor some of the spread is related to the δD-q slope, which reflects more subtle physical processes.For example, one of the steepest latitudinal gradients in δD is featured by the Hadley Centre Atmospheric Model (HadAM), which is associated with the steepest δD-q slope (Fig. 8b).The reverse is true for isoGSM (Isotopes-incorporated Global Spectral Model).SWING2 simulations also exhibit a large spread in their representation of the seasonal cycle of δD at Kourovka, with seasonal magnitudes ranging from 40 to 160 ‰ (Fig. 8c).Again, LMDZ does not simulate the weakest seasonality among the models, suggesting that other models might share the same seasonal bias.There is no relationship between the amplitude in the seasonal cycle in δD and that in q among the different SWING2 models (Fig. 8c).This shows that the seasonality in δD reflects some physical processes that are not immediately detectable when looking at conventional variables, thus supporting the added value of water isotopic measurements to better evaluate these processes.

Specific humidity and isotopic composition of the water vapor
Continuous measurements of water vapor isotopic composition from Kourovka with the Picarro instrument allow us to analyze daily time series from April to December 2012.A comparison of observations with the results of LMDZ simulation of water vapor for 2012 at Kourovka is shown in Fig. 9.
For LMDZ, we use outputs from the closest grid point from Kourovka.Since LMDZ is nudged by reanalyses, it captures the daily variations in circulation, so that it is possible to make a day-to-day comparison.The LMDZ results correlate very well with the observations of both q and δD for the whole period (Table 5).LMDZ underestimates q from June to November and the difference between observations and model reaches a maximum (8 g kg −1 ) in August.This is consistent with the summer dry bias of LMDZ discussed earlier.A dry bias was also noticed in the free troposphere compared to the IASI (In-frared Atmospheric Sounding Interferometer) data (Pommier et al., 2014).The daily model-data correlation for q is better in winter and rises to 0.94.
LMDZ almost overestimates observed δD and δ 18 O, with the exception of December, where there is very low q.The fact that the LMDZ results are too enriched compared to observations is consistent with the fact that the LMDZ results are too enriched compared to SNIP observations (Figs. 4 and  5).Compared to both data sets, the LMDZ results are on average about 20 ‰ too enriched.
As previously mentioned in Sect.2.3.1, we decide not to use the data for d excess from April to the end of September due to a leakage in the Picarro calibration system.Confidence in d excess for observations is in doubt for this period, and we limit the analysis of these data to the last 3 months, October to December.For this period modeled d excess show fewer variations than observed values and remain in the range of 9 to 15 ‰.An underestimation of daily d excess variability by GCMs had already been noticed for other mid-latitude sites (Risi et al., 2012a).Despite the large differences in values, there is a good qualitative agreement of the results in November.For example, LMDZ captures the d excess increase in the end of October, the decrease in early November and the increase in late November.
The vapor is collected at 8 m above ground level.In the model, the first layer is 70 m thick and the model outputs represent vapor integrated over this depth.The difference in vertical footprint could explain some model-data discrepancy, though we expect this effect to be small.For an average δD gradient of −15 ‰ km −1 , the expected different between 8 m and 33 m (mean altitude of the LMDZ first level) is only 0.3 ‰.
To summarize, LMDZ captures well the temporal variations observed in ground-based measurements, especially in the winter.During the summer, the model underestimates q.All year long, δD is systematically too enriched.

Difference of isotopic composition between precipitation and water vapor
Fig. 10 shows the daily precipitation amount and difference in δD and d excess between precipitation and water vapor: The theoretical estimation of this dif-ference, assuming isotopic equilibrium between precipitation and vapor, is shown in green.Equilibrium fractionation coefficients between vapor and liquid water or ice are calculated following Majoube (1971a, b); Merlivat and Nief (1967).We take into account kinetic effects during snow formation, following Jouzel and Merlivat (1984).
Precipitation is snow, except during the first 3 days.Its amount is well captured by the model (Fig. 10a).
We find significant correlation between the model and observations for δD p − δD v (Table 5).This good agreement is probably due to the fact that most of the precipitation is snow.The isotopic composition of snow is easier to simulate than that of rain because it is less affected by post-condensational processes.Observed precipitation is on average about 80 ‰ more enriched in δD than the water vapor at the surface (Fig. 9a, red dots).The LMDZ results are consistent with these observations, though there are some model-data differences for individual events (Fig. 10a, blue dots).If the snow was in equilibrium with the surface water vapor, it would be about 115 ‰ more enriched (Fig. 10a, green curve).The δD p − δD v of 80 ‰ can be explained by the fact that snow forms from vapor at a higher altitude.To check this quantitatively, we calculate the daily condensation altitude as the average altitude weighted by the condensation rate.In LMDZ, the condensation occurs at an altitude of 2 km on average.δD decreases with altitude and the vertical gradient between the surface and 2 km is about −15 ‰ km −1 .Therefore, most of the snow forms on average from water vapor which is about 30 ‰ more depleted than at the surface.The fact that the observed snow is about 35 ‰ less enriched than predicted based on equilibrium with the surface water vapor is thus consistent with the fact that the snow is formed from water vapor at 2 km.In addition, daily variations in δD p − δD v reflect variations in condensation altitude: in LMDZ, the correlation between δD p − δD v and the condensation altitude for the entire period is −0.51.
In Sect.3.1, we showed that LMDZ simulates a precipitation that is systematically too enriched.The fact that LMDZ correctly simulates δD p − δD v confirms that the cause of the overestimation of the precipitation δD (both at Kourovka and on all GNIP and SNIP sites) is the overestimation of the δD of the vapor from which the precipitation is formed.
Observed d p − d v is 3 ‰ on average (Fig. 10b).LMDZ simulates d p − d v values of −2 ‰ on average.If the snow was in equilibrium with the surface water vapor, d p − d v would be about −5 ‰.The d excess increases with altitude and the vertical gradient between the surface and 2 km is 1 ‰ km −1 .Taking this effect into account, d p − d v should be about -3 ‰.This theoretical estimate is very consistent with what is simulated by LMDZ.Why observed d p − d v is 3 ‰ rather than −3 ‰ could be due to microphysical processes or post-condensational processes.However, the large spread of d p −d v values prevent us from concluding for sure that observations are inconsistent with LMDZ and with the theoretical estimate.5 Processes controlling water vapor and δD

Understanding the simulated evolution in humidity and δD
We showed that LMDZ reproduces well, at least qualitatively, the seasonal (Sect.3.2, Fig. 6) and daily (Sect.4.1, Fig. 9) variations in q and water vapor δD in the lower troposphere.Therefore, we can make use of LMDZ simulations to investigate in more detail the physical processes controlling these variations.

Method based on the tendency analysis
The different processes affecting q and δD at the surface in the model in general are schematically illustrated in Fig. 11.Below we focus in more detail on each of them.
1. Vertical and horizontal advection by the large-scale winds (red).Since we focus on surface q and δD, only the horizontal component of the advection is nonzero.Horizontal advection may moisten or dehydrate the air depending on the horizontal q gradients and on the direction of the wind.Similarly, it may enrich or deplete the water vapor.

Deep convection (green)
. This represents the effect of vertical motions and phase changes in convective systems.At the surface, convection dehydrates the air through subsident motions such as unsaturated downdrafts (Zipser, 1977).This subsidence also has a depleting effect on water vapor (Risi et al., 2008(Risi et al., , 2010b)).Convection may also moisten the lower levels through rain re-evaporation (Folkins and Martin, 2005).In this case, its effect on δD can be either enriching or depleting (Worden et al., 2007;Risi et al., 2010a;Field et al., 2010).
www.atmos-chem-phys.net/14/9807/2014/Atmos.Chem.Phys., 14, 9807-9830, 2014 3. Large-scale condensation and re-evaporation (blue).This represents the effect of phase changes occurring outside of the convective systems, such as in frontal systems or stratiform clouds.Near the surface, it may moisten the BL through the evaporation of the precipitation formed by such systems or clouds.As for the evaporation of convective rain, this can either enrich or deplete the water vapor.Occasionally, it may dehydrate and deplete the vapor through in situ condensation in fogs.
4. Surface evaporation and BL processes (cyan).Our model diagnostics do not allow for the direct separation of these two effects despite their different effects on q and on δD of near-surface vapor.Surface evaporation moistens the BL.Since we assume no fractionation during land surface evaporation, this always has an enriching effect (Risi et al., 2013).In contrast, BL mixing is associated with vertical redistribution of moisture, which has a dehydrating and depleting effect on the lower levels.To qualitatively separate the effect of surface evaporation and BL processes, we use the fact that surface evaporation (which is a model output) and BL mixing (which is expected to be more active in summer, especially during the warmest days) have opposite effects.For example, if the "surface evaporation and BL processes" are more moistening when both surface evaporation and BL mixing are stronger, then we suggest that surface evaporation drives the moistening effect.In contrast, if the "surface evaporation and BL processes" are more moistening when both surface evaporation and BL mixing are weaker, then we suggest that the weaker BL mixing drives the moistening effect.
In LMDZ, the temporal derivatives of q and δD are computed as the sum of different processes.Their contribution is shown in Fig. 12 and in Table 6.

Contribution to humidity variations
At the seasonal scale, the increase of 13 g kg −1 in q from spring to summer (Fig. 8a) is associated with a positive and increasing contribution in the "surface evaporation and BL processes" component, with a maximum of 6 (g kg −1 ) day −1 in July (Fig. 12a).Since we expect BL mixing to be stronger in summer, surface evaporation likely drives the moistening in spring.This increased moistening is partly compensated by an increased dehydration caused by summertime deep convection and large-scale advection (about 4 (g kg −1 ) day −1 for both).In autumn, the decrease in q is associated with a negative contribution of the "surface evaporation and BL processes" component.Since the BL is expected to be most active in summer, the decrease in surface evaporation likely drives the dehydration in autumn.Advection processes also take part in this decrease with 1 (g kg −1 ) day −1 .At the daily scale, in summer, q variations correlate the best with the "horizontal advection" component (r = 0.40, Table 6).In addition, increases in q are associated with large values in the "surface evaporation and BL processes" component.Therefore, surface evaporation also drives the daily variability of q in the summer.This is consistent with the important role of continental recycling on q variations diagnosed by Risi et al. (2013).This consistency could be explained by the fact that the model is the same.
In winter, variations in q are also mainly associated with variations in the large-scale advection component and "surface evaporation and BL processes".This is confirmed by the good correlation between q and "advection" and "surface evaporation and BL" tendencies and (r = 0.65 and 0.67 respectively, Table 6).The importance of advection is consistent with the predominance of synoptic disturbances which modify the direction of the winds.Since evaporation is low during winter, "surface evaporation and BL processes" are likely associated mostly with BL processes.This temporal derivative can be decomposed into the tendencies from different processes (illustrated in Fig. 11): horizontal advection (red); deep convection (green); large-scale condensation and reevaporation (blue); surface evaporation and boundary-layer dynamics (cyan).We used a 5-day filter to focus on seasonal and synoptic variations.Table 7. Model-data difference in near-surface air specific humidity in g kg −1 and its different contributions associated with temperature (T ) or with relative humidity (RH).The sign refers to the modeldata difference.

Season
DJF JJA Mean q 0.52 −2.0 Mean q attributed to T 0.33 2.3 Mean q attributed to RH 0.13 −3.3 Correlation between q and T 0.91 −0.52 Correlation between q and RH 0.68 0.80 Correlation between q and q −0.71 −0.77

Contribution to δD variations
For δD, the contributions of the different processes to the time derivative mirror those for q (Fig. 11b): most moisten-ing processes act to enrich the water vapor, and vice versa.However, the relative amplitudes of the different processes are different.This supports the idea that δD provides some independent information compared to q.At the seasonal scale, surface evaporation contributes to the increase in δD in spring and to its decrease in summer, consistent with what we explained for q.
In the winter, "surface evaporation and BL processes" make the biggest contribution to δD variations: they are well correlated, with r = 0.59 (Table 6).However, in early spring, large-scale advection also plays a significant role in the enrichment of the water vapor (20-30 ‰ day −1 ).At the daily level, surface evaporation drives the spikes of δD in summer with maximum in July with 30 ‰ day −1 , and large-scale advection is mainly responsible for the δD fluctuations in winter (r = 0.44), consistent with our explanation for q.In summer, the large-scale condensation also plays a slightly enriching role, due to the re-evaporation of rain in the BL.In spring and autumn, large-scale condensation strongly depletes the water vapor on some days.This is consistent with the formation of fog (Noone et al., 2013).These fog events occur on days with strong surface evaporation but relatively stable and shallow BL.On these days, the depletion by largescale condensation compensates for the enrichment by surface evaporation.Thus we suggest that the BL dynamics and the large-scale advection control the daily variability of q and isotopic composition especially in winter.In summer, surface evaporation also plays an important role.

Interpreting model-data differences
The goal of this section is to interpret model-data differences in water isotopic composition, and to investigate what we can learn from water isotopic measurements about the representation of processes in models.We investigate first the modeldata differences in q, and then the model-data differences in the δD-ln(q) relationship and how they contribute to modeldata differences in δD.

Model-data differences in specific humidity
LMDZ tends to estimate q correctly during winter, and underestimate it during summer (Fig. 9a).The specific humidity (q) can be expressed as a function of temperature (T ) and relative humidity (RH): q = RH • q s (T ), with q s being the specific humidity at saturation.At the large scale, T reflects mainly dynamical and radiative processes, whereas RH reflects mainly dynamical, surface and cloud processes.T and RH may also reflect small scale processes that are not represented on the coarse grid of LMDZ.
Table 7 shows that on average in winter, LMDZ overestimates q because it overestimates both the temperature (accounting for 72 % of the overestimation of q) and the relative humidity, compared to meteorological data.LMDZ overestimates q the most on days when it overestimates T the most.This occurs mainly during the driest days in the observations.One hypothesis is that LMDZ overestimates the advection of warm and moist air from the southwest during these days.
In summer, Table 7 shows that on average, LMDZ underestimates q because it strongly underestimates the relative humidity, although it overestimates the temperature.LMDZ underestimates q the most on days when it is too warm.Also, LMDZ underestimates q the most during the moistest days in observations.One hypothesis is that LMDZ overestimates the advection of warm and dry air from the south.Another hypothesis is that LMDZ underestimates the surface evapotranspiration.More evapotranspiration would both moisten the atmosphere and cool the surface.
In the next subsection, we investigate whether water isotopic measurements can help us test these hypotheses.

Model-data differences in the δD-ln(q) relationship
Figure 13a shows daily vapor δD as a function of ln(q).The surface observations cluster along a line, consistent with δD varying along a Rayleigh distillation line (Sect.2.4).The correlations are high (Table 8).We compare the observed δD-ln(q) slope with that predicted by Rayleigh distillation (Eq.1).We assume that the initial temperature is 25 • C, the initial RH is 70 % and that the initial δD is predicted by the Merlivat and Jouzel (1979) closure using these temperature and RH conditions.The following results are not strongly sensitive to these assumptions.We calculate the fractionation coefficient using the same temperature as the one used to calculate the saturation specific humidity.We calculate that the theoretical δD-ln(q) slope should be 80 to 120 ‰.The observed slopes are twice lower than expected from the theory: 40 to 60 ‰.This suggests that δD does not follow a pure distillation line.It might also be affected by mixing processes (Hurley et al., 2012).If δD is affected by mixing processes, then δD is expected to vary linearly with 1 / q.To test this hypothesis, we calculate the correlations of δD as a function of 1 / q (Table 8).Correlations are slightly lower but still very high.Therefore, the influence of mixing processes on δD cannot be completely excluded.
Figure 13 and Table 8 show that LMDZ captures the δDln(q) slopes reasonably well in summer and in winter.This suggests that the Rayleigh distillation processes, or mixing processes (if they exist), are well represented.However, δD values are systematically offset: LMDZ is systematically more enriched by about 15 ‰ than observations, for a given q.This suggests that the origin of water vapor that undergoes the distillation is not properly represented.
We analyze the relationship between model-data differences in δD and in ln(q) (Fig. 13b).In winter, there is no systematic relationship between model-data differences in δD and those in ln(q).This suggests that the too-large LMDZ enrichment is not related to a Rayleigh distillation that is too weak.Rather, some processes in LMDZ, independent of the  .refers to the model-data difference.Symbols refer to the season: December-January-February (filled squares), June-July-August (empty squares) and March-April-May-September-October-November (crosses).(c) Comparison of the simulated δD-ln(q) relationships in JJA in the nearsurface vapor for the different days at Kourovka (green squares), for different latitudes at the longitude of Kourovka (magenta crosses) and for different longitudes at the latitude of Kourovka (blue stars).
Table 9. Characteristic of the daily δD-ln(q) distribution in JJA for the LMDZ simulation, for the different LMDZ-ORCHIDEE sensitivity tests listed in Table 1 and for in situ observations.The "δD values for ln(q) = 2.5" are calculated as 2.5 • a + b where a is the slope (given in the first column) and b is the intercept of the linear regression.This gives an idea of the systematic offsets between the different simulations when plotted in the δD-ln(q) diagram.The uncertainty on the slopes is about 6 ‰.

Model
Slope of δD vs. JJA-mean q Minimum δD value for version ln(q) (‰) (g kg −1 ) q value ln(q) = 2. Rayleigh distillation, are misrepresented, and they lead to biases in δD that are not directly related to biases in q.We could speculate on different hypotheses, such as the origin of water vapor issues or cloud processes.However, these are difficult to test.In contrast, in summer, there is a systematic relationship between model-data differences in δD and those in ln(q).When LMDZ has the largest enrichment bias in δD, LMDZ simulates the largest q (the closest to observations).When LMDZ has the most depleted δD (the closest to observations), LMDZ has the largest dry bias in q.Furthermore, the slope of the relationship between δD and is similar to that between δD and q (Table 8).On days when LMDZ is the driest, the slope is the same as other days.Therefore, the model-data differences in q and in δD are probably due to the same processes as those controlling the daily variability of q and δD in nature, i.e., surface evaporation, BL processes and large-scale advection (Sect.5.1).
We notice that the simulated δD-ln(q) slope for the daily data in summer is the same as for the north-south gradient in δD (Fig. 13c).This suggests two possibilities: (1) the δD-ln(q) slopes are not a good way to discriminate between different processes, or (2) the processes controlling the daily variability in δD are the same as those controlling the northsouth gradient in δD.The slope for the east-west gradient is significantly different from the slope for the north-south gradient (Fig. 13c, blue), which may discard the first possibility.Section 3.3.4suggested that the north-south gradient in δD was controlled by purely atmospheric processes.This would discard issues in the surface evaporation as the cause of the dry bias.This would leave BL processes and large-scale advection.
To summarize, analyzing observations and model-data differences in δD-ln(q) plots reveals some information about the sources of model biases, though it does not allow us to pinpoint them with certainty.The model is systematically too enriched all year long probably due to issues in the composition of the moisture at its source.In winter, the δD bias is modulated on a daily basis by processes that are independent from the Rayleigh distillation.In summer, the δD bias is modulated by the same processes as those explaining the daily variability and the north-south gradient in observations, i.e., BL processes and large-scale advection.

Impact of the representation of the land surface
To check whether the enriched bias in LMDZ could be due to the fact that we neglect fractionation during bare soil evaporation, we compare LMDZ simulations with fractionation during bare soil evaporation activated or disabled.When fractionation during bare soil evaporation is activated, δD increases mainly during the driest days of winter.This is because the driest days are associated with less transpiration and more bare soil evaporation into air masses on their way to Kourovka.Since the bare soil evaporation has a lower δD (Risi et al., 2013), more bare soil evaporation leads to a lower δD.Compared to the data, this is not an improvement since the simulated δD-ln(q) slope in winter is already too large in winter.
To qualitatively assess the impact of the representation of the land surface, we compare the sensitivity tests with LMDZ-ORCHIDEE.We find that the representation of the land surface has an important impact on the simulation of q at Kourovka in summer (Table 9).The minimum daily q value varying from 4.2 to 6.0 g kg −1 depending on land surface parameters.This suggests that the land surface representation is important to properly simulate the Siberian climate.
In contrast, the δD-ln(q) slopes are very robust across all sensitivity tests, within the uncertainty.Again, there are two possible interpretations to this robustness in slopes: (1) the δD-ln(q) diagram is not sufficient to distinguish between different processes; (2) the intensity of the moistening and dehydrating processes vary across the different simulations, but keep the same proportions, so that the slopes are the same.Identifying the correct interpretation would require a detailed analysis of the model tendencies for the LMDZ-ORCHIDEE sensitivity tests for which all daily outputs are not available.
The different representations of the land surface result in δD shifts of less 10 ‰.This is much smaller than the modeldata difference of about 30 ‰.Therefore, we suggest that purely atmospheric processes are to cause for the enriched bias of LMDZ.

Conclusions
In this paper, we evaluate the humidity and isotopic composition of water vapor and precipitation simulated by the LMDZ GCM over Siberia using several data sets: TES and GOSAT satellite observations of tropospheric water vapor, GNIP and SNIP networks of precipitation, and daily, in situ measurements of water vapor and precipitation at the Kourovka site in Western Siberia.
LMDZ captures the spatial, seasonal and daily variations reasonably well, except for a few features.LMDZ overestimates the δD in precipitation compared to the precipitation networks, by about 20 ‰.Consistently with this result, LMDZ overestimates the δD observed in both the vapor and precipitation at Kourovka.Based on the analysis of δD-Q diagrams, this bias is most likely associated with a problem with the composition of the vapor at the moisture source.LMDZ slightly underestimates the latitudinal gradient in δD compared to satellite data sets.LMDZ also underestimates the seasonality at Kourovka compared to satellite data sets, but not at the surface compared to the in situ data.Finally, LMDZ captures some aspects of the spatial and daily variations in d excess.
The performance of LMDZ in capturing the spatial and seasonal variations in δD is consistent with other state-ofthe-art models participating in the SWING2 intercomparison project.There is a large spread in the simulation of the latitudinal and seasonal variations in δD, which is not explained by the spread in the simulation of q.This confirms the added value of deuterium measurements compared to q measurements only, though we still need to make progress to better understand why there are systematic biases in δD in models and why different models simulate different δD variations.
Using LMDZ to investigate the processes controlling the daily variability, we find that two processes dominate.First, the large-scale advection determines the origin of the moisture, enriched in the south and west and depleted in the north and east.Second, surface evaporation and BL processes determine the proportion of the vapor coming from enriched surface evaporation or from the more depleted free troposphere.
Most GCMs suffer from a warm and dry bias over midlatitude continents in summer.LMDZ shares this same bias.In addition to a systematic bias towards too enriched δD values, LMDZ exhibits the strongest dry bias on days when it simulates the most depleted δD.Moreover, the slope of δD biases vs. ln(q) biases is consistent with the observed and simulated δD-ln(q) slopes at the daily scale.The slope is also consistent with the observed and simulated δD-ln(q) slopes for latitudinal variations and inconsistent with those for longitudinal variations.This suggests that the same processes that explain the joint δD-q variability at the daily scale and for latitudinal variations explain the biases in q.We suggest that the cause of the dry bias could be either a problem in the large-scale advection bringing too much dry and warm air from the south, or excessive BL mixing.
This paper shows the potential of combining δD and q measurements in δD-q or δD-ln(q) diagrams, to interpret model-data and model-model differences in terms of the representation of physical processes.However, even using such diagrams, it is difficult to discriminate for sure between Rayleigh lines and mixing lines.In addition, different kinds of δD-ln(q) regressions may have the same slope: daily data in summer, daily data in summer for the driest days, and latitudinal gradient in summer.Observations, LMDZ simulations and different sensitivity tests on the land surface may feature large differences in q but still share the same δDln(q) slope.The similarity between different observational, simulated and theoretical slopes leaves us with some ambiguity that we cannot resolve here.Are the δD-ln(q) diagrams a sufficient framework to interpret isotopic data, or do we need to develop a more sophisticated theoretical framework?Or does the similarity between the different slopes shows that the same processes are at play in the observations, in all simulations and for both daily and latitudinal variations?
To resolve this issue, more work would be needed to develop new theoretical frameworks, or to better constrain the existing δD-ln(q) framework.A better spatial coverage of water vapor isotopic measurements would be useful for documenting the composition of the water vapor from the different air masses that are being mixed.Measuring vertical gradient in isotopic composition would be useful for documenting the composition of the evapotranspiration.Furthermore, since the summer dry bias develops as air masses move over the Korouvka region towards the southeast, daily data over stations along a southeast-northwest transect would be useful to document how the dry bias is being formed during synoptic events.Finally, more precise measurement and a better understanding of d excess could add an additional constraint to the system of isotopic equations that has so far suffered from too many unknowns compared to the number of equations.

Fig. 2 .
Fig. 2. Annual mean maps of observed (left) and simulated (right) precipitable water (a-d) and δD in the total column water vapor (e-g) as observed by GOSAT (a-b, e-f), and by TES (c-d, g-h).The black dot indicates the location of Kourovka.Here we plot δD values without any subtraction because LMDZ happens to show values similar to those of GOSAT and TES.The black dot indicates the location of Kourovka.

Figure 2 .
Figure 2. Annual mean maps of observed (left) and simulated (right) precipitable water (a-d) and δD in the total column water vapor (e-g) as observed by GOSAT (a-b, e-f), and by TES (c-d, g-h).The black dot indicates the location of Kourovka.Here we plot δD values without any subtraction because LMDZ happens to show values similar to those of GOSAT and TES.The black dot indicates the location of Kourovka.

Fig. 3 .
Fig. 3. Annual mean maps of observed (left) and simulated (right) δD in the surface precipitation (a-b) and d-excess in the surface precipitation (c-d).The black dot indicates the location of Kourovka.

Figure 3 .
Figure 3. Annual mean maps of observed (left) and simulated (right) δD in the surface precipitation (a-b) and d excess in the surface precipitation (c-d).The black dot indicates the location of Kourovka.

Fig. 4 .
Fig. 4. Annual mean north-south transects of observed (black) and simulated (green) precipitable water observed by GOSAT (a) and by TES (b), δD in the total column water vapor as observed by GOSAT (c) and by TES (d), δD in the surface precipitation (e) and d-excess in the surface precipitation (f).Transects are taken around the longitude of Kourovka: 50 • E-70 • E. For satellite datasets, the annual-mean δD averaged over 50 • E-70 • E-30 • N-80 • N is subtracted to focus on the latitudinal variations.LMDZ outputs are collocated with each dataset and averaging kernels are applied to compare with satellite datasets.The dashed line indicates the latitude of Kourovka.

Figure 4 .
Figure 4. Annual mean north-south transects of observed (black) and simulated (green) precipitable water observed by GOSAT (a) and by TES (b), δD in the total column water vapor as observed by GOSAT (c) and by TES (d), δD in the surface precipitation (e) and d excess in the surface precipitation (f).Transects are taken around the longitude of Kourovka: 50-70 • E. For satellite data sets, the annual-mean δD averaged over 50-70 • E-30-80 • N is subtracted to focus on the latitudinal variations.LMDZ outputs are collocated with each data set and averaging kernels are applied to compare with satellite data sets.The dashed line indicates the latitude of Kourovka.

Fig. 5 .
Fig. 5. Annual mean east-west transects of observed (black) and simulated (green) precipitable water observed by GOSAT (a) and by TES (b), δD in the total column water vapor as observed by GOSAT (c) and by TES (d), δD in the surface precipitation (e) and d-excess in the surface precipitation (f).Transects are taken around the latitude of Kourovka: 50 • N-64 • N.For satellite datasets, the annual-mean δD averaged over 20 • E-120 • E-50 • N-64 • N is subtracted to focus on the longitudinal variations.LMDZ outputs are collocated with each dataset and averaging kernels are applied to compare with satellite datasets.The dashed line indicates the longitude of Kourovka.

Figure 5 .
Figure 5. Annual mean east-west transects of observed (black) and simulated (green) precipitable water observed by GOSAT (a) and by TES (b), δD in the total column water vapor as observed by GOSAT (c) and by TES (d), δD in the surface precipitation (e) and d excess in the surface precipitation (f).Transects are taken around the latitude of Kourovka: 50-64 • N.For satellite data sets, the annual-mean δD averaged over 20-120 • E-50-64 • N is subtracted to focus on the longitudinal variations.LMDZ outputs are collocated with each data set and averaging kernels are applied to compare with satellite data sets.The dashed line indicates the longitude of Kourovka.
APR MAY JUN JUL AUG SEP OCT NOV DEC JAN FEB MAR APR MAY JUN JUL AUG SEP OCT NOV APR MAY JUN JUL AUG SEP OCT NOV DEC

Figure 6 .
Figure 6.Seasonal cycles of observed (black) and simulated (green) precipitable water observed by GOSAT (a) and by TES (b), δD in the total column water vapor as observed by GOSAT (c) and by TES (d), δD (e) and d excess (f) in the surface precipitation.We average over the region of Kourovka (50-64 • N-50-70 • E).For satellite data sets, the annual-mean δD is subtracted to focus on the seasonal variations.For GNIP/SNIP data sets, the spatial standard deviation is also plotted as an envelope.LMDZ outputs are collocated with each data set and averaging kernels are applied to compare with satellite data sets.

Fig. 7 .
Fig. 7. JJA-mean model-data difference for precipitable water (a-b) and total-column water vapor δD (c-d) compared to GOSAT (left) and TES (right).LMDZ outputs are collocated with each dataset and averaging kernels are applied to compare with satellite datasets.The black dot indicates the location of Kourovka.Superimposed are the simulated surface winds.They are the same for all figures.

Figure 7 .
Figure 7. June-July-August (JJA) mean model-data difference for precipitable water (a-b) and total-column water vapor δD (c-d) compared to GOSAT (left) and TES (right).LMDZ outputs are collocated with each data set and averaging kernels are applied to compare with satellite data sets.The black dot indicates the location of Kourovka.Superimposed are the simulated surface winds.They are the same for all figures.

Fig. 8 .
Fig. 8. (a) Annual mean north-south transects of water vapor δD at 600 hPa simulated by different SWING2 models.Transects are taken around the longitude of Kourovka: 50 • E-70 • E. The annual-mean δD averaged over 50 • E-70 • E-30 • N-80 • N is subtracted to focus on the latitudinal variations.The spatial averaging are the same as in Fig. 3. (b) Water vapor δDD at 600 hPa as a function of humidity at 600 hPa.(c) Seasonal amplitudes (June-July-August minus December-January-February) of δD at 600 hPa in the region of Kourovka (50 • N-64 • N-50 • E-70 • E), as a function of the seasonal amplitude in humidity at 600 hPa.

Figure 8 .
Figure 8.(a) Annual mean north-south transects of water vapor δD at 600 hPa simulated by different SWING2 models.Transects are taken around the longitude of Kourovka: 50-70 • E. The annual-mean δD averaged over 50-70 • E-30-80 • N is subtracted to focus on the latitudinal variations.The spatial averaging are the same as in Fig. 3. (b) Water vapor δDD at 600 hPa as a function of humidity at 600 hPa.(c) Seasonal amplitudes (June-July-August minus December-January-February) of δD at 600 hPa in the region of Kourovka (50-64 • N-50-70 • E), as a function of the seasonal amplitude in humidity at 600 hPa.

Figure 9 .
Figure 9. Daily mean time series of humidity (a), water vapor δD (b) and d excess (c) for observations in Kourovka (red) and simulated by LMDZ (blue).

Fig. 10 .
Fig. 10.Daily mean difference between precipitation and water vapor δD (a) and d-excess (b) for observations in Kourovka (red) and LMDZ (blue).The green line represents a theoretical estimation of this difference.

Figure 10 .
Figure 10.Daily mean difference between precipitation and water vapor δD (a) and d excess (b) for observations in Kourovka (red) and LMDZ (blue).The green line represents a theoretical estimation of this difference.

Fig. 11 .
Fig. 11.Schematic representation of the hydrological cycle which includes physical process: advection, deep convection, large-scale condensation, re-evaporation, surface evaporation and boundary-layer processes.

Figure 11 .
Figure 11.Schematic representation of the hydrological cycle which includes physical process: advection, deep convection, large-scale condensation, re-evaporation, surface evaporation and boundary-layer processes.

Fig. 12 .
Fig. 12. Temporal derivative (black line) of (a) humidity and (b) δD at the surface simulated by LMDZ.This temporal derivative can be decomposed into the tendencies from different processes (illustrated in Fig.11): horizontal advection (red); deep convection (green); large-scale condensation and reevaporation (blue); surface evaporation and boundary-layer dynamics (cyan).We used a 5-day filter to focus on seasonal and synoptic variations.

Figure 12 .
Figure 12.Temporal derivative (black line) of (a) humidity and (b) δD at the surface simulated by LMDZ.This temporal derivative can be decomposed into the tendencies from different processes (illustrated in Fig.11): horizontal advection (red); deep convection (green); large-scale condensation and re-evaporation (blue); surface evaporation and boundary-layer dynamics (cyan).We used a 5-day filter to focus on seasonal and synoptic variations.
Fig.13.(a) δD as a function of ln(q), in the observed data (red) and in the model (green).(b) ∆δD as a function of ∆ln(q).∆ refers to the model-data difference.Symbols refer to the season: December-January-February (filled squares), June-July-August (empty squares) and March-April-May-September-October-November (crosses).(c)Comparison of the simulated δD-ln(q) relationships in JJA in the nearsurface vapor for the different days at Kourovka (green squares), for different latitudes at the longitude of Kourovka (magenta crosses) and for different longitudes at the latitude of Kourovka (blue stars).

Figure 13 .
Figure13.(a) δD as a function of ln(q), in the observed data (red) and in the model (green).(b) δD as a function of ln(q).refers to the model-data difference.Symbols refer to the season: December-January-February (filled squares), June-July-August (empty squares) and March-April-May-September-October-November (crosses).(c) Comparison of the simulated δD-ln(q) relationships in JJA in the nearsurface vapor for the different days at Kourovka (green squares), for different latitudes at the longitude of Kourovka (magenta crosses) and for different longitudes at the latitude of Kourovka (blue stars).

Table 1 .
List of sensitivity tests using the LMDZ-ORCHIDEE model.

Table 4 .
Type of relationship Data set δD − Q properties in the data δD − Q properties in LMDZ Characteristics of the spatial distribution in q, δD and in d excess of near-surface vapor, for the LMDZ simulation and for the different LMDZ-ORCHIDEE sensitivity tests listed in Table1.At Kourovka, we calculate an average in the 52-62 • , 55-65 • E domain.The latitudinal gradient is quantified as the difference between the average in the 30-40 • N, 55-65 • E domain and in the 70-80 • N, 55-65 • E domain.The longitudinal gradient is quantified as the difference between the average in the 52-62 • N, 20-30 • E domain and in the 52-62 • N, 110-120 • E domain.

Table 6 .
Daily correlation coefficient (r) and its significance (p value) between temporal derivatives of q or δD and tendencies from different processes simulated by LMDZ: horizontal advection (adv), deep convection (conv), large-scale condensation and re-evaporation (cond), surface evaporation and boundary-layer dynamics (evap).N/S means "not significant".

Table 8 .
Relationships between δD and q and between δD and ln(q).The sign refers to the model-data difference.N/S means "not significant".