Stable water isotopes in the coupled atmosphere – land surface model ECHAM 5-JSBACH

In this study we present first results of a new model development, ECHAM5-JSBACH-wiso, where we have incorporated the stable water isotopes H 2O and HDO as tracers in the hydrological cycle of the coupled atmosphere–land surface model ECHAM5-JSBACH. The ECHAM5-JSBACH-wiso model was run under present-day climate conditions at two different resolutions (T31L19, T63L31). A comparison between ECHAM5-JSBACH-wiso and ECHAM5-wiso shows that the coupling has a strong impact on the simulated temperature and soil wetness. Caused by these changes of temperature and the hydrological cycle, theδ18O in precipitation also shows variations from −4 ‰ up to 4 ‰. One of the strongest anomalies is shown over northeast Asia where, due to an increase of temperature, the δ18O in precipitation increases as well. In order to analyze the sensitivity of the fractionation processes over land, we compare a set of simulations with various implementations of these processes over the land surface. The simulations allow us to distinguish between no fractionation, fractionation included in the evaporation flux (from bare soil) and also fractionation included in both evaporation and transpiration (from water transport through plants) fluxes. While the isotopic composition of the soil water may change for δ18O by up to +8 ‰, the simulatedδ18O in precipitation shows only slight differences on the order of ±1 ‰. The simulated isotopic composition of precipitation fits well with the available observations from the GNIP (Global Network of Isotopes in Precipitation) database.


Introduction
Since Dansgaard (1964) explored the coherence between the isotopic composition of H 2 16 O, H 2 18 O, and HDO in precipitation and climate variations, stable water isotopes have proven to be a useful tool for understanding climate variations and climate changes in the past.The composition of stable water isotopes as recorded in various paleoclimate archives (e.g., in ice cores, sediment cores, corals, tree rings, or speleothems) have been used to reconstruct temperature and other climate changes of the past.This is possible as the stable water isotopes differ by their mass and symmetry of their molecules.As a result, they behave differently at any phase transition of a water mass within the hydrological cycle on Earth.While the heavier molecules H 2 18 O and HDO tend to stay in the liquid or solid phase, the lighter H 2 16 O molecules evaporate more easily.The strength of this partitioning effect, called fractionation, depends on the surrounding environmental conditions, with temperature as one of its key influencing parameters.
However, the interpretation of the isotope proxy data (usually expressed in a δ-notation) is often not straightforward because the proxy data include a mixture of fractionation processes occurring during evaporation (from bare soil or open water bodies) and transpiration (through plants) of liquid water, mixing of water masses of different origin, and fractionation during condensation processes, leading to the final isotopic composition of precipitation.Furthermore, the measured isotopic signal may also be affected by local postdepositional surface processes, e.g., for terrestrial archives by river runoff or percolation through the soil, or for ice cores by wind erosion or sublimation.
An enormous benefit of modeling stable water isotopes is the ability to directly compare field data to modeled isotope data.Thus, the models can be evaluated with present-day observational data found for example in the GNIP (Global Network of Isotopes in Precipitation) database (IAEA/ WMO, 2006).Furthermore, the interpretation of the measured variations of isotopes can be supported by model simulations.Studies like those of Jouzel et al. (2000), Vuille and Werner (2005), Herold and Lohmann (2009), and Risi et al. (2010) show that the interpretation of proxy data benefits from the addition of isotope modeling.
Over land surfaces two main processes exist which include a phase transition of water masses: evaporation and transpiration.Whereas isotope fractionation occurs during an evaporation process, it is often assumed that the transpiration is a non-fractionating process (see Gat, 1996).Many of the presently existing GCMs enhanced with isotopes do not consider such difference between the evaporation and transpiration flux but simply assume that the whole evapotranspiration from land surface is a non-fractionating process (see, e.g., Hoffmann et al., 1998, for a more detailed discussion of this issue).So far, only very few GCM studies, e.g., Aleinov and Schmidt (2006), have started to investigate fractionation processes over land.
In this study, we present the first results of a newly developed isotope scheme within the ECHAM5-JSBACH model (named ECHAM5-JSBACH-wiso hereafter).The model is built from two separate components, the atmosphere model ECHAM5 (Roeckner et al., 2003) and the land surface scheme JSBACH (Jena Scheme for Biosphere-Atmosphere Interaction in Hamburg; Raddatz et al., 2007).The atmosphere isotope processes in this coupled model are almost identically implemented as in the stand-alone ECHAM5wiso model version (Werner et al., 2011), while the isotopic diagnostics within land surface processes are a novel development for JSBACH.With this setup it is possible to distinguish between the two partial fluxes of evapotranspiration, evaporation and transpiration, and separately incorporate the relevant fractionation processes for both fluxes.
We focus in our study on two questions: first, what are the implications of using ECHAM5-JSBACH-wiso instead of ECHAM5-wiso?Here we examine key variables of JSBACH, which can influence the atmospheric water cycle in ECHAM5, and the related changes of the isotopic composition of precipitation.Next, we analyze the sensitivity of the isotope results of ECHAM5-JSBACH-wiso to different assumptions regarding the fractionation processes over land.In general, any isotopic fractionation during evaporation consists of two parts: an equilibrium fractionation occurring between the liquid water and a thin, saturated vapor layer above the water mass, plus a kinetic fractionation process occurring during the diffusion of the water molecules from the saturated vapor layer into the undersaturated free atmosphere (Gat, 1996).For the equilibrium fractionation we perform sensitivity studies to distinguish between three different approaches.First, we assume that no fractionation during evapotranspiration occurs at all, similar to the approach used in the ECHAM5-wiso model (Werner et al., 2011).Second, we assume that fractionation only occurs during evaporation from bare soil but not during transpiration.Last, we consider that fractionation processes take part during both evaporation and transpiration of water from land surface.For the impact of the kinetic fractionation factor, we additionally analyze two different formulations given by Merlivat and Jouzel (1979) as well as Mathieu and Bariac (1996).
In the following section we give a detailed description of the ECHAM5-JSBACH-wiso model.Furthermore we explain the performed set of simulations as well as the selection of observational data for evaluating the model results.The comparison of ECHAM5-JSBACH-wiso and ECHAM5-wiso follows in Sect.3.1.In Sect.3.2 we investigate the sensitivity of the impact of fractionation over land, and distinguish between the equilibrium fractionation and the relevance of the kinetic fractionation factor.The final section of this manuscript includes the conclusion and an outlook.

