Water isotope variations in the global ocean model MPI-OM

. The stable water isotopes H 182 O and HDO are incorporated as passive tracers into the oceanic general circulation model MPI-OM, and a control simulation under present-day climate conditions is analyzed in detail. Both δ 18 O and δ D distributions at the ocean surface and deep ocean are generally consistent with available observations on the large scale. The modelled δ D- δ 18 O relations in surface waters slightly deviates from the slope of the global meteoric water line in most basins, and a much steeper slope is detected in Arctic Oceans. The simulated deuterium excess of ocean surface waters shows small variations between 80 ◦ S and 55 ◦ N, and a strong decrease north of 55 ◦ N. The model is also able to capture the quasi-linear relationship between δ 18 O and salinity S , as well as δ D and S , as seen in observational data. Both in the model results and observations, the surface δ - S relations show a steeper slope in extra-tropical regions than in tropical regions, which indicates relatively more addition of isotopically depleted water at high latitudes.


Introduction
The stable water isotopes H 18 2 O and HDO (D = deuterium) have been considered as a useful tool in ocean and climate research since the pioneering work by Epstein and Mayeda (1953), Dansgaard (1953) as well as Friedman (1953). It was recognized that the observed ratios of water isotopes in different water samples are affected by temperature-dependent isotope fractionation occurring at any phase transition within the Earth's hydrological cycle. In addition, different freshwater fluxes into the ocean (e.g. precipitation, or river discharge) have distinctive isotopic signatures, which enables the use of water isotopes as tracers of different oceanic water masses, comparable but not identical to salinity (Craig and Gordon, 1965).
In modern ocean research, measurements of water isotopes have been used to obtain information on the origin and mixing pattern of water masses (e.g. Östlund and Hut, 1984;Fröhlich et al., 1987;Gat et al., 1996;Khatiwala et al., 1999). Oxygen isotopic ratios are also considered as an important proxy in many climate archives (e.g. marine and lacustrine sediment cores, ice cores, etc.), which among others, are used for reconstructing paleo-temperature (e.g. Broecker, 1986), paleo-salinity (Duplessy et al., 1991;Sarnthein et al., 1992), and past sea level variations (Sowers and Bender, 1995;Broecker and Henderson, 1998).
For quantitative analyses, the isotopic composition of a water sample is commonly described in a permil-based delta notation (δ 18 O, δD). The interrelationship between hydrogen and oxygen isotope ratios is defined as deuterium excess d-excess (d-excess = δD-8*δ 18 O, Dansgaard, 1964). It is an indication of kinetic (non-equilibrium) fractionation effects occurring during evaporation of water from oceanic areas (Dansgaard, 1964). The d-excess records from polar ice cores have been used to reconstruct ocean surface temperature at the evaporative source region (e.g. Ciais et al., 1994;Vimeux et al., 2001;Masson-Delmotte et al., 2005;Stenni et al., 2010). As the d-excess value is determined by the air-sea interaction processes, the ocean surface conditions strongly affect d-excess values in vapor and rainfall (Uemura et al., 2008). Consequently, the d-excess variations in surface ocean water should go in parallel with d-excess variations in vapor and rainfall as well.
Since several years, general circulation models (GCMs), which include explicit stable water isotopes diagnostics, are established tools to investigate the relationship between water isotopes and various climate variables. The water isotopes 810 X. Xu et al.: Water isotope variations in the global ocean model MPI-OM have been incorporated in both atmospheric GCMs (first realizations by e.g. Joussaume et al., 1984;Jouzel et al., 1987;Hoffmann et al., 1998) and oceanic GCMs (Schmidt, 1998;Paul et al., 1999;Wadley et al., 2002), as well as in coupled atmosphere-ocean models (Schmidt et al., 2007;Tindall et al., 2010). These previous modeling studies stressed the importance of understanding the spatial and temporal isotopic variations for a quantitative interpretation of their relationship to climate changes. For oceanic GCMs, in particular, the studies showed also the potential of δ 18 O to characterize different water masses.
Here, we report first results with a new isotope-enhanced version of the ocean general circulation model MPI-OM (Marsland et al., 2003). Our study focuses on a simulation of the present-day global oceanic distribution of δ 18 O and δD as well as the deuterium excess. The model results are analyzed for different oceanic basins to understand the processes that lead to the simulated isotopic distribution of both oxygen-18 and deuterium at the ocean surface. A comparison of the simulation results with available observations examines the model performance in capturing the main characteristics of the present-day oceanic water isotopic distribution as well as the relationship of salinity versus δ 18 O and δD in different ocean basins.

