Modelling stable water isotopes : Status and perspectives

Modelling of stable water isotopes H2 18O and HDO within various parts of the Earth’s hydrological cycle has clearly improved our understanding of the interplay between climatic variations and related isotope fractionation processes. In this article key principles and major research results of stable water isotope modelling studies are described. Emphasis is put on research work using explicit isotope diagnostics within general circulation models as this highly complex model setup bears many resemblances with studies using simpler isotope modelling approaches.


Introduction
Since several decades, the use of the two stable water isotopologues H 2 18 O and HDO (hereafter simply called stable water isotopes, SWI) has been established as a well-recognized, most useful tool in climate research.It was already in the year 1964, when Dansgaard showed in his pioneering work the general influence of climatic variations on the isotopic composition of precipitation [1].Dansgaard explained successfully how and why in many regions of the Earth the isotopic composition of precipitation is linearly related to the local temperature at the precipitation site (the so-called "temperature effect").He also noticed the breakdown of this effect for tropical regions, where the temperature dependency of the isotopic composition of precipitation is replaced by a much weaker dependency on the total amount of fallen precipitation (the so-called "amount effect").Further regional effects, which may overlap with the temperature effect, are a gradual rainout over land areas (the "continental effect") and a strong dependency between the isotopic composition of precipitation and the altitude of the location in alpine regions (the "altitude effect").
In 1961, Craig suggested to define the isotopic composition of mean ocean water as a useful reference standard for the description of the isotopic composition in any other water sample [2].His approach was accepted by the International Atomic Energy Agency (IAEA), Vienna.Since then the isotopic composition of a sample can be described in a δ-notation as the deviation from the Vienna Standard Mean Ocean Water (V-SMOW).For H 2 18 O, the δ 18 O value of a water sample is calculated as δ 18 O = (([H 2 18 O]/[H 2 16 O]) Sample /([H 2 18 O]/[H 2 16 O]) V-SMOW − 1) • 1000.For HDO, δD is calculated in an identical manner.Craig also analysed in another landmark paper in 1961 the relationship between H 2 18 O and HDO in the ocean and the marine atmosphere [3].He introduced the concept of a "Global Meteoric Water Line" as well as the "Deuterium Excess" (defined as d = 8 • δD − δ 18 O) to describe the relation between the two water isotopologues.Craig discovered that the Deuterium Excess signal is mainly influenced by climatic conditions during evaporation, but not during condensation processes.In the following years, the variability of SWI in precipitation and its dependency on surface temperatures, precipitation amount, evaporation rates and other climate variables was investigated in numerous studies for different regions and time scales (e.g.[4,5]).

EPJ Web of Conferences
The wealth of detailed insights into the various fractionation processes occurring in the Earth's hydrological cycle were possible as in 1961 the IAEA together with the World Meteorological Organisation (WMO) initiated the Global Network for Isotopes in Precipitation (GNIP) [6].Since then, monthly values of surface temperature, precipitation amount as well as the δ 18 O and δD values of precipitation have been reported for more than 700 stations worldwide.While some of the stations have only been active for a few years, the present GNIP station network still includes a substantial number of active stations.Among the longest data series with more than 40 years of continuing measurements are the stations of Vienna, Reykjavik, Ankara and Manaus.Recently, the IAEA started to broaden its monitoring programme by a sibling product to the GNIP database, the Global Network of Isotopes in Rivers (GNIR) database.Until now, another pendant of the GNIP network for regular sampling and measuring the isotopic composition of ocean water masses does not yet, exist.However, since several years Schmidt and co-workers from the NASA Goddard Space Institute, New York, provide a database containing available published δ 18 O measurement of different oceanic samples to the scientific community [7].
Since the spearheading studies of Dansgaard and others, SWI have been shown to be a most valuable tool for paleoclimate studies.They represent a common signal among many different types of climate archives on past changes of the Earth's hydrological cycle.SWI have been measured in ice cores, lake and marine sediments, speleothems, tree-ring and peat-bog cellulose, groundwater and many other paleoclimate archives.However, the physical interpretation of the δ-signal in the various archives is often not straightforward, as the measured signal is a complex integrative composition of different fractionation processes.Especially for terrestrial archives, fractionation during precipitation and evaporation are often mixed with local hydrological (e.g.runoff, infiltration) and biological processes.
Modelling of stable water isotopes within the various parts of the Earth's hydrological cycle might help to improve our understanding of the δ-signal in various archives, and this article strives to describe the key principles and some major research results of SWI modelling.In the next section different conceptual ideas of isotope modelling are presented.We mainly focus on the implementation of SWI into the hydrological cycle of general circulation models (GCM) as this most complex model approach includes many details of isotope modelling also found in simpler isotope modelling studies.The third section describes in brief some key scientific results based on past SWI modelling studies.Of course, such a compendium of achievements will never be complete and is highly biased by a personal choice of results.However, it might nevertheless be helpful informing the readership about the broad range of topics, which contain results of SWI modelling studies.The last section of this article will give a short outlook on potential future progress in the research field of stable water isotope modelling.