Model description
ECHAM5 is an AGCM, developed mainly at the Max Planck Institute for Meteorology, Hamburg, that consists of a spectral, dynamical core based on the equations of conservation of momentum, mass and energy.This set of equations is completed by the hydrostatic equation, the continuity equation, and a prediction equation for the surface pressure (Roeckner et al., 2003).The hydrological cycle in the model consists of the formulations for evaporation of ocean water, evapotranspiration of terrestrial water, two schemes for the formation of large-scale and convective clouds, as well as an independent advective transport of vapor, liquid and frozen water within the atmosphere.A detailed description of the physics of the model as well as changes to the earlier model version can be found in Roeckner et al. (2003).
For the coupled ECHAM5-JSBACH model, the JSBACH routines calculate the terrestrial boundary conditions for the ECHAM5 over the land surface for each time step.This includes a simulation of the exchange of energy, water, and momentum between the land surface and the atmosphere.JSBACH is based on the ECHAM3 surface hydrology (DKRZ, 1992), which is also used by ECHAM5, and the biosphere model "Biosphere Energy Transfer and Hydrology scheme", called BETHY (Knorr, 2000).The basic idea of the model structure is a partitioning of the land surface.Each grid cell includes 8 tiles, which represent the fraction covered by one of the plant functional types (PFTs), distinguishing between tropical and nontropical as well as deciduous and evergreen trees, deciduous and evergreen shrubs, C3 grasses, and C4 grasses, as well as seasonally bare soil and permanently bare soil, i.e., desert (Raddatz et al., 2007).The simulated vegetation is based on temporal change of growing, natural mortality, and disturbance mortality (e.g., wind, fire).The modeling of vegetation and its dynamics are explained in detail by Brovkin et al. (2009).
In ECHAM5-JSBACH the same land hydrology model is used as in ECHAM5.The model comprises three surface water reservoirs: a snow layer (sn), water at the skin layer of the canopy or bare soil (wl), and a soil water layer (ws).These three types are each represented by a single layer bucket model, and each of them has a prescribed maximum field capacity.The snow reservoir is filled by snowfall and depleted by snowmelt or sublimation.The skin layer and the soil layer are filled by rainfall and snowmelt in the following order: first the skin layer is filled until its water holding capacity is exceeded, and secondly the non-intercepted water fills the soil reservoir.The modeled depletion of the skin layer can only occur by evaporation; the depletion of the soil water reservoir occurs by evapotranspiration.There is no exchange between these two reservoirs.If the soil water reservoir is saturated, surface runoff occurs.Drainage occurs independent of the new precipitation, and it is calculated if the amount of soil water reaches 5 % or more of the maximal soil water capacity.The runoff and drainage scheme is based on examination of variations of the field capacity for soil water over the land surface (Dümenil and Tondini, 1992).Furthermore, lakes are prescribed by a lake mask; to calculate the evaporation over larger lakes (i.e., grid cells with a lake fraction greater than 50 %) the same scheme as for the ocean is used.A more detailed description of the land hydrology model can be found in Roeckner et al. (2003).
As in the stand-alone atmosphere model ECHAM5-wiso the water isotope tracers in ECHAM5-JSBACH-wiso are implemented parallel to the normal water cycle.Fractionation of H 2 18 O and HDO versus H 2 16 O occurs during every phase change.Aside from fractionation during evapotranspiration from the land surface, all fractionation processes in ECHAM5-JSBACH-wiso are implemented in an identical manner to ECHAM5-wiso.For evaporation over the ocean, we use the bulk formula described by Hoffmann et al. (1998).This equation includes the dependence of the isotope evaporation flux on the isotopic compositions of water vapor close to the ocean surface, evaporation temperature, relative humidity, and wind speed at the ocean surface (Hoffmann et al., 1998).The implementation of fractionation processes inside the cloud schemes, specifically during cloud formation, are described in detail by Werner et al. (2011).Furthermore, as in ECHAM5-wiso we use the assumption that convective showers generate primarily large raindrops equilibrating isotopically to only 45 % as they fall through an undersaturated atmosphere, and that large-scale clouds generate smaller rain drops equilibrating nearly completely (95 %) with their surrounding (see Hoffmann et al., 1998, for details).
The water isotope tracers are almost passive in the land surface scheme JSBACH.So for example, during surface runoff and drainage the stable water isotopes are a completely passive tracer and are following the normal water.The runoff is calculated as a composition of precipitation and snowmelt.The same is valid for the calculation of its isotope ratio.The drainage has the isotopic composition of the soil water.We also assume no fractionation during snowmelt.Thus, the melt water has the same isotopic composition as the snow.The melt water is added to the skin reservoir and the soil reservoir, respectively.After these reservoirs are filled the residual melt water is added to the runoff.
The only exception is the evapotranspiration.In order to calculate the evapotranspiration in ECHAM5-JSBACH, each grid cell is divided into four cover fractions: the fraction C sn covered by snow, the fraction (1 − C sn )C wl covered with water in the skin reservoir, the fraction (1 − C sn )(1 − C wl )C veg covered by vegetation, and the fraction covered by bare soil.The complete evapotranspiration flux is calculated by the weighted sum of these four fractions.The evapotranspiration from the surface to the atmosphere is implemented with a negative sign convention.In order to incorporate the stable water isotopes in JSBACH we follow the same method.
Water sublimates from snow at a potential evaporation rate, which is given by the following equation: with q sat as the saturation-specific humidity at the corresponding surface temperature, q vap as the humidity of the air level directly above surface, v h as the horizontal wind speed at the surface, C V as the drag coefficient for water flux, and ρ as the density of air.Since the diffusion rate in the ice crystal structure is very low, we assume no fractionation occurs during sublimation, which leads to the model assumption that the evaporative flux from snow has the same isotopic composition as the snow itself (Here and in the following paragraph we use I x for the amount of an isotopic species and R x for the ratio of a isotope species with respect to H 2 O, with x ∈ H 2 16 O, H 2 18 O, HDO .)This assumption leads to the following equation for the isotope flux during snow sublimation: Analogously to Eq. ( 1) evaporation from the skin layer (wl) in ECHAM5-JSBACH is calculated as (3) The skin layer wl is modeled as a thin layer of water, which in general evaporates completely within a few model time steps.If this entire water reservoir evaporates, the total flux has an identical isotopic composition as the water source and no fractionation occurred.As this study focuses on annual mean changes, we assume for simplification that no fractionation occurs during evaporation from skin layer at any time step, which is expressed as In ECHAM5-JSBACH, the following equation is used for the evaporation from bare soil (bs): with h ws as the relative humidity of the soil surface (DKRZ, 1992).
To calculate the fractionation during evaporation over land surface the same bulk formula is used as described by Hoffmann et al. (1998).So, to calculate the fractionation during evaporation, we use the equilibrium fractionation factor α x (T ), with T as surface temperature, obtained from Majoube (1971), which results in the temperature dependency of the isotopic composition of evaporation, and a factor for kinetic fractionation (α k ).Furthermore the mixing ratio of the water isotopes in the layer above the surface q x vap and the isotopic ratio of the saturation mixing ratio q x qsat analogue to Eq. ( 5) are needed.While q x vap = R x vap q vap is known in the atmosphere component of the model, q x qsat can be calculated with q x qsat = R x qsat q sat .Here R x qsat is the isotopic ratio of the saturation-specific humidity.If we use the equilibrium fractionation factor, R x qsat can be expressed by using the isotopic ratio of soil water with R x qsat = R x ws α x (T ) .So, the evaporation from land surface enhanced with fractionation is described by The term α k in Eq. ( 6) includes the non-equilibrium fractionation effects, taking into account the kinetics during the diffusion of vapor from a thin layer just above the soil water into the free atmosphere.For the calculation of the kinetic fractionation two different approaches are tested.First, we assume that the same kinetic fractionation factor as for evaporation over the ocean can be used over land as well: . Here V s is the horizontal wind speed on the surface and λ describe the ratio of the isotope molecular diffusivity in air.In this approach α k is dependent on the molecular and turbulent resistance of water vapor and has been described in detail by Merlivat and Jouzel (1979).For this approach, typical values of α k for H 18 2 O range between 0.994 and 0.998, and according to Eq. ( 7) the values for HDO are slightly larger.
The second approach is presented by the study Mathieu and Bariac (1996), where α k is calculated as the nth power of the molecular diffusivity ratio in air: with d v (d x v ) as the vapor diffusivity in air (vapor diffusivity of the isotopic species x).The exponent n includes the influence of the turbulent and molecular resistance, and we use, as suggested by Riley et al. (2002), n = 0.67, which results in α k = 0.981 for H 2 18 O and α k = 0.984 for HDO.The impact of these two different kinetic fractionation factors on the isotopic composition of the different modeled water reservoirs is analyzed and discussed in detail in Sect.3.2.3.
Additionally to Eq. ( 6) we implement a second approach for evaporation from bare soil, based on the assumption that no fractionation during evaporation over land surfaces occurs.This leads to the modified formulation of Eq. ( 6) with Ẽx bs = R x ws ρ C V |v h | q vap − h ws q sat .This setup is identical to the implementation of ECHAM5-wiso and allows a comparison between this two models.
Since the hydrology inside the plants is not described by ECHAM5-JSBACH, the transpired water is modeled as a potential transpiration flux: The factor S −1 is the transpiration efficiency, which includes among others the stomatal resistance of canopy.A detailed description can be found in DKRZ (1992).Gat (1996) has shown that there is no fractionation between isotopes as roots take up water.This leads to the model assumption that the isotopic composition (R x veg ) inside the plants is identical to the isotopic composition of the soil water (R x ws = R x veg ).If we assume no fractionation occurring during transpiration, the transpiration isotope flux is calculated as follows: To estimate the potential maximum fractionation effect for the combined evapotranspiration flux over land surface, we perform an additional sensitivity study.Here we assume that the equilibrium fractionation occurs during both evaporation and transpiration.As the JSBACH model does not resolve the hydrology inside the plants and does not simulate the amount of leaf water, we assume that the whole amount of transpired water can fractionate.This leads to the altered Eq. (10): ws α x (T ) q sat .We are aware that this sensitivity study does not mimic the natural process of isotope changes during transpiration (e.g., described by Sachse et al., 2012).Nevertheless we rate it as a useful for estimating the upper limit of isotope changes related to the simulated evapotranspiration in ECHAM5-JSBACH.
Dew formation occurs in ECHAM5-JSBACH if the vapor of the lowest model layer q vap is larger then the saturationspecific humidity q sat .For this case, we assume the equilibrium fractionation between the dew and the surrounding vapor.

