The Water Isotopic Version of the Land-Surface Model ORCHIDEE: Implementation, Evaluation, Sensitivity to Hydrological Parameters

Land-Surface Models (LSMs) exhibit large spread and uncertainties in the way they partition precipitation into surface runoff, drainage, transpiration and bare soil evaporation. To explore to what extent water isotope measurements could help evaluate the simulation of the soil water budget in LSMs, water stable isotopes have been implemented in the ORCHIDEE (ORganizing Carbon and Hydrology In Dynamic EcosystEms: the land-surface model) LSM. This article presents this implementation and the evaluation of simulations both in a stand-alone mode and coupled with an atmospheric general circulation model. ORCHIDEE simulates reasonably well the isotopic composition of soil, stem and leaf water compared to local observations at ten measurement sites. When coupled to LMDZ (Laboratoire de Météorologie Dynamique-Zoom: the atmospheric model), it simulates well the isotopic composition of precipitation and river water compared to global observations. Sensitivity tests to LSM (Land-Surface Model) parameters are performed to identify processes whose representation by LSMs could be better evaluated using water isotopic measurements. We find that measured vertical variations in soil water isotopes could help evaluate the representation of infiltration pathways by multi-layer soil models. Measured water isotopes in rivers could help calibrate the partitioning of total runoff into surface runoff and drainage and the residence time scales in underground reservoirs. Finally, co-located isotope measurements in precipitation, vapor and soil water could help estimate the partitioning of infiltrating precipitation into bare soil evaporation.


Introduction
Land-surface models (LSMs) used in climate models exhibit a large spread in the way they partition radiative energy into sensible and latent heat [1,2] precipitation into evapo-transpiration and runoff [3][4][5], evapo-transpiration into transpiration and bare soil evaporation [6,7], and runoff into surface runoff and drainage [8][9][10]. This results in an large spread in the predicted response of surface temperature [11] and hydrological cycle [12,13] to climate change [11] or land use change [14,15]. Therefore, evaluating the accuracy of the partitioning of precipitation into surface runoff, drainage, transpiration and bare soil evaporation (hereafter called the soil water budget) in LSMs is crucial to improve our ability to predict future hydrological and climatic changes.
The evaluation of LSMs is hampered by the difficulty to measure over large areas the different terms of the soil water budget, notably the evapo-transpiration terms and the soil moisture storage [16,17]. Single point measurements of evapo-transpiration fluxes [18] and soil moisture [19] are routinely performed within international networks, but those measurements remain difficult to upscale to a climate model grid box due to the strong horizontal heterogeneity of the land surface [20,21]. Spatially-integrated data such as river runoff observations are very valuable to evaluate soil water budgets at the regional scale [22,23], but are insufficient to constrain the different terms of the water budget. Additional observations are therefore needed.
In this context, water isotope measurements have been suggested to help constrain the soil water budget [24,25], its variations with climate or land use change [26], and its representation by large-scale models [27,28]. For example, water stable isotope measurements in the different water pools of the soil-vegetation-atmosphere continuum have been used to quantify the relative contributions of transpiration and bare soil evaporation to evapo-transpiration [29][30][31][32], to infer plant source water depth [33], to assess the mass balance of lakes [34][35][36] or to investigate pathways from precipitation to river discharge [37][38][39][40]. These isotope-based techniques generally require high frequency isotope measurements and are best suitable for intensive field campaigns at the local scale. At larger spatial and temporal scales, some , where R sample and R SMOW are the isotopic ratios of the sample and of the Vienna Standard Mean Ocean Water (V-SMOW) respectively [50,51]. To first order, variations in δD are similar to those in δ 18 O but are 8 times larger. Deviation from this behavior can be associated with kinetic fractionation and is quantified by deuterium excess (d=δD-8.δ 18 O [50,52]). Hereafter, we note δ 18 O p , δ 18 O v , δ 18 O s , δ 18 O stem and δ 18 O river the δ 18 O of the precipitation, atmospheric vapor, soil, stem, river water respectively. The same subscripts apply for d.

The LMDZ model
LMDZ is the atmospheric GCM (General Circulation Model) of the IPSL (Institut Pierre Simon Laplace) climate model [53,54]. We use the LMDZ-version 4 model [49] which was used in the International Panel on CLimate Change's Fourth Assessment Report simulations [55,56]. The resolution is 2.5 ° in latitude, 3.75 ° in longitude and 19 vertical levels. Each grid cell is divided into four sub-surfaces: ocean, land ice, sea ice and land (treated by ORCHIDEE) (Figure 1a). All parameterizations, including ORCHIDEE, are called every 30 min. The implementation of water stable isotopes is similar to that in other GCMs [57,58] and has been described in [59,60]. LMDZ captures reasonably well the spatial and seasonal variations of the isotopic composition in precipitation [60] and water vapor [61].

The ORCHIDEE (ORganizing Carbon and Hydrology In Dynamic EcosystEms: the land-surface model) model
The ORCHIDEE model is the LSM component of the IPSL climate model. It merges three separate modules: (1) SECHIBA (Schématisation des EChanges Hydriques a l'Interface entre la Biosphère et l'Atmosphère [44,62]) that simulates land-atmosphere water and energy exchanges, (2) STOMATE (Saclay-Toulouse-Orsay Model for the Analysis of Terrestrial Ecosystems [45]) that simulates vegetation phenology and biochemical transfers ; and (3) LPJ (Lund-Postdam-Jena [63]) that simulates the vegetation dynamics. Water stable isotopes were implemented in SECHIBA, and we use prescribed land cover maps so that the two other modules could be de-activated.
Each grid box is divided into up to 13 land cover types: bare soil, tropical broad-leaved ever-green, tropical broad-leaved rain-green, temperate needle-leaf ever-green, temperate broad-leaved ever-green, temperate broad-leaved summer-green, boreal needle-leaf ever-green, boreal broad-leaved summer-green, boreal needle-leaf summer-green, C3 grass, C4 grass, C3 agriculture and C4 agriculture. Water and energy budgets are computed for each land cover type. Figure 1b illustrates how ORCHIDEE (ORganizing Carbon and Hydrology In Dynamic EcosystEms: the land-surface model) represents the surface water budget. Rainfall is partitioned into interception by the canopy and through-fall rain. Through-fall rain, snow melt, dew and frost fill the soil. The soil is represented by two water reservoirs: a superficial and a bottom one [64,65]. Taken together, the two reservoirs have a water holding capacity of 300 mm and a depth of 2 m. Soil water undergoes transpiration by vegetation, bare soil evaporation or runoff. Transpiration and evaporation rates depend on soil moisture to represent water stress in dry conditions. Runoff occurs when the soil water content exceeds the soil holding capacity and is partitioned into 95% drainage and 5% surface runoff [66]. Snowfall fills a singlelayer snow reservoir, where snow undergoes sublimation or melt. By comparison, when not coupled to ORCHIDEE, the simple bucketlike LSM in LMDZ makes no distinction neither between bare soil evaporation and transpiration nor between surface runoff and drainage [67].
Surface runoff and drainage are routed to the coastlines by a water routing model [68]. Surface runoff is stored in a fast ground water reservoir which feeds the stream reservoir with residence time of 3 days. Drainage is stored in a slow ground water reservoir which feeds the stream reservoir with residence time of 25 days. The water in the stream reservoir is routed to the coastlines with a residence time of 0.24 days.