Past approaches
SWI modelling studies have nowadays an almost similar long history as the observational studies on stable water isotope variability in the climate system.Dansgaard already used in 1964 a simple Rayleigh-type model approach to explain the various climatic effects on SWI [1].Since then, such Rayleigh-type models have been extensively used in various research projects, e.g.[8,9].To allow a more realistic description of the precipitation formation over Antarctica, Ciais et al. developed the Mixed Isotope Cloud Model (MCIM).This model belongs also to the Rayleigh model family but fully accounts for mixed cloud processes between 0 • C and −30 • C [10].An alternate model approach was presented by Fisher, who constructed a zonally averaged global oxygen isotope model to simulate isotope values in high-latitudinal regions [11].Based on the classical approach of trajectory models to follow the path of an air parcel back and forward in time, Helsen et al. combined the MCIM model with the Royal Netherlands Meteorological Institute trajectory model for detailed studies of the isotopic composition of present-day Antarctic precipitation [12].
The implementation of SWI into general circulation models representing a most complete threedimensional state of the atmosphere was initiated in 1984 by Joussaume et al. using the LMD model [13] and in 1987 by Jouzel et al. using the NASA GISS model [14].For ocean GCM, Schmidt was the first, implementing SWI into the GISS ocean model [15].In contrary to simpler isotope model approaches, the implementation of SWI into general circulation models has the key advantage that the model system is closed and all relevant parameters determining the strength and evolution of isotopic fractionation are simulated by the model itself.This allows in-depth analyses of the dependence of the isotopic composition of various water components in the modelled hydrological cycle on climate-related variables, e.g.surface temperature and precipitation amount.In addition, air masses with different isotopic compositions can be advected and mixed by the dynamical core of a GCM.Furthermore, by using different boundary conditions for GCM model experiments, possible changes of any isotope-climate-relation, e.g. the much-used δ-temperature-relation, can be explored in time and space.Since the pioneering work by Jouzel, Joussaume, Koster and others, water isotopes have been implemented in a growing number of GCM and numerous studies have clearly demonstrated the potential of this approach for an improved interpretation of observed water isotope variability in terms of climate and environmental change (see Section 3).The disadvantages of this modelling technique is the increased complexity of incorporating SWI into a GCM hydrological cycle as well as the required substantial computing power for performing simulations.
In the following text section we will briefly describe the basic principles of modelling SWI within GCM.Many details of the presented concepts are also valid for simpler model concepts.