Simulation setup
All simulations are run under present-day conditions with a prescribed vegetation distribution over a simulation period of 10 yr after a spin-up period of 2 yr.We distinguish between the model resolutions T31L19 (horizontal grid size 3.8 • × 3.8 • , 19 vertical model levels) and T63L31 (1.8 • × 1.8 • , 31 levels).The simulations are performed with AMIP-conform present-day boundary conditions including prescribed climatological sea surface temperatures and sea ice cover for the period 1979-1999 (see Taylor et al., 2000).The lower oceanic boundary condition for the atmospheric 18 O isotopic composition is based on the dataset described by LeGrande and Schmidt (2006).This is a global gridded dataset for sea surface water and sea ice.As no equivalent dataset is available for the composition of HDO we use as the lower oceanic boundary condition for the isotopic composition of deuterium the observed relation for meteoric water on a global scale (Craig and Gordon, 1965) and assume δD = 8 • δ 18 O for sea surface water and sea ice.
To evaluate the sensitivity of the fractionation processes over land we use a set of present-day simulations with various fractionation schemes implemented.The fractionation process over land will be varied between no fractionation (simulation named noF), fractionation occurring during evaporation only (FE), and the idealized setup with fractionation occurring during both evaporation and transpiration (FET).These three cases are all performed without any additional kinetic fractionation (α k = 1).
In order to investigate the influence of the kinetic fractionation of terrestrial evaporation on the isotopic composition of the different water reservoirs we use the FE fractionation scheme extended by two different calculations of the kinetic fractionation factor α k .The first setup, called FEK openwater , uses the same kinetic fractionation factor over land surface as over the ocean (Eq.7).The second setup calculates α k in dependence on the diffusion resistance (Eq.8) and is called FEK diffres .
For a comparison of the ECHAM5-JSBACH-wiso results with the stand-alone ECHAM5-wiso model, we use two comparable present-day ECHAM5-wiso control simulations in T31L19 and T63L31 resolution, from Werner et al. (2011).