Ocean model
We employ the Max-Planck-Institute Global Ocean/Sea-Ice Model MPI-OM (Marsland et al., 2003), which solves the primitive equations on an Arakawa C-grid (Arakawa and Lamb, 1977) while the vertical discretization is on z-level. The model has a free surface and includes subgrid-scale parameterizations for the bottom boundary layer flow across steep topography, horizontal and vertical viscosity, vertical and isopycnal diffusivity, and convection. As a sea ice component, a viscous-plastic rheology sea ice model (Hibler, 1979) is embedded. This sea ice model includes thermodynamic sea ice growth and melt, as well as the thermohaline coupling to the ocean model by brine rejection.
The model is implemented on a bipolar orthogonal spherical coordinate system, with the poles located over Antarctica and Greenland, respectively. The second location avoids a grid singularity in the Arctic Ocean while the horizontal resolution is high in the deep water formation regions of the Arctic and the northern North Atlantic Oceans.
For this study, the MPI-OM grid configuration GR30 with a formal horizontal resolution of 3 • × 1.8 • , 40 unequal vertical layers, and a time step length of 8640 seconds is chosen. Initial conditions for marine temperature and salinity are interpolated from climatological fields given by Steele et al. (2001). Sea surface salinity is restored with a time scale of 39 days.
The atmospheric forcing is derived from output of an atmospheric general circulation model (AGCM) presentday control simulation in a T31 spectral resolution (about 3.75 • × 3.75 • ), with 19 vertical levels (ECHAM5-wiso; Werner et al., 2011). Mean daily values of heat, freshwater and momentum fluxes at the air-sea interface are obtained from a 10 yr simulation period. All forcing fields from the ECHAM5-wiso simulation are bilinearly interpolated to the MPI-OM grid. The MPI-OM present-day climate control simulation has been run for 3000 yr into a quasi-steady state. For the analyses we used the mean state of the last 100 simulation years, only.