Implementation of H 2 18 O and HDO within general circulation models
For an atmosphere GCM, so far all existing SWI model developments followed an identical approach, first presented by Joussaume et al. [13] as well as Jouzel et al. [14].The primary idea of this approach is rather simple: While a default GCM just contains the "normal" water tracer H 2 16 O in the model code, for isotope modelling one simply has to triple this code to additionally simulate H 2 18 O and HDO (Figure 1).For many parts of the model code, e.g. the advection of water tracers by the dynamical core of the GCM, the two additional water tracers behave exactly like the default one.However, whenever there occurs a phase change in the model's hydrological cycle, one needs to include additional fractionation effects for both H 2 18 O and HDO.In general, such phase changes of water masses will occur only during evaporation processes from land and ocean surfaces as well as condensation processes during the formation of clouds and precipitation.The later process may include a transition from vapour to liquid (formation of rain), vapour to solid (formation of snow) or even liquid to solid (freezing of rain drops), and vice versa (e.g.re-evaporation of rain below the cloud base).
Although this conceptual approach of SWI modelling is rather simple, the actual incorporation of water isotopes in GCM is very time-consuming and complex.Modern GCM consist of hundreds of subroutines with thousands of lines of code, of which a third or more is related to the hydrological cycle.The physical changes have to be tracked in every single routine, every internal time step, and any of the phases of water (vapour, liquid, solid) implemented in the model code.Within this article, it will not be possible to describe all necessary coding in details.A good description is found, e.g., in [16].However, so far all the different SWI-enabled GCM have used the following common fundamentals: Evaporation from the ocean is in general treated as equilibrium fractionation, with a correction for wind-dependent kinetic effects [17].The required temperature-dependent fractionation strengths for the occurring phase changes are taken from lab experiments by Majoube [18,19].Over land surfaces, water can be uplifted into the atmosphere either by evaporation, too, or by transpiration from plants, which have taken up soil water through their roots.Under steady-state conditions, such plant transpiration is non-fractionating [20].Although water fluxes over land play a key role in the hydrological cycle, so far the soil moisture in isotope-enabled GCM is represented as a simple "bucket" scheme, which no realistic distinction between evaporation and transpiration fluxes.Therefore, the evapotranspiration flux in these models is considered non-fractionating (see [21] for a detailed discussion on this issue).Once water is uplifted as vapour into the air, it is transported horizontally and vertically without any fractionation effect.At some time steps later, when the environmental conditions, like relative humidity and temperature, are appropriate, condensation and the formation of precipitation within a model grid box might occur.Most GCM typically divide these processes into stratiform (large-scale) and convective cloud and precipitation formation.For condensation processes, in general an equilibrium fractionation is assumed, too.At low temperatures, e.g. during the formation of snow over the polar regions, an additional kinetic correction accounting for the slower diffusivity of heavier isotopes has to be applied [8].After the precipitation is formed, it typically falls through a zone of under-saturated air below the cloud base, which leads to non-negligible re-evaporation of the precipitation.As a consequence, a partial re-equilibration of rain droplets with surrounding water vapour occurs.The strength of this re-equilibration is mainly determined by the size of the droplets as well as their velocity of fall [21].
Although most isotope-enabled atmosphere GCM use very similar formulations of these key fractionation processes, the underlying physical parameterisations of the GCM (e.g.convective and stratiform cloud schemes, partitioning of different water reservoirs on land surfaces) are very different.Thus, the implementation of SWI into different atmosphere GCM enables an additional in-depth evaluation and comparison of the physical representation of the hydrological cycle within the different models.
For an ocean GCM, the principal implementation of SWI is slightly simpler than for an atmosphere GCM.In the ocean, no additional fractionation of the liquid water masses occur and both H 2 18 O and HDO can be treated as a passive tracer [15].Its concentration is changed by advection and convection of different water masses, only.The only major fractionation process to consider is the evaporation of surface water into the atmosphere, which leads to a slight enrichment of heavy SWI in the near-surface layer of the oceans.In high-latitude regions, minor additional fractionation of SWI occurs during the formation of sea ice.