Observational data
As observational data for evaluating the model results we choose the GNIP database.Since 1961, the International Atomic Energy Agency (IAEA) and the World Meteorology Organization (WMO) have collected monthly precipitation samples at more than 800 meteorological stations in 101 countries.Additional information and the available data can be found in IAEA/ WMO (2006).For this study we choose 248 GNIP stations where isotope data have been recorded for at least three consecutive years within the time period 1961 to 2008, and where at least 10 months of data per year are available.As a further restriction, we only use stations which provide a full monthly mean dataset, including values of 2 m air temperature (T 2 m ), precipitation amount (P ), and the isotopic composition of precipitation (δ 18 O P and δD P ).We are aware that a period of three years is perhaps too short to represent a long-term to represent a longterm climatological mean value at the stations' locations.On the other hand there are only 74 GNIP stations which have collected 10 yr or more of data.Since most of them are located in central Europe, many regions in Asia, Americas, Africa, and Australia would be underrepresented in such a limited dataset.Therefore we opted for a three-year time period in order to be able using a globally more representative sample distribution.

Impact of the coupling of ECHAM5 and JSBACH
In order to get an impression of how the overall model results change by using ECHAM5-JSBACH-wiso instead of ECHAM5-wiso we first compare the simulated surface temperature, precipitation amount, and soil wetness results of both models.All these variables are independent of the isotope diagnostic scheme, and differences between simulation results of both models are related to the changed representation of land surface processes in ECHAM5-JSBACH as compared to the stand-alone ECHAM5 model.Then, we take a look at the simulated distribution of δ 18 O in precipitation (hereafter named δ 18 O P ).As no fractionation for evaporation and transpiration processes has been assumed in the ECHAM5-wiso model by Werner et al., (2011), we use the analogous ECHAM5-JSBACH-wiso setup (noF) for this comparison. .These temperature changes are strongly related to the variation of the simulated surface albedo (Fig. 2a), which shows an increase over Antarctica as well as Greenland and a decrease over North America and Eurasia.For the finer-resolution T63L31 (not shown) most of the anomaly pattern are similar.Only in some regions, due to the finer description of regional attributes, deviations between the anomalies are detected.Simulated surface temperatures differ by 1 • C, or less.

Surface temperature, precipitation amount, and soil wetness
The simulated soil wetness differs between both models as well (Fig. 1d).The most notable changes are in the Amazon region, where an increase of 0.2 m is present, and in southern Africa, where a decrease of 0.25 m can be seen.There is also a clear increase in a range of 0.08 m to 0.15 m over the Sahara.Most locations displaying a decrease in soil moisture generally show also an increase of evapotranspiration, which can be linked to changes in the simulated surface temperature.The increase of soil moisture in the Amazon region and Saharan Africa can be directly linked to an increase of the prescribed maximum soil water capacity (Fig. 2b).This difference between ECHAM5 and ECHAM5-JSBACH was introduced to enable a more realistic simulation of vegetation coverage over the tropical regions (Hagemann et al., 1999).It was only introduced for T31L19 resolution.Consequently, no similar soil water anomalies are found in the corresponding T63L31 simulation.Furthermore, in T63L31 a slight increase of soil wetness is simulated over Australia.This could be related to the finer resolution of albedo which results in a temperature change.
The simulated mean annual precipitation amount (not shown) shows nearly the same pattern in both models.For resolution T63L19 less annual mean precipitation in the range of 30-60 mm month −1 (which corresponds to 0.5-4 % of the annual mean precipitation amount at the various locations) is simulated over central and southern Africa and over India.

Isotopic composition of precipitation and soil water
Figure 3 shows the simulated δ 18 O in precipitation (δ 18 O P ) using the noF setup (no fractionation during evaporation and transpiration) of ECHAM5-JSBACH-wiso for both model resolutions, T31L19 (Fig. 3a) and T63L31 (Fig. 3b).Both simulations show the typical δ 18 O P pattern described by Dansgaard (1964).We see a depletion from the tropics to the high latitudes (temperature effect) as well as a depletion from the oceans to the landmasses of North America and Eurasia (continental effect).A depletion of δ 18 O P above the mountain areas can also be identified (altitude effect), for example for the Andes.However, Fig. 3b also shows that the altitude effect is better represented in the higher model resolution T61L31.The root mean square error (RMSE) between the simulations and the GNIP data is 2.15 ‰ for T31L19 and 1.78 ‰ for T63L31, which shows that the simulated δ 18 O P values generally improve for a higher ECHAM5-JSBACHwiso model resolution.For the analogue simulations with the ECHAM5-wiso model the calculated RMSE with respect to the same set of GNIP stations is 2.25 ‰ for T31L19 and 1.89 ‰ for T63L31.Thus, both models show similar results for δ 18 O P on a global scale.
In order to further analyze the impact of the coupling of ECHAM5 with JSBACH for the simulation of stable water isotopes, we calculate the difference of δ 18 O P between ECHAM5-JSBACH-wiso and the ECHAM5-wiso simulations for both resolutions (Fig. 4).Due to the relatively short simulation period of 10 yr, we exclude in our analyses δ 18 O P changes in the range of −1 ‰ to +1 ‰, as such small differences might be caused by internal model variability only.For T31L19, the strongest differences with an increase up to  approx.+4 ‰ are located in the region of north Africa to the Arabian Peninsula.These anomaly can be related to a decrease in the amount of precipitation in ECHAM5-JSBACHwiso related to ECHAM5-wiso.Negative anomalies which are below −1 ‰ are only simulated in the high latitudes over Greenland, Antarctica and northeast Russia.Over Antarctica and Greenland the changes are most likely due to the different temperatures simulated in this region (see Fig. 1).Over northeast Russia the anomalies can be linked to an increase of precipitation.Largest differences between the resolutions are found for northeast Russia, wich is most likely linked to warmer temperatures and reduced regional precipitation in this simulation.
For a further model evaluation we investigate the relationship between δ 18 O P and temperature 2 m above the surface (δ 18 O P -T 2 m ) as well as δ 18 O P and the amount of precipitation (δ 18 O P -P ).For the δ 18 O P -temperature relationship we use those 186 GNIP stations where the annual mean temperature is below 20 • C. Figure 5 shows the simulated δ 18 O P -T 2 m relation for both ECHAM5-JSBACHwiso and ECHAM5-wiso.Both models show a similar δ 18 O P -T 2 m relation as derived from the GNIP data, but slightly overestimate δ 18 O P .The simulated strong correlation between δ 18 O P and T 2 m in ECHAM5-JSBACH-wiso is statistically significant for both model resolutions (Pearson Correlation coefficient: R 2 = 0.816 for T31 and R 2 = 0.845 for T63), similar to the observed correlation at the GNIP stations (R 2 = 0.909).As seen in Fig. 5, the simulated δ 18 O P -T 2 m relation also slightly improves for the finer model resolution T63L31.For the correlation of δ 18 O P and precipitation we choose the other 62 GNIP stations with a mean annual temperature above or equal to 20 • C. The simulated relation fits quite well to the observed relation for both model resolutions (Fig. 6) with a slight tendency to underestimate the δ 18 O P -P relation in the T31L19 resolution (both ECHAM5-wiso and ECHAM5-JSBACH-wiso).We refrain from a more quantitative analysis of the simulated δ 18 O P -T 2 m and δ 18 O P -P relation in this study as both the simulated and observed mean δ 18 O P , T , and P values may contain relatively large uncertainties due to the short simulation (10 yr) and GNIP observation (3 yr or more) period.
In summary, the analyses show that the coupling of the atmosphere model ECHAM5 with the surface scheme JSBACH has a detectable impact on the simulated temperature, evapotranspiration, and soil wetness.These changes are related to the alteration in the simulated surface albedo parameters and the prescribed maximum soil wetness.The simulated precipitation amount is less strongly influenced by the coupling.Since the isotopic composition of precipitation highly depends on these variables, the coupling of ECHAM5 with JSBACH also has a noticeable impact on the simulated δ 18 O P values in various regions.However, our analyses also reveal that the ECHAM5-JSBACH-wiso model is capable of simulating a global distribution of δ 18 O P in a good overall agreement with available observations from GNIP stations, similar to previous results retrieved with the stand-alone ECHAM5-wiso model.