Implementation of water stable isotopes in ORCHIDEE
We represent isotopic processes in a similar fashion as other isotopeenabled LSMs [69][70][71][72][73]. Some details of the isotopic implementation are described in Risi [74]. In absence of fractionation, water stable isotopes ( 16 2 H O , 18 2 H O , HDO, 17 2 H O ) are passively transferred between the different water reservoirs. We assume that surface runoff has the isotopic composition of the rainfall and snow melt that reach the soil surface. Drainage has the isotopic composition of soil water [24]. We calculate the isotopic composition of bare soil evaporation or of evaporation of water intercepted by the canopy using the Craig and Gordon equation [75] (Equation 3). We neglect isotopic fractionation during snow sublimation (Equation 2). We consider isotopic fractionation at the leaf surface (Equation 5) but we assume that transpiration has the isotopic composition of the soil water extracted by the roots (Equation 2).
In the control coupled simulation, we assume that the isotopic composition of soil water is homogeneous vertically and equals the weighted average of the two soil layers. However, transpiration, bare soil evaporation, surface runoff and drainage draw water from different soil water reservoirs whose isotopic composition is distinct [76][77][78]. Therefore, we also implemented a representation of the vertical profile of the soil water isotopic composition.

Stand-alone ORCHIDEE Simulations at MIBA (Moisture in Biosphere and Atmosphere: Network for Water Isotopes in Soil, Stem and Leaf Water) and Carbo-Europe Measurement Sites
First, we performed simulations using ORCHIDEE as a stand-alone model at ten sites. Using isotopic measurements in soil, stem and leaf water, simulations are evaluated at each site at the monthly scale. Sensitivity tests to evapo-transpiration partitioning and soil infiltration processes are performed.