Alternate 3-D water isotope modelling concepts
In addition to general circulation models, SWI were recently implemented in two regional climate models [22].The principles of implementation are identical to the GCM setup, but a regional climate model enables simulating the isotopic composition of different water masses for a selected area with a much finer spatial resolution (typically one order of magnitude finer as compared to a regular GCM simulation).However, a SWI-enabled regional model always needs suitable lateral isotopic boundary values, which have to be derived from an isotope-enabled GCM or an alternate data source.The advantage of the finer spatial resolution of regional models is accompanied by higher computational demands, which typically limits the utmost simulation period of such regional models to a few decades.Since a number of years, Earth System Models of Intermediate Complexity (EMIC) have become more and more popular and are nowadays an own class of model types.The conceptual idea of an EMIC is reducing the complexity of model physics found in a typical GCM to save computational costs.Such models can then either be run for a much longer simulation period than a GCM, or a multiensemble of shorter EMIC simulations might be used to study the effects of various parameterizations on selected climate periods and environmental processes.To our knowledge, water isotopes have so far be implemented into one EMIC, the CLIMBER-2 model, only.Roche et al. used SWI diagnostics within CLIMBER to analyse the oxygen-18 content of the water masses in the oceans for present-day and glacial maximum climate conditions [23].

Present achievements in modelling water isotopes
As explained in the last section, over the last three decades SWI have been built into a number of different models, ranging from simple Rayleigh models to complex GCM.For the technically most challenging group of isotope-enabled GCM, at present at least a dozen different models exist (Table 1).In this section we will give a brief overview of some key simulation results of these models.
Usually, any new model enhancement or change in a model setup has first to be evaluated against present-day observational data.For SWI, there is no exception of this scientific practice, and with the GNIP measurements an appropriate dataset for such model evaluation exists.It comes with no surprise that all published SWI model results do agree rather well with the GNIP data on a global scale.If a model would fail this test, it would probably not have passed any peer-review process for publication.Exemplarily for the simulation of the isotopic composition of precipitation using a SWI-enabled atmosphere GCM, a 20-year long simulation of the ECHAM4 model is shown (Figure 2).The horizontal model grid size for this simulation was about 1 • × 1 • (spectral T106 resolution).A comparison against a global compilation of GNIP measurements of δ 18 O in precipitation during the years 1961-1999 reveals a good agreement between simulation values and observational data.The ECHAM4 model clearly maps all known key effects of δ 18 O in precipitation: the general temperature-dependent depletion going from low-latitudes to high-latitudes, the very depleted isotope values observed in Greenland and East Antarctic snow, a typical continental effect over Western Europe and the North American continent, as well as altitude-related isotope depletion in alpine regions like the South American Andes.In more quantitative analyses, it has been shown that the ECHAM4 model is capable of simulating the observed δ 18 O-T Surface -gradient for various mid-to high-latitude areas as well as for the Polar Regions in a decent manner [21,24].The same is true for the relation between δ 18 O in precipitation and total amount of precipitation in tropical regions.Similar results have been reported for the majority of the models listed in Table 1.After the successful implementation of SWI into various GCM, different research studies have aimed on an improved understanding of the influence of various climate parameters and environmental processes on the isotopic composition of modern precipitation.Among others, it has been shown that a high variability of isotopes in precipitation exists on short-term synoptical time scales, which goes far beyond the variability imprinted in the monthly GNIP data set [21,25].On interannual time scales, precipitation variations might stronger influence the δ-signal than temperature variations [26,27].Furthermore, an imprint of interannual climate oscillations like ENSO and NAO on the δ 18 O signal in precipitation can be detected in various regions, both in observational data and simulation results [28][29][30].Some of the isotope-enabled GCM have performed so-called "tagging" experiments where the atmospheric transport of water masses from pre-defined evaporation regions can be traced [31].Such tagging experiments allow an explicit decomposition of precipitated water at a certain location into the contributions of different source regions and the influence of the different sources on the overall δ-signal in precipitation [24,32,33].
For SWI-enabled ocean GCM, comparisons with data showed that regional linear relationships between δ 18 O of seawater and salinity could be reasonably captured by the model.Another important result was the finding that the seasonal and long-term δ 18 O-salinity-relation at a fixed point is not generally equal to the nearby spatial relation at a fixed time [15].
Since the first isotope-enabled GCM became available, another focus of SWI simulations has been on various climates of the past.The inclusion of SWI within GCM for this purpose can be described as a typical example of "forward proxy modelling".It is possible to directly compare GCM simulation results with the measured isotopic data from climate archives, without any reconstruction of climate variables prior to the simulation-data comparison.The latter method is termed as "inverse proxy modelling" and is still the standard procedure for many other paleoclimate modelling studies.Like numerous regular GCM paleoclimate simulations, the first SWI-enabled GCM paleoclimate simulations have focussed on the climate optimum during the mid-Holocene (6 K, 6,000 years B.P.) and the climate of the Last Glacial Maximum (LGM, 21,000 years B.P.) Key results of the performed studies include general analyses of global δ-changes in precipitation under 6 K and LGM climate conditions [34,35], a most-likely explanation of the failure of SWI-based thermometry for glacialinterglacial temperature changes on Greenland [36,37], as well as an affirmation of isotope-based estimates of glacial-interglacial temperature changes derived from East Antarctic ice cores [38].Apart from climate simulations of the mid-Holocene and the LGM, SWI-enabled GCM simulations have also been performed and analysed for several other climate periods, like the middle Cretaceous [39], and warmer-than-present interglacials [40].It is certainly safe to say that these numerous studies on the isotopic signature of precipitation during past climate periods have considerably helped improving the interpretation of the observed isotope changes in various paleoclimate archives in terms of temperature and precipitation amount changes.
To foster the assessment of different SWI-enabled GCM studies and homogenize their evaluation against observational records, e.g., the GNIP data set, an isotope GCM inter-comparison study was initiated in 2004.The first Stable Water Isotope Intercomparison Group (SWING) workshop was hosted and funded by the Isotopic Hydrology Section of the IAEA, Vienna.As a result of this workshop, for the very first time three different SWI-enabled atmosphere GCM (ECHAM4, GISS-E, MUGCM) conducted two present-day climate simulations in an absolute identical manner.On a global scale, all three models simulate the pattern of d18O in precipitation in good agreement with the GNIP observations; however, on regional scales substantial differences between the three different GCM exist (Figure 3).It was shown that the presently available observational database of isotopes in precipitation is not always sufficient to assess the identified model differences in a more quantitative manner.Analyses also revealed that the three SWI-enabled GCM divergence more strongly for the simulation of the Deuterium Excess signal than for δ 18 O in precipitation (Figure 4).In the year 2008, a second SWING workshop was hosted by the IAEA to re-define SWI simulations suitable for an enhanced new model intercomparison study.At this workshop, 12 model groups participated (Table 1) and all models performed a common simulation for the present-day climate of the last 50 years.Three of the models (ECHAM4, GSM, LMDZ) have been nudged to reanalyses data.In these cases, the atmospheric circulation is forced to reproduce the observed weather while the (isotopic) water cycle is left unforced, which facilitates the direct comparison between observation and model output.In connection with the 3rd Paleoclimate Modelling Intercomparison Project (PMIP3), it is envisaged defining another SWING experimental setup, where participating SWI-enabled GCM will be compared under different paleoclimate conditions [22].