Modelling of isotope tracers H 18 2 O and HDO
The stable water isotopes H 18 2 O and HDO are incorporated as conservative passive tracers into MPI-OM. They are fully advected and mixed by the OGCM. The isotopic variations occurring in MPI-OM primarily depend on the temperaturedependent isotope fractionation during evaporation, on the isotopic composition of freshwater fluxes entering the ocean, and on oceanic advection and mixing of different water masses.
For an OGCM isotope simulation, the required isotopic composition of evaporation (E) and precipitation (P ) fluxes can either be taken from observational data or obtained from suitable atmospheric model output. Several previous studies applied a regression formula (Gat and Gonfiantini, 1981) to calculate the isotopic composition of precipitation (Schmidt, 1998;Wadley et al., 2002). This global regression is based on observed local temperatures and precipitation amounts over continents and some island stations. The high latitudes are underrepresented by this empirical approach, and the tropical amount effect is not well captured (Rozanski et al., 1993).
In this study, the daily isotopic content of P and E fluxes over ocean surface stem from the same ECHAM5wiso simulation as used for obtaining the default freshwater, heat and momentum fluxing forcing. This approach avoids the systematic uncertainties mentioned above and ensures maximum consistency between the prescribed climatological forcing and the related isotopic fluxes. ECHAM5wiso is a derivative of the ECHAM5 AGCM (Roeckner et al., 2003) enhanced by explicit water isotope diagnostics (Werner et al., 2011). The prescribed δ 18 O of sea surface water in ECHAM5-wiso is interpolated from the NASA gridded data set (LeGrande and Schmidt, 2006), and the δD values are obtained by assuming an oceanic d-excess value of 0 ‰ (δD = 8 * δ 18 O, Supplement Fig. S1). On a global scale, the ECHAM5-wiso simulation results are in good agreement with available observations of the isotopic composition of precipitation from the Global Network of Isotopes in Precipitation (GNIP), both on an annual as well as a seasonal time scale (Werner et al., 2011).
For the phase transition of water occurring during the formation and melting of sea ice, no change of δ 18 O and δD are assumed in our MPI-OM setup, as the isotopic fractionation between liquid water and ice is small compared to other fractionation processes (Craig and Gordon, 1965;O'Neil, 1968;Souchez et al., 2000).
The isotopic composition of river runoff is another boundary condition, which influences the marine isotopic values, particularly in coastal regions. In this study, the amount and isotopic composition of continental runoff are calculated from a bucket model, which drains the continental freshwater fluxes simulated by ECHAM5-wiso following the topographic slope. The discharge of continental net precipitation (P -E) towards the ocean depends on the slope of the underlying topography, and the calculated total river runoff is based on the assumption that the global water cycle is balanced, i.e. all net precipitation over land surfaces is transported to the ocean. In our procedure, lakes are not present and the hydrological cycle is closed. This may introduce small distortions in magnitude and location of the simulated river runoff compare to observations.
For the initial distribution of H 18 2 O and HDO in the ocean, a homogenous setup has been chosen with all values set to present-day δ 18 O and δD values of 0 ‰ with reference to the Vienna Standard Mean Ocean Water (Baertschi, 1976;Dewit et al., 1980). No surface restoring of water isotopes is applied.

Observational database of δ 18 O and δD
The Global Seawater Oxygen-18 Database (Schmidt et al., 1999) contains a collection of over 26 000 observations since about 1950, which offers a unique opportunity to compare the observed and modeled oceanic δ 18 O values. For the comparison with our model results, all data entries of the top 10 m-layer are taken as representative ocean surface water δ 18 O values. All measurements between 100 m and bottom are used for model-data comparison at subsurface and deep oceans.
Although δD observations of ocean waters are not as widespread as δ 18 O, there are some measurements available to validate our simulations. For this study, we consider global δD measurements from the GEOSECS expeditions (Östlund et al., 1987), which are complemented with observations from the Northwest Pacific (Horibe and Ogura, 1968), from the tropical West Pacific (Aharon and Chappell, 1986), the Indian Ocean (Duplessy, 1970;Delaygue et al., 2001), the Arctic Ocean (Friedman et al., 1961), the Baltic Sea (Ehhalt, 1969;Fröhlich et al., 1988), the Mediterranean (Gat et al., 1996), the Gulf of Mexico (Yobbi, 1992), and the Weddell Sea (Weiss et al., 1979). In the same way as applied for δ 18 O, the δD measurements above 10 m and below 100 m have been chosen. The chosen surface data have been averaged over the upper 10 m.

Simulated present-day distribution of water isotopes (H 18 2 O and HDO)
For a first assessment of our simulation results, we evaluate the global zonal-mean isotopic composition of δD and δ 18 O in ocean surface waters (Fig. 1). Due to the lager mass difference of the two stable isotopes relative to the mass of the element, zonal deuterium exhibits much larger absolute variations (−15 ‰ to +5 ‰) in the isotopic composition of surface waters as compared to δ 18 O (−3 ‰ to +1 ‰). However, as the same boundary fluxes (precipitation, evaporation, and river runoff) determine both H 18 2 O and HDO in ocean surface waters, strikingly similar features of the zonal mean isotopic composition can be seen in most regions. Both δ 18 O and δD decrease from mid-to-high latitudes as a result of the latitudinal temperature decrease and the rainout effect, similar to the trends observed in meteoric water samples (e.g. Craig and Gordon, 1965;Criss, 1999). There is a slight increase of the water isotopes concentration from equatorial regions to the subtropics, which can be related to a relative excess in evaporation over subtropical oceans. The most enriched areas are located around 30 • N, where the zonal mean δD and δ 18 O values reach almost +5 ‰ and +1 ‰, respectively. The most depleted regions are found north of 65 • N, with δD values less than −20 ‰ and δ 18 O values less than −1.5 ‰. The isotopic values in the Northern Hemisphere exhibit a much stronger heterogeneity as compared to the Southern Hemisphere.
For further analyses, the simulated global δ 18 O and δD distributions at the sea surface as well as the observations are shown in Fig. 2 (Schmidt et al., 1999). respectively), where evaporation is the dominating influence on the isotopic composition of seawater.
The spatial structures of δ 18 O and δD simulated by MPI-OM are similar to the available observations (Figs. 2b and 3b). Both model results and observational data show strongly depleted areas (δ 18 O < −5 ‰, δD < −20 ‰) in coastal regions with continental river discharge. For a more quantitative model-data comparison, we average all observations within a specific grid cell and compare the mean value with corresponding MPI-OM result (Fig. 4). Analyses of the correspondence between simulated and observed fields indicate a good agreement between these two fields.   similar features as δ 18 O, but larger variation. Its distribution should be very similar to δ 18 O after an appropriate scaling. Performing the same averaging method of observational values within different model grid cells as for the surface water data, we find that the model results below 100 m also agree well with the observations. The Rs are 0.76 and 0.82 for δ 18 O and δD, and the NRMSEs are 4.6 % and 13.2 %, respectively (Fig. 6).