Measurements used for evaluation
To first order the composition of all land surface water pools is driven by that in the precipitation [79]. Therefore, a rigorous evaluation of an isotope-enabled LSM requires to evaluate the difference between the composition in each water pool and that in the precipitation. Besides, to better isolate isotopic biases, we need a realistic atmospheric forcing. We tried to select sites where (1) isotope were measured in different water pools of the soil-plant-atmosphere continuum, during at least a full seasonal cycle and (2) meteorological variables were monitored at a frequency high enough (30 minutes) to ensure robust forcing for our model and (3) water vapor and precipitation were monitored to provide isotopic forcing for the LSM. Only two sites satisfy these conditions: Le Bray and Yatir. Relaxing some of these conditions, we got a more a representative set of ten sites representing diverse climate conditions (Table 1 and Figure 2).
Description of the ten sites: The ten sites belong to two kinds of observational networks: MIBA (Moisture Isotopes in the Biosphere and Atmosphere [80][81][82] or Carbo-Europe [83,84]. Le Bray site, in South-Eastern France, joined the MIBA and GNIP (Global Network for Isotopes in Precipitation) network in 2007. It is an even-aged Maritime pine forest with C3 grass understory that has been the subject of many eco-physiological studies since 1994, notably as part of the Carbo-Europe flux network [85]. In 2007 and 2008, samples in precipitation, soil surface, needles, twigs and atmospheric vapor were collected every month and analyzed for δ 18 O following the MIBA protocol [82,86]. This site was also the subject of intensive campaigns where soil water isotope profiles were collected between 1993 and 1997, and in 2007 [87].
The Yatir site, in Israel, is a semi-arid Aleppo pine forest. It is an afforestation growing on the edge of the desert, with mean-annual precipitation of 280 mm [88,89]. It has also been the subject of many eco-physiological studies as part of the Carbo-Europe flux network [89] and joined the MIBA network in 2004. It. In 2004-2005, samples of soil water at different depth, stems and needles were collected following the MIBA protocol. The water vapor isotopic composition has been monitored daily at the nearby Rehovot site (31.9 ° N, 34.65 ° E, [90]) and is used to construct the water vapor isotopic composition forcing. We must keep in mind however that although only 66 km from Yatir, Rehovot is much closer to the sea and is more humid than Yatir. The precipitation isotopic composition has been monitored monthly at  [93].
Isotopic measurements: Samples of soil water, stems and leaves were collected at the monthly scale. The MIBA and MIBA-US protocols recommend sampling the first 5-10 cm excluding litter and the Carbo-Europe protocol recommends sampling the first 5 cm [84], but in practice the soil water sampling depth varies from site to site. At some sites, soil water was sampled down to 1 m. For evaluating the seasonal evolution of soil water δ 18 18 O in soil, stem and leaf water vary from 0.1‰ to 0.2‰ depending on the sites and involved stable isotope laboratory (Table 2).
Meteorological, turbulent fluxes and soil moisture measurements: At most of the sites, meteorological parameters (radiation, air temperature and humidity, soil temperature and moisture) are continuously measured and are used to construct the meteorological forcing for ORCHIDEE.
Fluxes of latent and sensible energy are measured using the eddy co-variance technique and are used for evaluating the hydrological simulation. Gaps are filled using ERA-Interim reanalyses [94]. Soil moisture observations are available at most sites.

Simulation set-up
To evaluate in detail the isotope composition of different water pools, stand-alone ORCHIDEE simulations on the ten MIBA and Carbo-Europe sites were performed. We prescribe the vegetation type and properties and the bare soil fraction based on local knowledge at each site (Table 3).
ORCHIDEE offline simulations require as forcing several meteorological variables: near-surface temperature, humidity and winds, surface pressure, precipitation, downward longwave and shortwave radiation fluxes. At Le Bray and Yatir, we use local meteorological measurements available at hourly time scale. At other Figure 2: Location of the ten stations used in this study for single-point model-data comparison. The background represents the annual-mean precipitation from GPCP (Global Precipitation Climatology Project) to illustrate the diversity of climate regimes covered by the ten stations. Each station is described in more detail in Table 1 are not available on the other sites. Therefore, we create isotopic forcing using isotopic measurements in the precipitation performed on nearby GNIP or USNIP stations. To interpolate between the nearby stations, we take into account spatial gradients and altitude effects by exploiting outputs from an LMDZ simulation.

Simulated isotopic composition in soil, stem and leaf water:
The soil profile option is activated in all our stand-alone ORCHIDEE simulations. We compare the soil water samples collected in the first 15 cm of the soil (in the first 5-10 cm at many sites) to the soil water composition simulated in the uppermost layer.
The observed composition of stem water is compared to the simulated composition of the transpiration flux.
When comparing observed and simulated composition of leaf water, the Peclet effect, which mixes stomatal water with xylem water (Equation 8), is deactivated. Neglecting the Peclet effect may lead to overestimate of δ 18 O leaf values. sites, we use local meteorological measurements when available and combine them with ERA-Interim reanalyses at 6-hourly time scale for missing variables. At other sites, no nearby meteorological measurements are available and only ERA-Interim reanalyses [94] are used (Table 3).
At each site, we run the model three times over the first year of isotopic measurement (e.g., 2007 at Le Bray). These three years are discarded as spin-up. Then we run the model over the full period of isotopic measurements (e.g., 2007-2008 at Le Bray). We checked that at all sites, the seasonal distribution of δ 18 O, which is the slowest variable to spin-up, is identical between the last year of spin-up and the following year.
We force ORCHIDEE with monthly isotopic composition of precipitation and near-surface water vapor. Since we evaluate the results at the monthly time scale, we assume that monthly isotopic forcing is sufficient. At Le Bray and Yatir, monthly observations of isotopic composition of precipitation and near-surface water vapor are available to construct the forcing. Unfortunately, these observations   Table 3: Information on the offline simulations performed on the 10 sites listed in Table 1: meteorological forcing (6 hourly observations of temperature, humidity, winds, precipitation and radiative fluxes), isotopic forcing (monthly isotopic composition of the precipitation and near-surface water vapor), and prescribed vegetation type and LAI (leaf area index) properties. We give proportions (in %) of the total vegetated area, excluding bare soil. For example, if a given vegetation type covers 100% of the vegetated area and the bare soil fraction is 30%, then the vegetation type covers only 70% of the total area. Three kinds of meteorological forcing are possible: meteorological observations only (obs), meteorological observations filled with ERA-Interim for missing variables (obs_ERA) or ERA-Interim (ERA). Two kinds of isotopic forcing are possible: isotopic composition of precipitation and water vapor observed on the site (obs_iso), or interpolation between GNIP, USNIP or Carbo-Europe stations using the LMDZ atmospheric general circulation model. In the former case, the datasets used for prescribing the water vapor and precipitation isotopic composition forcing are mentionned. In the latter case, GNIP, USNIP or Carbo-Europe stations used to construct the interpolated precipitation isotopic composition forcing are listed. Impact of the temporal sampling: Over the ten sites, samples were collected during specific days and hours. This temporal sampling may induce artifacts when comparing observations to monthly-mean simulated ORCHIDEE values. For soil and stem water, the effect of temporal sampling can be neglected because simulated soil and stem water composition vary at a very low frequency. For leaf water however, there are large diurnal variations [95]. For example, if leaf water is sampled every day at noon when δ 18 O leaf is maximum, then observed δ 18 O leaf will be more enriched than monthly-mean δ 18 O leaf . The exact sampling time is available for Le Bray site only, where we will estimate the effect of temporal sampling.

Spatial heterogeneities:
We are aware of the scale mismatch between punctual in-situ measurements and an LSM designed for large scales (a typical GCM grid box is more than 100 km wide). However, for soil moisture it has been shown that local measurements represent a combination of small scale (10-100 m) variability [20,21] and a large-scale (100-1000 km) signal [96] that a large-scale model should capture [97]. The sampling protocol allows us to evaluate the spatial heterogeneities. For example at Le Bray, two samples were systematically taken a few meters apart, allowing us to calculate the difference between these two samples. On average over all months, the difference between the two samples is 3.5‰ for δ 18 18 O leaf . These error bars need to be kept in mind when assessing modeldata agreement.
Soil moisture: Soil moisture have a different physical meaning in observations and model. Soil moisture is measured as volumetric soil water content (SWC) and expressed in %. In ORCHIDEE, the soil moisture is expressed in mm and cannot be easily converted to volumetric soil water content: the maximum soil water holding capacity of 300 mm and soil depth of 2 m are arbitrary choices and do not reflect realistic values at all sites. In LSMs, soil moisture is more an index than an actual soil moisture content [3]. In this version of ORCHIDEE in particular, it is an index to compute soil water stress, but it was not meant to be compared with soil water content measurements. Therefore, to compare soil moisture between model and observations, we normalize values to ensure that they remains between 0 and 1. The observed normalized SWC is calculated as

Evaluation at measurement sites
In this section, we evaluate the simulated isotopic composition in different water reservoirs of the soil-vegetation-atmosphere continuum at the seasonal scale.
Hydrological simulation: Before evaluating the isotopic composition of the different water reservoirs, we check whether the simulations are reasonable from a hydrological point of view. ORCHIDEE captures reasonably well the magnitude and seasonality of the latent and sensible heat fluxes at most sites (Figures 3 and 4). At Le Bray for example, the correlation between monthly values of evapotranspiration is 0.98 and simulated and observed annual mean evapotranspiration rates are 2.4 mm/d and 2.0 mm/d respectively. However, the model tends to overestimate the latent heat flux at the expense of the sensible heat flux at several sites. This is especially the case at the dry sites Mitra and Yatir: the observed evapo-transpiration is at its maximum in spring and then declines in summer due to soil water stress. ORCHIDEE underestimates the effect of soil water stress on evapo-transpiration and maintains the evapo-transpiration too strong throughout the summer.
The soil moisture seasonality is very well simulated at all sites where data is available (Figures 3 and 4), except for a two-month offset at Yatir (Figure 3f).

Water isotopes in the soil water:
The evaluation of the isotopic composition of soil water is crucial before using ORCHIDEE to investigate the sensitivity to the evapo-transpiration partitioning or to infiltration processes, or in the future to simulate the isotopic composition of paleo-proxies such as speleothems [98].
In observations, at all sites, δ 18  However, the sites with a spring-summer enrichment in ORCHIDEE are not necessarily those with a spring-summer enrichment in observations. This means that ORCHIDEE misses what controls the inter-site variations in the amplitude of the δ 18 O s seasonality. The seasonality is not well simulated at Yatir. This could be due to the missed seasonality in soil moisture and evapo-transpiration. This could be due also to the fact that at Yatir ORCHIDEE underestimates the proportion of bare soil evaporation to total evapo-transpiration: less than 10% in ORCHIDEE versus 38% observed [89], which could explain why the spring enrichment is underestimated. Besides, ORCHIDEE does not represent the diffusion of water vapor in the soil, which could explain why the observed δ 18 O s decrease at Yatir in fall is missed.
When comparing the different sites, annual-mean δ 18 O s follows annual-mean δ 18 O p , with an inter-site correlation of 0.99 in observations. Therefore, it is easy for ORCHIDEE to capture the intersite variations in annual-mean δ 18 O s . A more stringent test is whether ORCHIDEE is able to capture the inter-site variations in annual-mean δ 18 O s -δ 18 O p . This is the case, with a correlation of 0.85 ( Figure 5a) between ORCHIDEE and observations. In ORCHIDEE (and probably in observations), spatial variations in δ 18 O s -δ 18 O p are associated with the relative importance of bare soil evaporation.
Water isotopes in the stem water: In observations, observed δ 18 O stem exhibits no seasonal variations distinguishable from month-tomonth noise (Figures 3 and 4). At Le Bray, Yatir, Mitra, Brloh, Hainich, observed δ 18 O stem is more depleted than the surface soil water. It likely  At Bily Kriz, observed δ 18 O stem is surprisingly more enriched than surface soil water. Several hypotheses could explain this result: (1) the surface soil water could be depleted by dew or frost at this mountainous, foggy site; (2) spruce has shallow roots and therefore sample soil water that is not so depleted; (3) the twigs that were sampled were relatively young so that evaporation from their surface could have occurred when they were still at tree; (4) twigs were sampled in sun-exposed part of the spruce crowns during sunny conditions, which could favor some evaporative enrichment. Additional measurements show a lower Deuterium excess in the stem water compared to the soil water, supporting evaporative enrichment of stems.
ORCHIDEE captures the fact that δ 18 O stem is nearly uniform throughout the year. As for soil water, it is easy for ORCHIDEE to capture the inter-site variations in annual-mean δ 18 (Figure 5b). It is not able to capture δ 18 O stem values that are either more enriched or more depleted than δ 18 O s . This could be due to the fact that ORCHIDEE underestimates vertical variations in soil isotopic composition. Also, ORCHIDEE is not designed to represent deep ground water sources or photosynthesizing twigs.
Vertical profiles of soil water isotope composition: At Le Bray, we compare our offline simulation for 2007 with soil profiles collected from 1993 to 1997 and in 2007 (Figure 6a-6b). The year mismatch adds a source of uncertainty to the comparison. In summer (profiles of August 1993 and September 1997), the data exhibits an isotopic enrichment at the soil surface of about 2.5‰ compared to the soil at 1 m depth (Figure 6a), likely due to surface evaporation [99]. Then, by the end of September 1994, the surface becomes depleted, likely due to the input of depleted rainfall. Previously enriched water remains between 20 and 60 cm below the ground, suggesting an infiltration through piston-flow [100]. ORCHIDEE predicts the summer isotopic enrichment at the surface, but slightly later in the season (maximum in September rather than August) and underestimates it compared to the data (1.5‰ enrichment compared to 2.5‰ observed, Figure 6b). The model also captures the surface depletion observed after the summer, as well as the imprint of the previous summer enrichment at depth. However, ORCHIDEE simulates the surface depletion in December, whereas the surface depletion can be observed sooner in the data, at the end of September 1994.
At Yatir, observed profiles exhibit a strong isotopic enrichment from deep to shallow soil layers in May-June by up to 10‰ ( Figure  6c). As for Le Bray, the model captures but underestimates this isotopic enrichment in spring and summer by about 3‰ (Figure  6d). This discrepancy could be the result of underestimated bare soil evaporation. Observed profiles also feature a depletion at the surface in winter that the model does not reproduce. This depletion could be due to back-diffusion of depleted vapor in dry soils [99,[101][102][103], a process that is not represented in ORCHIDEE but likely to be significant in this region. Soil evaporation fluxes measured with a soil chamber at Yatir shows that when soils are dry, there is adsorption of vapor from the atmosphere to the dry soil pores before sunrise and after sunset [104].
Water isotopes in leaf water: It is important to evaluate the simulation of the isotopic composition of leaf water by ORCHIDEE if we want to use this model in the future for the simulation of paleoclimate proxies such tree-ring cellulose [105,106], for the simulation of the isotopic composition of atmospheric CO 2 which may be used to partition CO 2 fluxes into respiration from vegetation and soil [107,108] or for the simulation of the isotopic composition of atmospheric O 2 which may be used to infer biological productivity [109,110].
In the observations, δ 18 O leaf exhibits a large temporal variability reflecting a response to changes in environmental conditions (e.g., relative humidity and the isotopic composition of atmospheric water vapor). At all sites except at Yatir, δ 18 O leaf is most enriched in summer than in winter, by up to 15‰ (Figures 3 and 4). This is because the evaporative enrichment is maximum in summer due to drier and warmer conditions. ORCHIDEE captures the maximum enrichment in summer. However, ORCHIDEE underestimates the annual-mean δ 18 O leaf at most sites ( Figure 5). This could be due to the fact that most leaf samples were collected during the day, when the evaporative enrichment is at its maximum, while for ORCHIDEE we plot the daily-mean δ 18 At Le Bray, if we sample the simulated δ 18 O leaf during the correct days and hours, simulated δ 18 O leaf increases by 4‰ in winter and by 10‰ in summer. Such an effect can thus quantitatively explain the model-data mismatch. After taking this effect into account, simulated δ 18 O leaf may even become more enriched than observed. This is the case at Le Bray, especially in summer. The overestimation of summer δ 18 O leaf could be due to neglecting diffusion in leaves or non-steady state effects.
Again, Yatir is a particular case. Minimum δ 18 O leaf occurs in spring-summer while the soil evaporative enrichment is maximum. In arid regions and seasons, leaves may close stomata during the most stressful periods of the day, inhibiting transpiration, and thus retain the depleted isotopic signal associated with the moister conditions of the morning [111,112]. ORCHIDEE does not represent this process and thus simulates too enriched δ 18 O leaf .

Summary:
Overall, ORCHIDEE is able to reproduce the main features of the seasonal and vertical variations in soil water isotope content, and seasonal variations in stem and leaf water content. Discrepancies can be explained by some sampling protocols, by shortcomings in the hydrological simulation or by neglected processes in ORCHIDEE (e.g., fractionation in the vapor phase).
The strong spatial heterogeneity of the land surface at small scales does not prevent ORCHIDEE from performing reasonably well. This suggests that in spite of some small-scale spatial heterogeneities at each site, local isotope measurements contain large-scale information and are relevant for the evaluation of large-scale LSMs.

Sensitivity analysis
Sensitivity to evapo-transpiration partitioning: Several studies have attempted to partition evapo-transpiration into the transpiration and bare soil evaporation terms at the local scale [29][30][31]113]. Estimating E/ET, where E is the bare soil evaporation and ET is the evapo-transpiration, requires measuring the isotopic composition of soil water, stem water and of the evapo-transpiration flux. The isotopic composition of the evapo-transpiration can be estimated through "Keeling plots" approach [114], but this is costly [29] and the assumptions underlying this approach are not always valid [115].
Considering a simple soil water budget at steady state and with vertically-uniform isotopic distribution, we show that although estimating E/ET requires measuring the isotopic composition of the where α eq and α k are the equilibrium and kinetic fractionation coefficients respectively.
Below, we show that this equation can apply to annual-mean quantities, neglecting effects associated with daily or monthly co-variations between different variables. We investigate to what extent this equation allows us to estimate the magnitude of E/I at local sites.
At the Yatir site, all the necessary data for equation 1 is available. An independent study has estimated E/I =38% [89]. Using annually averaged observed values (δ 18 O p =-5.1‰ and δ 18 O s =-3.7‰ in the the surface soil), we obtain E/I =46%. However, in ORCHIDEE, the annually averaged surface δ 18 O s is 0.8 lower when sampled at the same days as in the data. When correcting for this bias, we obtain E/I =28%. Observed E/I lies between these two estimates. This shows the applicability of this estimation method, keeping in mind that estimating E/I is the most accurate where E/I is lower.
When we perform sensitivity tests to ORCHIDEE parameters at the various sites, the main factor controlling δ 18 O s is the E/I fraction. This is illustrated as an example at Le Bray and Mitra sites (Figure 7). Sensitivity tests to parameters as diverse as the rooting depth or the stomatal resistance lead to changes in δ 18  To summarize, local observations of δ 18 O s -δ 18 O p could help constrain the simulation of E/I in models. This would be useful since the evapo-transpiration partitioning has a strong impact on how an LSMs represents land-atmosphere interactions [116].

Sensitivity to soil infiltration processes:
Partitioning between evapo-transpiration, surface runoff and drainage depends critically on how precipitation water infiltrates the soil [5,8,10], which is a key uncertainty even in multi-layer soil models where infiltration processes are represented explicitly [62]. It has been suggested that observed isotopic profiles could help understand infiltration processes at the local scale [100]. The capacity of ORCHIDEE to simulate soil profile allows us to investigate whether measured isotope profiles in the soil could help evaluate the representation of these processes also in largescale LSMs.
With this aim, we performed sensitivity tests at Le Bray. The simulated profiles are sensitive to vertical water fluxes in the soil. When the diffusivity of water in the soil column is decreased by a factor 10 from 0.1 to 0.01 compared to the control simulation, the deep soil layer becomes more depleted by about 0.7‰ (Figure 8) and the isotopic gradient from soil bottom to top becomes 30% steeper in summer, because the enriched soil water diffuses slower through the soil column.
Simulated profiles are also sensitive to the way precipitation infiltrates the soil. When precipitation is added only to the top layer (piston-flow infiltration) the summer enrichment is reduced by mixing of the surface soil water with rainfall, and it propagates more easily to lower layers during fall and winter. Conversely, when rainfall is evenly spread throughout the soil column (a crude representation of preferential pathway infiltration), the surface enrichment is slightly more pronounced and the deep soil water is more depleted by up to 0.8‰ in winter ( Figure 8). However, the observed surface depletion occurs in February with preferential pathways, compared to December in the piston-like in infiltration. The quick surface depletion observed after the summer suggests that infiltration is dominated by the pistonlike mechanisms.
To summarize, we show that vertical and seasonal variations of δ 18 O s are very sensitive to infiltration processes, and are a powerful tool to evaluate the representation of these processes in LSMs.

Global-scale Simulations Using the Coupled LMDZ-ORCHIDEE Model Simulation set-up
To compare with global datasets, we performed LMDZ-ORCHIDEE coupled simulations. In all our experiments, LMDZ three-dimensional fields of horizontal winds are nudged towards ECMWF (European Center for Medium range Weather Forecast) reanalyses [117]. This ensures a realistic simulation of the large-scale atmospheric circulation and allows us to perform a day-to-day comparison with field campaign data [60,118]. At each time step, the simulated horizontal wind field u  is relaxed towards the reanalysis following this equation: where obs u  is the reanalysis horizontal wind field, F  is the effect of all simulated dynamical and physical processes on u  , and τ is a time constant set to 1 h in our simulations [119].
To compare with global datasets, LMDZ-ORCHIDEE simulations are performed for the year 2006, chosen arbitrarily. We are not interested in inter-annual variations and focus on signals that are much larger. To ensure that the water balance is closed at the annual scale, we performed iteratively 10 times the year 2006 as spin-up. In these simulations, the Peclet and non-steady state effects are de-activated.
To compare with field campaign observations in 2002 and 2005, we use simulations performed for these specific years, initialized from the 2006 simulation. In these simulations, we test activating or deactivating the Peclet effect.
In all LMDZ-ORCHIDEE simulations, canopy-interception was de-activated (consistent with simulations that our modeling group performed for the Fourth Assessment Report).

Evaluation of water isotopes in leaf water at the diel scale during campaign cases
Daily data from field campaigns: Two field campaigns are used to evaluate the representation of δ 18   Because meteorological and isotopic forcing are not available for the entire year, we prefer to compare these measurements with LMDZ-ORCHIDEE simulations. At both sites, the simulated δ 18 O v and δ 18 O stem are consistent with those observed (model-data mean difference lower than 1.4‰ in Kansas and 0.4‰ at Hartheim), allowing us to focus on the evaluation of leaf processes.
Evaluation results: At the Kansas grassland site, δ 18 O leaf exhibits a diel cycle with an amplitude of about 10‰ [120]. LMDZ-ORCHIDEE captures this diel variability, both in terms of phasing and amplitude ( Figure 9). The model systematically overestimates δ 18 O leaf by about 4‰, in spite of the underestimation of the stem water by 1.4‰ on average. This may be due to a bias in the simulated relative humidity (LMDZ is on average 13% too dry at the surface, which translates into an expected enrichment bias of 3.9‰ on the leaf water assuming steady state based on Equation 7) or to uncertainties in the kinetic fractionation during leaf water evaporation.
At the Hartheim pine plantation, δ 18 O leaf is on average 8‰ more depleted for current-year needles than for 1-year-old needles. Also, the observed diel amplitude is weaker for current-year needles (5 to 8‰) than for 1-year-old needles (10 to 15‰). These observations are consistent with a longer diffusion length for current-year needles (15 cm) than for 1-year-old needles (5 cm) [121] and with a larger transpiration rate, leading to a stronger Peclet effect. When neglecting Peclet and non-steady state effects, ORCHIDEE simulates an average δ 18 O leaf close to that of 1-year-old needles, consistent with the small diffusion length and evaporation rate of these leaves. ORCHIDEE captures the phasing of the diurnal cycle, but underestimates the diel amplitude by about 4‰. This is probably due to the underestimate of the simulated diel amplitude of relative humidity by 20%. Accounting for Peclet and non-steady state effects strongly reduces both the average δ 18 O leaf and its diel amplitude (Figure 9), in closer agreement with current-year needles.
To summarize, ORCHIDEE simulates well the leaf water isotopic composition. The leaf water isotope calculation based on Craig et al. [75] simulates the right phasing and amplitude for leaves that have short diffusive lengths or low transpiration rates. Non-steady state and diffusion effects need to be considered in other cases. By activating or de-activating these effects, ORCHIDEE can simulate all cases.

Evaluation of water isotopes in precipitation
Precipitation datasets: To evaluate the spatial distribution of precipitation isotopic composition simulated by the LMDZ-ORCHIDEE coupled model, we use data from the Global Network for Isotopes in Precipitation (GNIP [122]), further complemented by data from Antarctica [123] and Greenland [124]. We also use this network to construct isotopic forcing at sites where the precipitation was not sampled, complemented with the USNIP (United States Network for Isotopes in Precipitation [125]) network. This good model-data agreement can be obtained even when we deactivate ORCHIDEE. When we use LMDZ in a stand-alone mode, in which the isotope fractionation at the land surface is neglected [60], the model-data agreement is as good as when we use LMDZ-ORCHIDEE. Therefore, fractionating processes at the land surface have a second  [122], Antarctica [123] and Greenland [124] data. The data is gridded over a coarse 7.5 × 6.5° grid for visualization purposes. b) Same as a) but for annual mean d p . c) Annual mean δ 18 O p simulated by coupled LMDZ-ORCHIDEE model for the control simulation. d) same as c) but for annual mean d p . order effect on precipitation isotopic composition, consistent with [28,[71][72][73].
To quantify in more detail the effect of fractionation at the land surface, we performed additional coupled simulations with LMDZ-ORCHIDEE. We compare the control simulation described above (ctrl) to a simulation in which fractionation at the land surface was de-  (Figure 11). In nofrac, the composition of bare soil evaporation equals that of soil water. Even when restricting the analysis to continental regions, the spatial correlations between the ctrl and nofrac simulations are 0.999 and 0.95 for δ 18 O p and d p respectively, and the root mean square differences are 0.27‰ and 1.1‰ for δ 18 O p and d p respectively. This confirms that fractionation at the land surface has a second-order effect on precipitation isotopic composition compared to the strong impact of atmospheric processes.
However, to second order, a detailed representation of fractionation at the land surface lead to a slight improvement in the simulation of δ 18 O p and to a significant improvement in that of d p . In ctrl, δ 18 O p is lower by up to 1.5‰ and d p higher by up to 5‰ than in nofrac over boreal continental regions such as Siberia, Canada and central Asia, consistent with the expected effect of fractionation at surface evaporation [42]. Taking into account fractionation at the land surface leads to a better agreement with the GNIP data over these regions, where δ 18 O p is overestimated by about 4‰ and d p underestimated by 4 to 7‰ when neglecting fractionation at the land surface. The effect of fractionation is maximal over these boreal regions because (1) the fraction of bare soil evaporation is maximal, (2) a significant proportion of evaporativelyenriched soil water is lost by drainage and (3) a larger proportion of the moisture comes from land surface recycling [48, 126,127]. Similar results were obtained with other models [128].
To summarize, LMDZ-ORCHIDEE simulates well the spatial distribution of precipitation isotopic composition, but this distribution is not a very stringent test for the representation of land surface processes in ORCHIDEE. In the next section, we argue that the distribution of river isotopic composition is a more stringent test.
Evaluation of water isotopes in river water: Large rivers integrate a wide range of hydrological processes at the scale of GCM grid boxes [22,23,[129][130][131]. Here we evaluate the isotopic composition of river water simulated by ORCHIDEE using data collected by the Global Network for isotopes in Rivers (GNIR [132,133]).
Observed annual mean δ 18 O river follows to first order the isotopic composition of precipitation [79], and is thus also well simulated by LMDZ-ORCHIDEE (Figure 12a and 12b), with a spatial correlation between measured and simulated δ 18 O river of 0.80 and a RMSE of 3.2‰ over the 149 LMDZ grid boxes containing data. Regionally however, the 18 O δ difference between precipitation and river water (δ 18 O river -δ 18 O p ) can be substantial and provides a stronger constraint for the model.
Over South America, Europe and some parts of the US, the river water is typically 1‰ to 4‰ more depleted than the precipitation (Figure  12a), because precipitation contributes more to rivers during seasons when it is the most depleted [134]. In contrast, over central Asia or northern America, river water is more enriched than precipitation, due to evaporative enrichment of soil water [79,134,135]. This is further confirmed by a simulation where fractionation at the land surface was neglected (not shown), for which the river water is in global average 5‰ more depleted.

Sensitivity to the representation of pathways from precipitation to rivers
At the local scale, water isotopes have already been used to partition river discharge peaks into the contributions from recent rainfall and soil water [37][38][39]. Given the property of rivers to integrate hydrological processes at the basin scales [22,23,[129][130][131], we now explore to what extent δ 18 O river could help evaluate pathways from precipitation to rivers in LSMs. We illustrate this using seasonal variations in δ 18 O river on two well established GNIR and GNIP stations in Vienna (Danube river) and Manaus (the Amazon) ( Figure 13). The seasonal cycle in δ 18 O river is attenuated compared to that in δ 18  To understand better what determines the attenuation and lag of the seasonality in δ 18 O river compared to that in δ 18 O p , we perform sensitivity tests to ORCHIDEE parameters. Parameters tested include the partitioning of excess rainfall into surface runoff and drainage and the residence time scale of different reservoirs (slow, fast and stream) in the routing scheme. River discharge is extremely sensitive to these parameters [136].
If all the runoff occurs as surface runoff (Figure 13), then the seasonal cycle of δ 18 O river is similar to that of δ 18 O p . This shows that the attenuation and lag of the seasonality in δ 18 O river compared to that in δ 18 O p are caused by the storage of water into the slow reservoir, which accumulates drainage water.
When the residence time scale of the slow reservoir is multiplied by 2 (i.e., the water from the slow reservoir is poured twice faster into the streams, Figure 12), the simulated lag of δ 18 O river at Vienna increases from 4 to 5 months (in closer agreement with the data). In contrast, the seasonal cycle in δ 18 O river is not sensitive to residence time scales in the stream and fast reservoirs, which are too short to have any impact at the seasonal scale.
To summarize, ORCHIDEE performs well in simulating the seasonal variations in δ 18 O river . In turn, δ 18 O river observations could help estimate the proportion of surface runoff versus drainage and calibrate empirical residence time constants in the routing scheme, offering a mean to enhance model performance. To test the effect of the assumption that the soil water isotopic composition is vertically constant, we applied Equation 1 using δ 18 O sδ 18 O p from a simulation with soil profiles activated. This assumption is a significant source of uncertainty on estimating E/I (Table 4). We also analyzed the effect of potential measurement errors in δ 18 O s, δ 18 O p , δ 18 O v temperature or relative humidity on the E/I reconstruction. Results are relatively insensitive to small errors in these measurements (Table  4). However, results are sensitive to the choice of the n exponent in the calculation of the kinetic fractionation α k ( Table 4): knowing the n exponent with an accuracy of 0.07 (e.g., estimated n ranges from 0.63 to 0.70) is necessary to estimate E/I with an absolute precision of 2%.
Finally, estimating E/I using equation 1 bears additional sources of uncertainty in that we cannot estimate using the ORCHIDEE model. These are related to all processes that ORCHIDEE does not simulate. For example, ORCHIDEE underestimates or mis-represents the vertical isotopic gradients in soil water at some sites and does not represent the effect of water vapor diffusion in the soil. These effects may disturb the proportionality between E/I and δ 18 O s -δ 18 O p in practical applications.
To summarize, co-located isotope measurements in precipitation, vapor and soil water could provide an accurate constrain on the proportion of bare soil evaporation to precipitation infiltration.

Conclusion and Perspectives
The ORCHIDEE LSM, in which we have implemented water stable isotopes, reproduces the isotopic compositions of the different water pools of the land surface reasonably well compared to local data from MIBA and Carbo-Europe and to global observations from the GNIP and GNIR networks. Despite the scale mismatch between local measurements and a GCM grid box, and despite the strong spatial heterogeneity in the land surface, the capacity of ORCHIDEE to reproduce the seasonal and vertical variations in the soil isotope composition suggests that even local measurements can yield relevant information to evaluate LSMs at the large scale.
We show that the simulated isotope soil profiles are sensitive to infiltration pathways and diffusion rates in the soil. The spatial and seasonal distribution of the isotope composition of rivers is sensitive to the partitioning of total runoff into surface runoff and drainage and to the residence time scales in underground reservoirs. The isotopic composition of soil water is strongly tied to the fraction of infiltrated water that evaporates through the bare soil. These sensitivity tests suggest that isotope measurements, combined with more conventional measurements, could help evaluate the parameterization of infiltration processes, runoff parameterizations and the representation of surface water budgets in LSMs.
Evaluating an isotopic LSM requires co-located observations of the isotope composition in precipitation, vapor and soil at least at the monthly scale. However, such co-located measurements are still very scarce, and most MIBA and Carbo-Europe sites are missing one of the components. Therefore, for LSM evaluation purpose, we advocate for the development of co-located isotope measurements in the different water pools at each site, together with meteorological variables. Our results suggest that isotope measurements are spatially relatively well representative and that even monthly values are already valuable to identify model bias or to estimate soil water budgets. Therefore, in the perspective of LSM evaluation, if a compromise should be made with sampling frequency and spatial coverage, we favor co-located

Evapo-transpiration partitioning
In this section, we generalize at the global scale our results on evapo-transpiration partitioning estimates.
We apply equation 1 to annual-mean outputs from a LMDZ-ORCHIDEE simulation. We compare E/I estimated from Equation 1 to E/I directly simulated by LMDZ-ORCHIDEE. The spatial pattern of E/I is remarkably well estimated by Equation 1 (Figure 14). The equation captures the maximum over the Sahara, Southern South America, Australia, central Asia, Siberia and Northern America. The isotope-derived spatial distribution of E/I correlates well with the simulated distribution (r=0.91). Average errors are lower than 50% of the standard deviation at the global scale. This confirms that covariation between the different variables at sub-annual time scales has measurements of all the different water pools at the monthly scale on a few sites representative of different climatic conditions, rather than multiplying sites where water pools are not all sampled. Additionally, at each observation site, collecting different soil samples a few meters apart is helpful to check that they are spatial representative. In the future, development in laser technology [137,138] will allow the generalization of water vapor isotope monitoring at the different sampling sites, which has long been a very tedious activity [90].
From the modeling point of view, kinetic fractionation processes during bare soil evaporation are a source of uncertainty, and a better understanding and quantification of this fractionation is necessary [103,139]. In addition, the accuracy of isotopic simulations by LSM is expected to improve as the representation of hydrological processes improves. In particular, given the importance of vertical water exchanges for the isotopic simulation, implementing water isotopes in a multi-layer hydrological parameterization with sufficient vertical resolution [69] is crucial. In the future, we plan to implement water isotopes in the latest version of ORCHIDEE, which is multi-layer and more sophisticated [140][141][142]. Finally, latest findings largely based on water isotopic measurements suggest that different water pools coexist within a soil column and that evaporation, transpiration, runoff and drainage tap from these different pools [77,143,144]. These effects are not yet represented explicitly in global LSMs. These effects were mainly evidenced based on isotope measurements, and in turn, their representation expected to significantly impact isotopic simulations. Such feedbacks between isotopic research and hydrological parameterization improvements should lead to LSM improvements    in the future. With this in mind, LSM inter-comparison projects would strongly benefit from including water isotopes as part of their diagnostics, in the lines of iPILSP (isotope counterpart of the Project for Intercomparison of Land-surface Parameterization Schemes [27]).

Representation of isotope fractionation during evaporation from land surface water pools
Processes for which we neglect fractionation: Snow sublimation is associated with a slight fractionation due to exchanges between snow and vapor in snow pores [115,145,146]. However, we assume that these effects are small enough to be neglected, as in other GCMs [58].
Water uptake by roots has been shown to be a non-fractionating process [147,148], but fractionation at the leaf surface during transpiration impacts the composition of transpired fluxes at scales shorter than daily [95,137]. As the application of ORCHIDEE in the context of our study focuses mainly on time scales of a month or longer, we assume here that the transpiration and stem water have the composition of soil water extracted by the roots.

Evaporation from bare soils and canopy-intercepted water:
We represent isotope fractionation during evaporation of soil and canopyintercepted water using the model of Craig [75]: at any time t, the isotopic composition of evaporation R E is given by: where R l and R v are the isotopic compositions of liquid water at the evaporative site and of water vapor respectively, h is the relative humidity normalized to surface temperature, α eq is the isotopic fractionation during liquid-vapor equilibrium [149] and α k is the kinetic fractionation during water vapor diffusion. The kinetic fractionation during soil evaporation is still very uncertain [103,150]. We use the very widespread formulation of [99,151]: where D and D i are the molecular diffusivities of light and heavy water vapor in air, respectively, and n is an exponent that depends on the flow regime (0.5, 0.67 and 1 for turbulent, laminar and stagnant regimes respectively) but remains difficult to estimate [103,150]. In this study, we take n =0.67 for both evaporation of soil and canopyintercepted water, corresponding to moist conditions in the case of soils [99]. However, we also tried 0.5 and 1.0 to estimate the range of uncertainty related to this parameter. The isotopic composition of precipitation is only slightly sensitive to the formulation of the kinetic fractionation: when n varies from 0.5 to 1, significant changes in δ 18 O p and d p are restricted to areas where bare soil covers more than 70%. Even in those case, changes in δ 18 O p and d p never exceed 2‰ and 7‰ respectively. The impact is slightly stronger on soils. Varying n from 0.5 to 1 leads to δ 18 O s variations of 2‰ in offline simulations on the Bray site, of the order of the observed average difference between two samples collected on the same day (2.2‰). In coupled simulations, the impact on δ 18 O s and d s reaches 8‰ and 20‰ respectively on very arid regions such as the Sahara.
To calculate the temporal mean isotopic composition of evaporation over the time step Δt, E R , we assume R v and h are constant throughout each time step. On the other hand, we allow the isotopic ratio of liquid water to vary over the simulation time step Δt following [151]. While assuming constant R l is a valid assumption for models with very short time steps [152], it is not the case in ORCHIDEE (Δt =30min). We then calculate E R as: where R l0 is the initial isotopic ratio of liquid water, f is the remaining liquid fraction in the water reservoir affected by isotopic enrichment, and β and γ are parameters defined by Stewart [151]: For canopy-intercepted water, the water reservoir is sufficiently small to assume that the water reservoir affected by isotopic enrichment is the total canopy-intercepted water. For soil evaporation on the other hand, we assume that the depth of the water reservoir affected by isotopic enrichment equals the average distance traveled by water molecules in the soil: where K D is the effective self-diffusivity of liquid water in the soil column. Neglecting the dispersion term, K D is given by Munnich et al. [100,147,[151][152][153]: where D m =2.5 9 2 10 / m s − ⋅ is the molecular liquid water selfdiffusivity [154,155], τ is the soil tortuosity and l θ is the volumetric soil water content. In the control simulation, we assume θ l .τ =0.1 leading to L =0.67 mm. This choice is consistent with a τ of 0.67 [151] and an average θ l of about 15%. At the Bray, measurements along profiles show θ l varying from about 5 to 30%. Since these values are difficult to constrain observationally and very variable spatially and temporally, sensitivity tests to θ l .τ are performed and described. We neglect the vapor phase in the soil and associated fractionation and diffusion processes [153].

Dew formation:
We assume fractionation during dew and frost formation following a Rayleigh distillation of the vapor in the lowest 10 hPa (~ 80 m) of the atmosphere. Since the atmospheric water vapor condenses in small proportion during frost and dew, this choice of the depth of atmosphere involved in the condensation has almost no impact on the composition of the dew and frost formed. Following common practice, we use equilibrium fractionation coefficient from Merlivat et al. [148,156,157] and the kinetic fractionation formation of [158] with λ =0.004, whose choice has very little impact on the results.
Leaf water evaporation: At isotopic steady state, the composition of water transpired by the vegetation is equal to that of the soil water extracted by the roots. In default simulations, we assume that isotopic steady state for plant water is established at any time and we diagnose the composition of the leaf water at the evaporation site, SS e R , by inverting the Craig and Gordon equation [75]: where R s and R v are the isotopic ratio in soil water and water vapor respectively, h is the relative humidity normalized to surface temperature, α eq is the isotopic fractionation during liquid-vapor equilibrium [148] and α k is the kinetic fractionation during water vapor diffusion. We take the same kinetic fractionation formulation as for the soil evaporation [150], with n =0.67 [31,69]. Leaf water compositions are significantly sensitive to parameter n , with variations of the order of 10‰ as n varies from 0.5 to 1. We assume that the leaf temperature used to calculate α eq is equal to the soil temperature, but results are very little sensitive to this assumption.
The isotopic composition of leaf water has been the subject of many observational and numerical modeling studies [86,[159][160][161]. Several studies have shown that the composition of the leaves is affected by mixing with xylem water and by non-stationary effects [161,162]. Nonsteady state effects are also incorporated in ORCHIDEE following [159]. The isotopic ratio in the leaf mesophyll SS L R is the result of the mixing between leaf water at the evaporative site and xylem water (Peclet effect): where f is a coefficient decreasing as the Peclet effect increases: E is the transpiration rate per leaf area, L eff is the effective diffusion length and W is the leaf water content per leaf volume (assumed equal to 10 3 kg/m 3 , order of magnitude in [121]). The Peclet number P can be tuned by changing L eff , that depends on leaf geometry and drought intensity (e.g., 7 to 12 mm in Cuntz et al. [161], 50 to 150 mm in Barnard et al. [121]). We take L eff =8 mm to optimize our simulation on Hartheim.
For some simulations, we account for the effect of water storage in leaves (leading to some memory in the leaf water isotopic composition) following Dongmann [163]. Assuming that W is constant, we calculate the leaf lamina composition R L as Farquhar [159]: and g is the sum of the total (stomatic and boundary layer) conductances. The isotopic composition of transpiration is then calculated so as to conserve isotope mass.

Representation of the vertical distribution of soil water isotopic composition
Principle: In control simulations, we assume that the isotopic composition of soil water is homogeneous vertically and equals the weighted average of the two soil layers. In addition, to test this assumption, we implemented a representation of the vertical distribution of the soil water isotopic composition: the soil water is spread vertically between several layers. The first layer contains a water height = D L K t ⋅ ∆ , where K D is the diffusivity of water molecules in water and Δt is the time step of the simulation, and the other layers contain a water height resol.L. The parameter resol can be tuned to find a compromise between vertical resolution and computational time. Layers are created from the top to bottom until all layers are full with water except the deepest one that contains the remaining soil water. For example, with L =0.67 mm, up to 16 layers can thus be created if the soil is saturated. Bare soil evaporation is extracted from the first layer. Transpiration is extracted from the different layers following a root extraction profile that reflects the sensitivity of transpiration to soil moisture [164]. Drainage takes water from the deepest layer. In the control simulation, rain and snow melt are added to the first layer (piston-like flow). In a sensitivity test, that can also be homogeneously distributed in the different layers, to crudely represent preferential pathways through fractures or pores in the soil.
At each time step, the soil water isotopic composition in each layer is re-calculated by taking into account the sources and sinks for each layer and ensuring that each layer remains full except the deepest one. Isotopic diffusion between adjacent layers is applied at each time step (Equation 6). The water budget of the total soil remains exactly the same as without vertical discretization.
Evaluation for an idealized case: The module representing vertical distribution of water isotopes in the soil is first evaluated for an idealized case when it is not yet embedded into ORCHIDEE.
First, we use a case in which the soil column evaporates at its top and is permanently refilled at the bottom by a water with δ 18 O of -8‰ [152]. The soil remains saturated, and we focus on the steady state reached after a few hundreds of days [152]. An analytical solution is available for this case [100,165]. The analytical solution and a much more sophisticated model of soil water isotopes (MuSICA [166]) yield very similar results (Figure 15a): the bottom of the soil is at -8‰ while the top of the soil is enriched up to 15‰. The soil module of ORCHIDEE is able to reproduce these results when the value of θ l .τ is set to be very low (0.001) and when the vertical resolution is sufficiently high (layers of 0.75 mm). Whatever the value for θ l .τ, ORCHIDEE results become less sensitive to the vertical discretization when layers are thinner than about 2 mm.
Second, we use a case in which the soil column, initially with a soil water of -8‰, evaporates at its top until the soil water content is only 20% [99]. The atmosphere has a relative humidity of 20% and a vapor δ 18 O of -15‰. The sophisticated models MuSICA and SiSPAT [152] feature a typical evaporative enrichment profile, with δ 18 O increasing from its initial value of -8‰ at the bottom to a maximum δ 18 O of 13‰ about 10 mm below the surface (Figure 15b). In the uppermost 10 mm, there is a slight depletion due to diffusion of water vapor into the soil column [101]. ORCHIDEE is not able to reproduce this vertical profile. First, since diffusion of water vapor in the soil is neglected, it is not able to simulate the depletion near the surface. Second, since θ l .τ is temporally and vertically constant in ORCHIDEE, it is not able to adapt to the drying of the soil. In the sophisticated model, as the soil dries, the soil water content l θ decrease, thus inhibiting vertical mixing of soil water and favoring strong isotopic gradients. In contrast in ORCHIDEE, θ l .τ remains constant at a value representative of a moister soil, thus favoring vertical mixing of soil water and leading to a nearly uniform enrichment with depth [167][168][169][170].
To summarize, our representation of isotopic vertical profiles in ORCHIDEE is probably most suited when soil moisture remains high and does not vary too strongly.

Calculation of isotopic forcing from LMDZ outputs and nearby GNIP or USNIP stations
When precipitation and water vapor isotopic observations are not available at a given site, we create isotopic forcing using isotopic measurements in the precipitation performed on nearby GNIP (Global Network for Isotopes in Precipitation [122]) or USNIP (United States Network for Isotopes in Precipitation [125]) precipitation stations. To interpolate between the nearby stations, taking into account spatial gradients and altitude effects, we use outputs from an LMDZ simulation.
Let's assume there are n GNIP or USNIP stations around the and where d i is the geographical distance between the site of interest and the GNIP or USNIP station, δ p,lmdz (s) is the precipitation isotopic composition simulated by LMDZ in the grid box containing the site s, δ p,lmdz (i) is the precipitation isotopic composition simulated by LMDZ in the grid box containing the GNIP or USNIP station, δ p,NIP (i) is the precipitation isotopic composition observed at the GNIP or USNIP station, Z site is the altitude of the site of interest, Z lmdz (s) is the altitude of the LMDZ grid box containing the site of interest and a s is the slope of the isotopic composition as a function of altitude simulated by LMDZ in the grid boxes containing and surrounding the site of interest [171]. The first term on the right hand side corresponds to the raw LMDZ output for the site of interest. The second term allows us to correct for the altitude effect. Since LMDZ is run at a 2.5 ° latitude * 3.75 ° longitude resolution, we cannot expect the average grid box size to be representative of the local altitude at the site. The third term allows us to correct for possible biases in LMDZ compared to GNIP and USNIP observations. Table 3 lists the GNIP and USNIP stations used to construct the forcing at each site of interest.
To calculate the isotopic composition of the water vapor, we assume that although LMDZ might have biases for simulating the absolute values of precipitation and water vapor composition, it simulates properly the precipitation-vapor difference [47,60]. Therefore, the isotopic composition of water vapor at the site of interest, δ v,site, is calculated as: δ v,site =δ p,site +δ v,lmdz (s)-δ p,lmdz where δ v,lmdz (s) is the isotopic composition of water vapor simulated by LMDZ in the grid box containing the site of interest.

A simple equation to relate the soil water isotopic composition to the surface soil water budget
To explore how the isotopic composition of soil water can help estimate terms of the soil water budget, we derive here a very simple theoretical framework.
We assume that the water mass balance is: where P is the precipitation, R the surface runoff, E is the bare soil evaporation, T the transpiration and D the drainage. Similarly, the isotopic mass balance is: P.R p =E.R E +T.R T +D.R D +R.R R (11) Where R p, R E , R T, R D and R R are the isotopic ratios of incoming water at the soil surface, bare soil evaporation, transpiration, drainage and surface runoff respectively.  [152]. a) The soil column evaporates at its top and is permanently refilled at the bottom by a water with δ 18 O =-8 ‰ b) The soil column is evaporated progressively until its soil water content is only 20%. See appendix 8.2 for more details. Simulations using the soil profile module of the isotopic version of ORCHIDEE (colors) with different parameters and vertical resolution are compared with the more sophisticated MuSICA and SiSPAT models and with an analytical solution. For θ l .τ =0.005, the vertical resolution for ORCHIDEE is 0.15 mm for the first layer and 0.75 mm (resol=5), 1.5 mm (resol=10), 3 mm (resol=20) or 6 mm (resol=40) for the other layers. For θ l .τ =0.01, the vertical resolution for ORCHIDEE is 0.21 mm for the first layer and 2.12 mm (resol=10) or 4.24 mm (resol=20) for the other layers. For θ l .τ =0.1, the vertical resolution for ORCHIDEE is 0.67 mm for the first layer and from 1.34 mm (resol=2) to 3.35 mm (resol=5) for the other layers. We assume that the bare soil evaporation isotope ratio depends on that of the soil (R s ) following the Craig [75] relationship (Equation 2) and that the transpiration composition is equal to that of the soil (R T =R s ), implying little vertical variations in soil water isotope ratios [172]. We assume that the isotopic composition of surface runoff is that of the incoming water (R R =R p ) and that the isotopic composition of drainage is that of the soil water (R D =R s ). In doing so, we neglect again vertical isotope variations in the soil and the temporal co-variation between R s, D and T. Combining equations for the mass balance of water (Equation 11) and of water isotopes (Equation 10) then yields: where I=P-R represents the incoming water that infiltrates into the soil. E/I represents the proportion of the infiltrated water which is evaporated at the soil surface.
The composition of the bare soil evaporation flux, R E , is a function of R s following the Craig [75] Therefore, E/I is a function of the isotopic difference between the soil water and the precipitation water, which is easy to observe on instrumented sites such as MIBA or Carbo-Europe sites.