Fractionation processes over land surfaces
In this section we investigate the sensitivity of the ECHAM5-JSBACH-wiso simulation results regarding different assumptions for both the equilibrium fractionation (Sect.3.2.1)as well as the kinetic fractionation (Sect.3.2.3)over land surface.All simulations in this part of our study are performed at resolution T31L19.

Equilibrium fractionation during evaporation and transpiration
When water evapotranspirates from the land surface, it can evaporate from bare soil or skin layer, sublimate from snow, or transpire through the vegetation.According to Wang and Dickson (2012), transpiration is the largest contribution to evapotranspiration on a global scale.This relevance of transpiration is also seen in the ECHAM5-JSBACH-wiso simulations.In Fig. 7 the modeled annual mean evapotranspiration flux from land surface (Fig. 7a) and the fraction of evaporation in relation to the total evapotranspiration flux (Fig. 7b) are shown.Especially in the (sub)tropical regions, transpiration is the dominant water flux from land surface to the atmosphere, while evaporation dominates over transpiration mainly in northern high-latitude regions as well as the Tibetan Plateau.
For the incorporation of stable water isotopes in GCMs or land surface schemes, various assumptions for the description of the equilibrium fractionation process during evapotranspiration have been utilized.In order to investigate the influence of fractionation processes over land surface on the isotopic composition of precipitation we assume two extreme cases for transpiration in our sensitivity studies: for one model setup (named FE) we assume isotope fractionation during evaporation processes only, and for another setup (FET) we assume isotope fractionation during the complete simulated evapotranspiration flux.As a third case comparable to many previous GCM studies we examine the case (noF) of no fractionation occurring during evaporation and transpiration at all.
Figure 8 shows the anomalies of the modeled annual mean δ 18 O ws between the FE-noF (Fig. 8a) and the FET-noF setup (Fig. 8c), as well as the modeled anomalies of δ 18 O P between the FE-noF (Fig. 8b) and the FET-noF setup (Fig. 8d).For the comparison of FE and noF (Fig. 8a) we see a relative stronger enrichment of δ 18 O ws in the FE setup from 0.5 ‰ to 4 ‰ in regions with a relatively high evaporation rate (Fig. 7).The fractionation during evaporation leads to a relative depletion of near-surface vapor in the FE setup as compared to the noF setup.This change in vapor leads to a slight depletion of δ 18 O P in the FE setup, compared to the noF one, ranging from −0.7 ‰ to −0.1 ‰ over the regions of North America, most parts of Eurasia, the Amazon region, and southern Africa.However, over the Tibetan Plateau and northeast Africa, the δ 18 O P in the FE setup is relatively more enriched than in the noF setup, with differences in the range of 0.1 ‰ to 0.9 ‰.These enrichments are most likely a result of recycling of relatively enriched local soil water.
The anomaly plot of the isotopic composition of soil water of FET-noF (Fig. 8c) reveals a stronger enrichment of δ 18 O ws for the FET setup relative to the noF setup, as well as to the FE setup, over the tropics and mid-latitudes.The range of this enrichment is 0.2 ‰ to 10 ‰.Only at northeast Russia is a slight depletion of δ 18 O ws of approx.0.1 ‰ in FET setup compared to noF setup shown, which can be linked to the depletion of precipitation in this area.When using the FET setup instead of the noF one, a relatively stronger enrichment of modeled annual mean δ 18 O P in the range of 0.1 ‰ to 0.8 ‰ is detected over the region stretching from north Africa via the Arabian Peninsula to the Tibetan Plateau, southern Africa, Central America, the Amazon region and northern Australia.This enrichment is most likely caused by the recycling of the modeled enriched soil water due to the relatively high evapotranspiration rate at these areas.Furthermore, a stronger depletion of δ 18 O P from −1 ‰ to −0.2 ‰ is modeled over North America and northern Eurasia, where the strongest anomaly is shown over northeast Russia.
Next, we analyze how accurately the different setups FE, FET, and noF simulate δ 18 O (δD) values in precipitation as compared to the various present-day GNIP observations.Table 1 show the calculated correlation between simulated and observational values.For this calculation, we use the set of 248 GNIP stations described in Sect.2.3 and distinguish again between GNIP data of stations with a mean annual temperature T ≤ 20 • C (shown in Fig. 9a) and those stations with a mean annual temperature T > 20 • C (shown in Fig. 9b).For all three model setups, the calculated correlation between simulated and observational values is significant for δ 18 O P and δD P (see Table 1) and very similar for all setups.However, Fig. 9a also shows that all three simulations overestimate δ 18 O P for most of these GNIP stations.A slightly different result is found for GNIP stations with T > 20 • C (Fig. 9b).For these stations, δ 18 O P is in numerous cases underestimated in all setups.