Future developments
It seems very likely that the future progress in SWI modelling will be mainly driven by recent scientific advances in two different fields: the continuing trend of building more integrative models of the Earth System, as well as some major technological advancements in the fields of water isotope measurements.
Over the last two to three decades, the merely physical-based atmospheric and oceanic GCM setups have more and more evolved into complex integrative models of the Earth System including also biological and chemical processes relevant for climate change.Accordingly, nowadays state-of-the-art model approaches are named Earth System Models (ESM).Apart from atmospheric and oceanic physical processes, such ESM may also include, among others, chemical processes within the troposphere and stratosphere, a dynamical representation of land vegetation processes, multi-level soil moisture schemes, as well as key biogeochemical processes to simulate a closed carbon cycle within the ESM.Dynamical ice models are not yet part of most ESM, but will most certainly be included in future model releases, too.For the modelling of SWI, a few research groups have only recently started to incorporate H 2 18 O and HDO in such ESM [41,42].However, it seems safe to assume that SWI will be incorporated into several other ESM in the very near future.With a wealth of new observational data sets (see below) validation of the present-day water cycle within an ESM might become even more relevant for its general evaluation.For example, first attempts using SWI for an improved description of tropical convection processes have been started [43].For paleoclimate simulations using ESM, a combination of forward modelling of different climate proxies, e.g.water (H 2 18 O, HDO) and carbon ( 13 C, 14 C) isotopes, seems feasible.With the anticipated further increase of computational speed and resources, a dynamic ESM simulation covering a full glacial-interglacial cycle should be possible soon.Such simulations will reveal further insights into temporal and spatial isotope-temperature-relationships and the effects of climate change imprinted in the various δ-paleorecords.
While for many years the isotopic composition of a water sample could only be measured by high-precision mass spectrometry, the situation has drastically changed during the last years.With the development of Cavity Ring-Down Spectroscopy methods it is now possible to measure SWI in liquid water samples with the required high accuracy in much faster time.Furthermore, isotopic measurements of water vapour samples are possible with the same technique.Thus, in the future any isotope model evaluation can not only be based on the GNIP data set, but might also be evaluated against a wealth of new data of H 2 18 O and HDO in the atmosphere.These new ground-based measurement techniques of SWI are accompanied by new isotope data derived from space-born measurements.At present, two such data sets are available [44,45].While the Tropospheric Emission Spectrometer dataset measured aboard the Aura spacecraft contains measurements of the period October 2004 to March 2005 [44], the SCIAMACHY dataset from the ENVISAT satellite includes the integrated δD signal of the total atmospheric water column on a global scale for the period 2003-2005 [45].Once again, it is safe to assume that both the newly developed ground-based isotope measurements and the various upcoming satellite products will enable a much-improved SWI model evaluation.
To conclude, the ongoing fast development of more integrative ESM as well as the imminent sum of newly derived observational data will most certainly enable a more quantitative assessment of different SWI model approaches in the near future.This represents an exciting challenge for the related scientific community and will further both broaden and deepen our insights in the various isotope-related processes of the Earth's hydrological cycle.