Simulated present-day δD-δ 18 O relationship and d-excess variations
Based on precipitation measurements, a linear relation between δD and δ 18 O in precipitation was first reported by Dansgaard (1964), with δD = 8*δ 18 O. This relation can be explained by a simple Raleigh-type model of isotopic fractionation occurring during evaporation and condensation processes. As illustrated in Fig. 7, a similar linear relationship between δD and δ 18 O at sea surface is found in most basins. The simulated deuterium and oxygen-18 composition in both basins are indeed highly correlated (r > 0.98), and the slopes of the linear fit between δD and δ 18 O slightly deviate from the factor 8 at most basins (Atlantic Ocean: 7.76;  Based on the simulated water isotopic compositions, the modeled ocean surface d-excess values are calculated. For ocean regions, the d-excess in the surface waters may be affected by evaporation, precipitation, as well as mixing of different water masses. We compile existent observational data where both δD and δ 18 O measurements exist at the same location and calculate for these location the d-excess of ocean surface waters (Fig. 8). The modeled δD-δ 18 O relationship agrees with the observations, but a clear negative offset is detected for the simulated d-excess values. Both observations and simulation show similar ranges of d-excess variations with an increase in d-excess for water masses with stronger δ 18 O-depletion.
The simulated global zonal mean d-excess values of sea surface waters vary between −2 ‰ and −13 ‰, with  rather small variations (−4 ± 1 ‰) between 80 • S and 55 • N (Fig. 9). These negative values are in general agreement with the enriched (positive) deuterium excess values of atmospheric water vapor simulated (and observed) in the marine boundary layer. A sharp decrease is detected at around 55 • N, where the d-excess drops from approximately −3.5 ‰ to −13 ‰.
In Fig. 10, we illustrate the spatial variations of the dexcess values. There are few marginal seas (Baltic Sea, Bay of Bengal, James Bay) which exhibit a positive d-excess, most likely strongly influenced by river runoff. The variances of d-excess in major basins are small, but there are significant differences between the different ocean basins. The Pacific and Indian Oceans are relatively enriched in d-excess as compared to the Atlantic Ocean, and the most depleted regions are located in the Arctic Oceans. The water transported by the Labrador Current and the East Greenland Current contains much lower d-excess values than water of the Gulf Stream and the North Atlantic Drift. For the Arctic Ocean, the inflow from the North Atlantic and Pacific preserves its dexcess value from its original region. The Pacific water flow through the Bering Strait has higher d-excess, and the North Atlantic current brings relatively lower d-excess values. In the Southern Hemisphere, the Antarctic Circumpolar Current traps the depleted water in the Southern Ocean.

Simulated present-day relations between H 18 2 O, HDO and salinity
Because similar physical processes (precipitation, evaporation, and river runoff) affect the salinity and the isotopic composition of ocean waters, a positive relationship between salinity (S) and water isotopes can be expected, particularly for ocean surface waters.
As illustrated in Fig. 11, the model is able to represent the general features of the observed δ 18 O-S relation. The modelled δ 18 O-S slope at tropical oceans is  For δD, there is a limited number of measurements, and even less investigation on the δD-S relation (Baltic Sea: Fröhlich et al., 1988; Mediterranean Sea: ∼1.72 ‰ PSU −1 , Gat et al., 1996). This sparse data set is not enough for a robust data-model comparison. Therefore, we have selected the same grid cells used in the δ 18 O-S analyses to determine the δD-S relations in tropical and extra-tropical oceans (Fig. 12). In agreement with the similarity of the basin-scale variations of deuterium and oxygen-18, we find coincident relations of δ 18 O and δD with salinity. For δD, a slope of ∼1.06 ‰ PSU −1 is found in tropical regions, and a much steeper slope of ∼4.09 ‰ PSU −1 is simulated in extratropical regions. At both tropical and extra-tropical oceans, δD exhibits larger absolute variation compare to δ 18 O with the same range of salinity, and shows relatively weaker correlation to salinity.

Discussion
We implemented the passive isotopic tracers H 18 2 O and HDO into the ocean model MPI-OM. In most regions, the simulated water isotope composition of surface waters agrees well with observations. The correlation between model results and observations is higher for δD (R = 0.96) than for δ 18 O (R = 0.75). However, this is probably biased due to the smaller number of δD observations compared to δ 18 O. If the model-data comparison for δ 18 O is restricted to the same locations used in the δD analyses, a similar high correlation coefficient (R = 0.97) is calculated. The largest differences are found in some marginal seas, where river runoff has a strong effect, such as in the Baltic Sea, the Kara Sea and Laptev Sea. The deviations between data and model are likely due to the unresolved scales and processes at those marginal areas, as a result of the coarse resolution in our model simulation. Because we excluded the small fractionation during sea ice formation, they will be slightly over-depleted in isotopes at sea ice formation area. But this effect is rather minor (Craig and Gordon, 1965;O'Neil, 1968;Souchez et al., 2000). Vertically, the number of observations from deeper water layers is insufficient to draw a precise global picture of the vertical distribution of water isotopes. However, the simulated vertical values of δD and δ 18 O agree also well with the observations. There are a small number of δD values completely off the 1:1 line in Fig. 6. These values are located at the costal line of the Weddell Sea. The depletion in our simulation is probably mainly due to the runoff calculation, which considers all the snowfall as river runoff into the ocean according to the mass balance. The main discrepancies in δ 18 O as seen in the model-data comparison are found mainly in the Gulf of St. Lawrence. This mismatch may be due to our coarse model resolution, which is unable to resolve the physical processes in this region. The modeled water isotope distribution at Atlantic meridional section is consistent with the different water masses. Our model results compare reasonably well with previous results from other isotope included OGCM model studies, either by ocean only (Schmidt, 1998;Paul et al., 1999;Wadley et al., 2002) or coupled versions (Schmidt et al., 2007).
Since precipitation and evaporation are the main water fluxes affecting the isotopic composition of ocean surface water, a similar linear relation between δD and δ 18 O as the global meteoric water line (Craig, 1961) is expected. Our simulation represents this similar relation in most oceans and consistent with the observed slope 8.1 ± 0.4 (Friedman, 1953) in the ocean as well. In addition, the unique δDδ 18 O relation in the Mediterranean Sea found by Gat et al. (1996) is not reproduced by our model. In agreement with recent observations (Cox et al., 2011), the simulated Mediterranean relation follows the general δD-δ 18 O relation. The reasonable δD-δ 18 O relation enables a realistic d-excess surface field, which is necessary to model d-excess. However, there are only a few studies on the d-excess in ocean waters, which makes it rather difficult to evaluate the modelled d-excess variations in sea surface waters. Based on the approach that the atmospheric d-excess signature is dominated by non-equilibrium isotope fractionation during evaporation processes from marine surface waters (Dansgaard et al., 1964;Gat et al., 1994), the remaining ocean surface water should have a negative d-excess value. In our simulation results, we find such negative d-excess values in most ocean regions. Water masses with positive d-excess values are found in regions that are controlled by the river runoff, preserving a positive d-excess signal found in precipitation. Compared to the available observations, our model simulates in general too negative d-excess values. This mismatch might be related to our chosen model setup: The used isotope forcing fluxes of precipitation and evaporation stem from an ECHAM5-wiso simulation, where a constant d-excess value of 0 ‰ was assumed for all ocean surfaces (Werner et al., 2011). This inconsistency can be only solved in a fully coupled ocean-atmosphere system and should result in slightly higher d-excess values of ocean surface waters.
For the relation between δ 18 O and salinity S in marine surface waters, different gradients have been obtained in various model studies. These studies suggest an extra-tropical region of 0.5-0.6 ‰ S −1 and a tropical slope of 0.2 ‰ S −1 (Schmidt, 1998;Schmidt et al., 2007;Delaygue et al., 2000;Wadley et al., 2002). Our results are in agreement with these findings. Both deuterium and oxygen-18 exhibit a linear relationship with salinity, and both δ 18 O-S and δD-S gradients vary when different ocean basins are chosen for the calculation. Thus, for a more precise interpretation of δ-changes in terms of salinity variations, basin-scale relationships should be applied. However, the surface salinity is not a free variable in our model setup due to the surface salinity restoring, which is necessary to maintain the general MPIOM ocean state in agreement with the observations. Our relation between δ 18 O and S is therefore not purely determined by the model, but also by the salinity data used for restoring. These inconsistencies in the δ-S relation will be overcome in future work using a fully coupled ocean-atmosphere set up.

Conclusions
In this study, the water isotope tracers H 18 2 O and HDO have been successfully implemented into the ocean GCM MPI-OM to explicitly simulate the present-day oxygen and hydrogen isotopic composition of seawater. In addition, a secondorder isotope parameter, the deuterium excess of sea surface water, has also been analyzed on a basin scale. Our model captures the observed basin-scale features and meridional gradients of surface ocean water. The model simulates a more enriched Atlantic in comparison to the Pacific Ocean, the Mediterranean Sea as the most enriched oceanic basin, and the Arctic estuaries as the isotopically most depleted regions. Regional deficiencies showing up in certain marginal sea are related to relatively coarse horizontal resolution of the model setup. Our model also successfully simulates the observed vertical distribution and water-isotopic composition of Atlantic water masses. Moreover, the simulated d-excessδ 18 O relations agree with the observations, with an increase of the d-excess for water masses stronger depleted in δ 18 O. A reverse relationship is detected at high latitudes, where the δD-δ 18 O relations show a steeper slope (>8). The simulated δ 18 O-salinity relationships of tropical and extra-tropical oceans are also in good agreement with observations, with a flatter δ 18 O-S slope at low latitudes and steeper slope at high latitudes. The modelled δD-salinity relationship shows similar features as δ 18 O.
The first evaluation of the newly developed stable water isotope module within the MPI-OM ocean model has been focused on present-day climate conditions. As a logical next step, the model will be applied to different paleoclimate conditions in a fully coupled setup to improve our understanding of past marine isotopic changes and its use for paleo-climate reconstructions.