Seasonal changes
In order to get a more detailed picture regarding the modeled isotope variations, we analyze the seasonal cycle of the simulations using the FE, FET and noF model setups.For this purpose we choose nine GNIP stations from different geographical positions where the seasonal cycle of vegetation, amount of precipitation, surface temperature and the influence of evaporation over land strongly vary and compare the ECHAM5-JSBACH-wiso results to these GNIP data.
The first two stations are located on islands, where the influence of evaporation from the land surface is negligible in comparison to evaporation from the surrounding ocean.The station Reykjavik is chosen to represent the high northern latitudes, and the GNIP station in Jakarta represents the tropics.Since the only distinguishing factor between the three model setups is the fractionation of evapotranspiration over land, one can assume that the model behaves similarly in all implementations for the selected islands.For Reykjavik (Fig. 10a) all simulations reveal a correct seasonal timing of temperature, precipitation, and δ 18 O P .While the simulated δ 18 O ws shows an enrichment of +2 ‰ in the FE and FET setups in comparison to the noF setup, the simulated δ 18 O P is very similar in the three setups.For Jakarta (Fig. 10b) the simulated evaporation and transpiration from land as well as the simulated soil wetness are undefined.For the surface temperature, there is a good agreement between the simulated and observed values in Jakarta, while the simulated precipitation is strongly overestimated in the period April till July.The δ 18 O P has a correct timing of the seasonal cycle, but slightly too-enriched values in fall.For all three model setups the simulated δ 18 O P is very similar.
Because some of the strongest depletion in δ 18 O P between the different ECHAM5-JSBACH-wiso sensitivity experiments takes place in North America and Eurasia (as seen in Fig. 8) we choose three stations of these regions for comparison: Vienna, Ottawa, and Yakutsk.At all these locations, strong seasonal variations of vegetation and temperature exist, but the amplitude of the temperature variations varies strongly.At Vienna (Fig. 11a), the simulated temperature fits well with the observations, but the simulated precipitation shows an overestimation during the winter and spring.For all three model setups, the δ 18 O P shows the correct seasonality but a slight offset in the range of +1 ‰ to +2 ‰ as compared to the GNIP values.Only in spring and summer do the three simulations differ in a range of ±1 ‰.The simulated temperature also fits well in Ottawa (Fig. 11b); (a) δ 18 O ws FE−noF 0246810 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 72 74 76 78 80 82 84 86 88 90 92 94 96 98 2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8 8.0 8.2 8.4 8.6 8.8 9.0 9.2 9.4 9.6 9  however all simulation setups the seasonality of precipitation.For δ 18 O P , the simulations results have a correct seasonal timing, but all simulations overestimate the seasonal δ 18 O P amplitude, especially in summer.For Yakutsk, all simulations reveal a correct timing of the seasonality for temperature, precipitation, and δ 18 O P (Fig. 11c).While the seasonal amplitude of temperature and δ 18 O P agrees well with the GNIP observations, the ECHAM5-JSBACH-wiso model simulates too much summer precipitation in this region.For the noF, FE, and FET model setups, the simulated δ 18 O P is very similar except for the summer; there a difference up to 1 ‰ is detected between the simulations.
By comparison of the simulated soil wetness for the three GNIP stations Vienna, Ottawa, and Yakutsk, differences in the amplitude can be detected.The calculated amplitudes of the seasonal changes of the soil water bucket depth are approximately 11 cm for Vienna, 8 cm for Ottawa, and 3 cm for Yakutsk.Furthermore, the time interval in which transpiration takes place varies for these three stations: the longest, from March to November, is simulated for Vienna; a similar range (April to November) is simulated for Ottawa; and for Yakutsk only an interval from June to October is calculated.
To analyze the model performance in arid areas or areas with strong seasonal precipitation changes, we examine the stations Alexandria (Fig. 12a), Bamako (Fig. 12b), Kinshasa (Fig. 12c), and Addis Ababa (Fig. 12d).Alexandria is located in a very arid area with a dry season between May and September.This dry season is well simulated in ECHAM5-JSBACH-wiso, but the winter precipitation in the model is underestimated.Both temperature and δ 18 O P agree well with the GNIP observations with a slight overestimation of the simulated δ 18 O P .Furthermore, the ECHAM5-JSBACH-wiso simulates a very thin soil bucket depth (approx.0.05 m) as well as a very small evapotranspiration flux.While for FE and noF the simulated δ 18 O ws is nearly the same, the δ 18 O ws in FET setup is approx.0.5 ‰ heavier.For Alexandria all three simulations show the same weak seasonality for the isotopic composition of soil water.For Bamako (Fig. 12b observations.The simulated δ 18 O P values are approximately the same for all implementations, with too-depleted δ 18 O P values in the dry season between January and May as compared to the GNIP data.The peak of the summer depletion is simulated with a delay of one month.For the soil water bucket depth, the ECHAM5-JSBACH-wiso model simulates a strong seasonality (from 0.4 m during the dry season to 0.85 m during the wet season), but the related δ 18 O ws values of the noF and FE setup display weak seasonal variations.Additionally, these two simulations have more or less the same δ 18 O ws values.Only for the FET setup are strong seasonal changes of δ 18 O ws simulated.Similarly to the situation at Bamako, the monthly temperature and precipitation model results for Kinshasa (Fig. 12c) fit well to the observations.One major exception is an underestimation of the modeled precipitation amount in November.The simulation results reveal also a strong seasonality of the soil water bucket depth (from 0.20 m during the dry season to 0.40 m during the wet season).Again, the simulated δ 18 O ws values for FE and noF are more or less the same with a weak seasonal cycle, while the FET results show a strong seasonal cycle, inversely correlated to the seasonality of ws.Furthermore, the modeled δ 18 O ws values for the FET setup are stronger enriched by 3-8 ‰ when compared to the noF or FE setup.These differences of the noF or FE, and FE setup, in combination with the amount of evaporation, are directly imprinted in the simulated δ 18 O P values at the location Kinshasa.At Addis Ababa (Fig. 12d), the simulated temperature is strongly overestimated by +5 • to +12 • C. Modeled precipitation values have a correct seasonal timing, but the amount of summer precipitation is underestimated.The simulated soil wetness also shows a strong seasonality, which lags the seasonal cycle of precipitation by 2 months.The modeled δ 18 O ws values are almost constant in the noF and FE setup, but the FE setup is slightly more enriched.While the FET setup shows seasonal changes in δ 18 O ws inversely correlated to the seasonal cycle of soil wetness, the simulated seasonal cycle of δ 18 O P in all model setups is more or less the same, but does not agree with the GNIP observations.The performed sensitivity studies reveal that the various simulation results with the ECHAM5-JSBACH-wiso model are in relatively good agreement with the GNIP observations.In Fig. 9, it is shown that in many cases with a surface temperature T ≤ 20 • C the model rather overestimates the isotopic values in precipitation, while in cases with higher surface temperature (T > 20 • C) the simulated values are more often underestimated.The incorporation of fractionation effects during evaporation and transpiration in FE and FET setup does not lead to substantial improvements for δ 18 O P as compared to the noF setup and the observations (Fig. 9, Table 1).
Part of the model mismatch is probably related to the rather simple one-layer bucket model of soil water, implemented in the coupled ECHAM5-JSBACH model.When using a simple bucket model for the soil water, the whole soil water reservoir does have an identical isotopic composition.Any vertical moisture dynamics and changes of the isotopic composition with the soil moisture depth are neglected.But it is well known from observations (see, e.g., Allison and Hughes, 1983;Hsieh et al., 1998) that strong vertical isotope gradients in soil can exist.Enrichment does mainly occur in the upper soil layers, while water in deeper soil layers, which can be used for plant transpiration, is more depleted.Thus, a one-layer bucked model will most likely result in toodepleted isotope values of evaporated and too-enriched isotope values of transpired water.Furthermore, in a previous study, Schulz et al. (2000) analyzed the results of coupling the ECHAM model with various land surface schemes of different complexity.They showed that a bucket model tends to calculate higher evapotranspiration amounts than more complex schemes.Such overestimation will result in a too-strong influence of the isotopic composition of the soil water on the atmospheric isotopic composition and, consequently, on the isotopic values simulated in precipitation.

Sensitivity of kinetic fractionation
In order to examine the influence of the kinetic fractionation coefficient α k of terrestrial evaporation on the isotopic composition, we use the model setup FE (equilibrium fractionation occurring during evaporation only) extended by two calculations of the kinetic fractionation: for the first model setup (named FEK openwater ) we assume the same kinetic coefficient as over the ocean (see Eq. 7), which is presented in the study given by Merlivat and Jouzel (1979).The second setup (FEK diffres ) is based on the study given by Mathieu and Bariac (1996), where α k is calculated as the nth power of the molecular diffusivity ratio (see Eq. 8).As the third setup of the analyses we use FE, which has no kinetic fractionation included.Figure 13 shows the anomalies of the modeled annual mean δ 18 O ws between the FEK openwater -FE (Fig. 13a) and the FEK diffres -FE setup (Fig. 13c), as well as the modeled anomalies of δ 18 O P between the FEK openwater -FE (Fig. 13b) and the FEK diffres -FE setup (Fig. 13d).When using the FEK openwater or the FEK diffres setup instead of the FE one, a relative enrichment of δ 18 O ws of soil water in the range of 0.5 ‰ to 2.8 ‰ is detected in the areas where the evaporation is relatively high (Fig. 7).While FEK diffres leads to enrichment of δ 18 O ws only, the setup FEK openwater simulates positive as well as negative anomalies.Both setups, FEK openwater and FEK diffres , simulate the strongest impact of the kinetic fractionation in the northern high latitudes.This enrichment of soil water leads to a relative depletion of near-surface water vapor; as a result a stronger depletion of δ 18 O P is simulated for the setups including kinetic fractionation compared to the FE setup.The anomaly of FEK diffres -FE, with a depletion of −0.2 ‰ to −0.02 ‰, is stronger than the difference of FEK openwater and FE with a depletion in the range of −0.05 ‰ to −0.02 ‰.
Furthermore, we compare the simulated δ 18 O P values as well as the simulated relation of δD P and the deuterium excess (defined as dex P = δD P − 8 • δ 18 O P ) with the observational data.For these studies we use again those 246 GNIP stations described in Sect.2.3.Figure 14a depicts a comparison of the simulated annual mean δ 18 O P values with the observations.For all three model setups, the simulated δ 18 O P fits well with the observational values, but all three simulations overestimate the δ 18 O P for most of these GNIP stations.Moreover, Fig. 14a also shows that the calculated δ 18 O P is indistinguishable from the setups FE, FEK openwater , and FEK diffres .Figure 14b shows the simulated relation of δD P and the deuterium excess in precipitation (dex P ).It can be seen that the simulated δD P − dex P relation behaves very similarly for all three setups and shows a similar distribution in comparison to the GNIP data.(b) δ 18 O P FE openwater −FE 0.00 0.08 0.16 0.24 0.32 0.40 0.48 0.56 0.64 0.72 0.80 0.88 0.96 1.04 1.12 1.20 1.28 1.36 1.44 1.52 1.60 1.68 1.76 1.84 1.92 2.00 2.08 2.16 2.24 2.32 2.40 2.48 2.56 2.64 2.72 2.80 2.88 2.96 3.04 3.12 3.20 3.28 3.36 3.44 3.52 3.60 3.68 3.76 3.84 3.92 4.00 4.08 4.16 4.24 4.32 4.40 4.48 4.56 4.64 4.72 4.80 4.88 4.96 5.04 5.12 5.20 5.28 5.36 5.44 5.52 5.60 5.68 5.76 5.84 5.92 6.00 6.08 6.16 6.24 6.32 6.40 6.48 6.56 6.64 6.72 6.80 6.88 6.96 7.04 7.12 7.20 7.   The performed sensitivity test for the kinetic fractionation factor α k reveals that the setups FE, FEK diffres , and FEK openwater of the ECHAM5-JSBACH-wiso model simulate a different isotopic composition of the soil water.The simulations have shown that the setup FEK diffres leads to the strongest fractionation in terms of δ 18 O P as well as in terms of d-excess (not shown).However, the simulations of δ 18 O P as well as at the simulation of the δD P − dex P relation show no substantial difference between FE, FEK diffres , and FEK openwater .

Conclusions
In this study we show first simulation results of stable water isotopes successfully implemented in the coupled atmosphere-land surface model ECHAM5-JSBACH.The ECHAM5-JSBACH-wiso model is able to simulate the isotopic composition of precipitation (δ 18 O P and δD P ) in a good manner as compared to the stand-alone ECHAM5-wiso model.Furthermore we demonstrate that the relation between simulated temperature and δ 18 O P and between precipitation and δ 18 O P , respectively, is simulated in good agreement with the observations.An analysis of the impact of the coupling of ECHAM5 and JSBACH reveals that the simulated land surface temperature and surface albedo are remarkably influenced by the coupled setup and lead to some substantial regional changes of the hydrological cycle between the model ECHAM5-JSBACH and the stand-alone ECHAM5 model.This results in differences of the modeled soil wetness and evapotranspiration fluxes between the two models.
To investigate the importance of isotope fractionation processes over land surfaces, we use three different model setups.Our studies show that all three setups give relatively similar results.The simulations including fractionation over land result in a slightly higher depletion of δ 18 O in precipitation of up to −1 ‰ for both the FE and FET setup.For the FET setup, an enrichment of δ 18 O P on the same order of magnitude can occur for some (sub)tropical regions.As we assume an unrealistic fractionation of the total transpired water in our FET sensitivity studies, these enrichment effects are most likely much smaller (or even not existing at all) in reality.The inclusion of fractionation processes over land does not lead to substantial improvement of the simulated δ 18 O P in ECHAM5-JSBACH-wiso (Table 1) on a global scale.As in most places δ 18 O P is first and foremost controlled by atmospheric processes, we expect that even with a highly realistic land surface scheme the simulated δ 18 O P would not substantially improve at such locations.
In contrast to the minor simulated changes of δ 18 O in precipitation between the different model setups study) are simulated for the soil water reservoir in ECHAM5-JSBACH.At present, it is not possible to evaluate these simulated soil water changes by direct observations.A potential model-data comparison is hampered by the simple soil water scheme of ECHAM5-JSBACH.It is well known that the isotopic composition of soil moisture can strongly vary with depth.But since in ECHAM5-JSBACH a one-layer bucket model is used, it is not possible to simulate a vertical isotope profile within the soil.Thus, we cannot yet compare our simulation results with available isotopic vertical profiles within the soil.But such a comparison should be possible for a more complex multi-layer soil scheme, which might be implemented in a future JSBACH model release.Furthermore, recently started networks for isotopes in the biosphere, like MIBA (Moisture Isotopes in the Biosphere and Atmosphere) or BASIN (Biogeosphere-Atmosphere Stable Isotope Network), which monitor the isotopic compositions of soil water, have been operational for only a limited time so far.The available data do not yet represent long-term annual mean values.But such data will hopefully become available during the next few years and will then allow a much more profound model-data comparison of the isotopic composition of soil moisture on a global scale.
In the future, we are planning a set of Holocene simulations with the ECHAM5-JSBACH-wiso model which will distinguish between prescribed and dynamic vegetation.By using the ECHAM5-JSBACH-wiso model with dynamical vegetation we are able to investigate the feedback mechanisms between the hydrological cycle and the vegetation during the past.Moreover, the new isotope diagnostics will give the opportunity to compare the simulated isotopic composition of ECHAM5-JSBACH-wiso with available proxy data to improve our understanding of past hydrological changes.Furthermore, since the ocean model MPI-OM has also been enhanced with stable water isotopes (see Xu et al., 2012), we will be able to run simulations with a full coupled atmosphere-ocean-land-surface GCM including isotopes in the future.