Fig. 1 .
Fig.1.The principal implementation of stable water isotopes into the hydrological cycle of an atmospheric general circulation model.Two additional tracers, H 218 O and HDO, are introduced into in parallel to the default model's water H 216 O.In addition, isotope fractionation effects during any phase change of a water mass, e.g. during evaporation and condensation processes, have to be implemented in the model code.

Fig. 2 .
Fig. 2. Comparison of δ18O in precipitation measurements archived in the GNIP database for the period 1963-1999 versus an ECHAM4 T106 simulation of δ18O in precipitation covering the period 1979-1998.Long-time annual mean values of the monthly GNIP data and ECHAM5 simulation are plotted.For comparison reasons, the GNIP data has been spatially interpolated in the vicinity the different station values.

Fig. 3 .
Fig. 3. Comparison of global maps of mean δ 18 O in precipitation based on GNIP measurements and simulation results of three different SWI-enhanced atmosphere GCM (ECHAM4, GISS-E, MUGCM).All three simulations were forced with identical present-day boundary conditions to enable an in-depth comparison of the different GCM results.

Fig. 4 .
Fig. 4. Comparison of zonal means of δ 18 O in precipitation and Deuterium Excess values d based on GNIP measurements and simulation results of three SWI-enhanced atmosphere GCM (ECHAM4, GISS-E, MUGCM) taking part in the first phase of the SWING intercomparison project.

Table 1 .
[22]view of SWI enabled GCM and institutes participating in the 2nd phase of the SWING project (table adapted from[22]).