FigureFig. 1 .Fig. 1 .
Figure1a and cshow the mean annual temperature and soil wetness as simulated by ECHAM5-JSBACH-wiso for the model resolution T31L19.The corresponding anomaly as compared to the comparable ECHAM5-wiso simulation

Fig. 5 .Fig. 5 .
Fig. 5. Comparison of the simulated δ 18 OP −T2m relation of ECHAM5-JSBACH-wiso (noF ) with ECHAM5-wiso observed for the resolutions (a) T31L19 and (b) T63L31.For comparison with the observed relation, we use data from those GNIP stations, where the annual mean temperature is below 20 • C.

Fig. 6 .Fig. 6 .
Fig. 6.Comparison of the simulated δ 18 OP −P relation of ECHAM5-JSBACH-wiso (noF ) with ECHAM5-wiso for the resolutions (a) T31L19 and (b) T63L31.For comparison with the observed relation, we use data from those GNIP stations, where the annual mean temperature is above or equal 20 • C. (Please note that the linear fits of ECHAM5-JSBACH-wiso experiment (green line) and ECHAM5-wiso experiment (red line) are almost identical and strongly overlap in the plot.

Fig. 7 .
Fig. 7. Panel (a) shows the annual mean amount of evapotranspiration from land surface, and (b) the fraction of evaporation expressed as percental amount of total evapotranspiration.

Fig. 10 .
Fig. 10.Seasonal cycles of surface temperature T , precipitation amount P , isotopic composition of precipitation δ 18 O P , isotopic composition of soil water δ 18 O ws , depth of soil water bucket reservoir ws, evapotranspiration from land surface ET, and fraction of evaporation E for the locations (a) Reykjavik and (b) Jakarta.The dotted lines represent the observational GNIP values (left = black, right = red).For the simulations the black/red lines represent the simulated T , P , ET, ws and the fraction of evaporation.The simulated δ 18 O values in precipitation and the soil water bucket reservoir are the yellow (noF), green (FE) and blue (FET) lines.

Table 1 .
Pearson Correlation Coefficient R and root mean square error (RMSE) of observed and by-ECHAM5-JSBACHwiso-simulated δ 18 O p (δD P ) values.