Stable water isotopes in the MITgcm ( checkpoint 64 w )

Völpel et al. present in their manuscript first results of a newly implemented stable water isotope (SWI) diagnostics within the ocean GCM MITgcm. Their evaluation of this model enhancement focuses on modelling results of H218O in a simulation under preindustrial climate conditions. This evaluation contains both a model-data comparison using measurements of δ18O in seawater and in planktic foraminifera as well as a brief analysis of the simulated δ18O-salinity relationship in different water bodies.


Introduction
Stable water isotopes (H 2 16 O, H 2 18 O and HD 16 O = HDO) are widely used tracers of the hydrological cycle (Craig and Gordon, 1965;Gat and Gonfiantini, 1981) and can be used to determine the origin and mixing pattern of different water masses (e.g.Jacobs et al., 1985;Khatwala et al., 1999;Meredith et al., 1999).Due to differences in their physical and chemical properties, stable water isotopes undergo fractionation processes at any phase transition within the hydrological cycle (Craig and Gordon, 1965).This leads to distinctive isotopic signatures for different freshwater fluxes, which are commonly expressed as δi (i = 18 O or D) with reference to the Vienna Standard Mean Ocean Water (VSMOW) standard and given as: where R is the ratio of the abundance of the heavier water isotope H 2 18 O (HDO) to the abundance of the lighter isotope H 2 16 O and   = 2005.2• 10 −6 for δ 18 O (Baertschi, 1976) and 155.95 • 10 −6 for δD (de Wit et al., 1980).
Stable water isotopes have been used as an important proxy in a wide range of climate archives e.g. to infer past temperatures from ice cores (Johnsen et al., 2001) or the amount of monsoonal rainfall from speleothems (Fleitmann et al., 2003).As an indirect record, stable water isotopes are preserved in carbonates (CaCO 3 ) from marine species such as planktonic and benthic foraminifers.Due to the temperature-dependent fractionation effect that occurs during the formation of CaCO 3 , the oxygen isotopic composition of foraminiferal CaCO 3 (δ 18 O c ) is a function of both the ambient temperature and the isotopic composition of the seawater (δ 18 O w ) in which the calcification takes place (Emiliani, 1955).Hence, δ 18 O c records form sediment cores provide information on water mass changes.
During the last few decades, stable water isotopes have been incorporated more extensively in general circulation models (GCMs).First realizations were made in atmospheric GCMs (AGCMs - Joussaume et al., 1984;Jouzel et al., 1987) and more than a decade later in oceanic GCMs (OGCMs - Schmidt, 1998;Paul et al., 1999;Wadley et al., 2002, Xu et al., 2012), which mainly focused on the linear relation between δ 18 O w and salinity, since they are affected by similar physical processes.
This topic is of significant interest in paleoceanography, because it is likely that changes in advection and freshwater budgets as well as the source of precipitation have altered the relationship (Rohling and Bigg, 1998).Using real freshwater flux boundary conditions in conjunction with the nonlinear free-surface (Huang, 1993) is essential to simulate this relationship properly.The Massachusetts Institute of Technology general circulation model (MITgcm) offers this very opportunity and further provides the adjoint method to perform data assimilation (Errico, 1997).
Here, we present first results of the implementation of stable water isotopes in the MITgcm, by performing an equilibrium pre-industrial (PI) simulation and comparing it to available observations and reconstructions.

Ocean Model
The MITgcm ("checkpoint" 64w) was employed, solving the Boussinesq, hydrostatic Navier-Stokes equations and using a nonlinear free-surface configuration (Marshall et al., 1997;Adcroft et al., 2004b).A cubed-sphere grid was used which provided a nearly uniform resolution and avoided pole singularities (Adcroft et al., 2004a).It consisted of 6 faces, each of which comprised 32 x 32 grid cells, resulting in a horizontal spatial resolution of approximately 2.8°.There were 15 vertical levels, ranging in thickness from 50 m at the surface to 690 m at the seafloor, giving a maximum model depth of 5200 m.
Associated with the non-linear free-surface is the possible vanishing of the upper layer.To avoid this problem, the rescaled vertical coordinate z* was employed (Adcroft and Campin, 2004).Furthermore, the shaved cell formulation was used, which reduced the representation error of the bathymetry (Adcroft et al., 1997).The model was coupled to a dynamicthermodynamic sea ice model with a viscous-plastic rheology (Losch et al., 2010).Isopycnal diffusion and eddy-induced mixing was parameterized with the GM/Redi scheme (Redi, 1982;Gent and McWilliams, 1990).Background vertical diffusivity for tracers was uniform at 3 • 10 -5 [m 2 s -1 ], and for the equation of state the polynomial approximation of Jackett and McDougall (1995) was used.Advection of tracers was computed using thirdorder advection with direct space-time treatment (Hundsdorfer and Trompert, 1994).
Atmospheric forcing (air temperature, specific humidity, zonal and meridional wind velocity, wind speed, (snow-) precipitation, incoming shortwave and longwave radiation as well as river runoff -12 months climatological means) was obtained from the PI ocean state estimate by Kurahashi-Nakamura et al., (submitted), which was based on the protocol of the Coordinated Ocean-ice Reference Experiments (COREs) project (Griffies et al., 2009).They optimized the forcing fields to reconstruct tracer distributions that were consistent with observations.Air-sea fluxes were internally computed in the model following the bulk forcing approach by Large and Yeager (2004).
Our simulation was initialized with present-day salinity and temperature distributions (Levitus et al., 1994 andLevitus andBoyer 1994, respectively) and spun up from the state of rest.We used asynchronous time stepping to accelerate computation with a time step of 1 day for the tracer equations and 20 min.for the momentum equations.
Furthermore, we globally balanced the freshwater flux by annually adjusting the precipitation field (Appendix A).We compiled the code using the GNU Fortran compiler gfortran version 5.3.0 and performed the simulation on 6 cores of a processor of type Intel Xeon E5-2630 v3.The simulation was integrated for 3000 years (1000 model years took ~ 7.5 CPU hours) to reach a quasi-steady state, continued for a further 3000 years with stable water isotopes as passive tracers.For analysis the average of the last 100 years was used.

Implementation of water isotopes
We implemented the stable water isotopes H 2 16 O, H 2 18 O and HDO as conservative, passive tracers in the ocean component of the MITgcm.Isotopic variations at the sea surface were driven by evaporation (E), precipitation (P) and river runoff (R), while advection, diffusion and convection affected the distribution in the interior of the ocean.Monthly climatological means of the isotopic content of precipitation and water vapor were available from the National Center for Atmospheric Research Community Atmosphere Model including a water isotope scheme (NCAR IsoCAM - Tharammal et al., 2013).Note that the prescribed atmospheric forcing fields obtained from the PI ocean state estimate by Kurahashi-Nakamura et al., (submitted) and the corresponding isotopic fluxes are not entirely consistent and might introduce an error in our model simulation.
However, to minimize the uncertainty we only took the ratio of the isotopic content of precipitation and water vapor and applied it to the corresponding atmospheric forcing fields.The isotopic composition of river runoff affects the isotopic composition of ocean water (δ 18 O w and δD w ) particularly in coastal regions.Since there was no land model in the MITgcm to calculate the amount and isotopic composition of continental runoff, we assumed that it equals the isotopic composition of the local precipitation at the river mouth and again applied it to the runoff forcing field.Fractionation during evaporation, taking both equilibrium effects and kinetic effects into account, was treated explicitly in the MITgcm.The formulation for the isotopic composition of evaporation E i [mol m -2 s -1 ] is Here, q i is the specific humidity [kg kg -1 ] multiplied by the isotopic ratio derived from the NCAR IsoCAM and is the tracer specific humidity [kg kg -1 ] in thermodynamic equilibrium with the liquid at the ocean surface (Merlivat and Jouzel, 1979).While is the local sea surface humidity [kg kg -1 ] with q sat being the saturation specific humidity [kg m -3 ] and ρ air being the atmospheric density [kg m -3 ], (5) is the local sea surface mass ratio with c being the concentration [mol m -3 ] and M the molar mass [g mol -1 ] of the respective stable water isotope.The equilibrium fractionation factor α l-v between liquid water and water vapor has been found empirically as a function of temperature and was given by Majoube (1971): Due to different molecular diffusivities of the isotopes, kinetic fractionation occurs.The kinetic fractionation factor K depends on wind speed U [m s -1 ] through the roughness of the air-sea interface (Merlivat and Jouzel 1979;Jouzel et al., 1987): The kinetic fractionation factor was used to calculate the isotopic profile coefficient Γ i following: where ρ is the air density and C E is the transfer coefficient for evaporation.
Fractionation during the formation of sea ice was neglected, because it is very small compared to other fractionation processes and thus only leads to minor effects on δ 18 O w and δD w (Craig and Gordon, 1965).Due to the absence of isotopes in the sea ice we approximated the isotopic surface flux F i [mol m -3 s -1 ] by: with A ice being the ice covered area fraction.Based on this approximation, there was no isotopic surface flux in areas covered by sea ice unless they were influenced by river runoff.Within the MITgcm, processes that affected the stable water isotopes were taken care of by the "gchem" and "ptracers" packages.While the "gchem" package added F i to the passive tracer surface tendency gPtr [mol m -3 s -1 ] the "ptracers" package mainly accounted for advection and diffusion of the isotopes.Furthermore, due to the freshwater flux that effectively changed the water column height, an additional tracer flux    [mol m -2 s -1 ] associated with this input/output of freshwater F w [kg m -2 s -1 ] was calculated following with x being an units conversion factor.   was then additionally added to the tracer surface tendency within the "ptracers" package with z [m] being the surface grid cell thickness.
In the MITgcm, the stable water isotopes were not treated as ratios, but as individual concentrations.Therefore, we initialized the ocean with homogenous concentrations of H 2 16O, H 2 18 O and HDO matching present-day δ 18 O w and δD w values of 0‰ with reference to the VSMOW.The ratios were calculated during the analysis of the results.
Furthermore, similar to the freshwater flux, a correction factor for the tracer specific precipitation was applied, whereby the respective global tracer concentration in the ocean was conserved (cf.Appendix A).

δ 18 O w data
The Goddard Institute for Space Studies (GISS) Global Seawater Oxygen-18 Database v1.21 comprises over 26,000 seawater δ 18 O values collected since about 1950 (Schmidt et al., 1999) and therefore offers an opportunity to validate the modeled oceanic δ 18 O values.
For comparison, we interpolated the GISS samples to the nearest tracer grid point of our model grid using inverse distance weighting.We excluded any data point with applied correction, from enclosed lagoons, representing estuarine or river data from near the coast or heavily influenced by meltwater, which means that we rejected all data points flagged as G. H, I, J, L and X in the GISS database (see Schmidt et al., 1999 for details -23,232 data points remained).We could not expect our model to reproduce such conditions, based on our relatively coarse grid resolution.
Since the GISS data usually represent samples taken at a certain time during the year, we did not compare them to simulated annual mean isotope values.Instead, we used a long-term mean monthly value of the specific month, when the GISS sample was measured.to determine this dependency.Since water samples are reported relative to the VSMOW standard and carbonate samples relative to the Vienna PeeDee Belemnite (VPDB) standard, the δ 18 O w values need to be converted by subtracting -0.27 ‰ (Hut, 1987).

Temperature and salinity distribution at the sea surface
We compare the simulated SST and sea surface salinity (SSS) to the World Ocean Atlas 2013 (WOA13 - Locarnini et al., (2013)), respectively (Fig. 1a, b and 2a, b).In most regions of the world ocean, SST differences are around 1 °C or even less (root mean square error (RMSE) = 1.18 °C) and therefore in good agreement with the data.Larger differences are mainly located in regions of coastal and equatorial upwelling, in the Gulf Stream and around Indonesia.
A different picture emerges for the SSS anomaly.While most parts of the surface ocean are slightly too fresh, especially the Mediterranean Sea, Bay of Bengal, Hudson Bay and north of Iceland, both the Arctic Ocean and the east coast of North America are too salty.Nevertheless, we obtain a RMSE of 0.45 psu without using any salinity restoring.

Stable water isotope distribution in ocean water
In the following, we will only focus on the comparison for δ 18 O, since the number of measurements for δD is rather small.
The surface distribution of δ 18 O w simulated by the MITgcm gradually decreases from the mid-latitudes to high latitudes (Fig. To take a closer look at the model-data fit, we interpolated the GISS data to the nearest tracer grid point (see section 2.3.1).
The separation of the model-data comparison into different ocean basins (Atlantic Ocean -Fig.5a, Pacific Ocean -Fig.5b, Arctic Ocean -Fig.5c and Indian Ocean -Fig.5d) points to deviations that mainly occur in higher latitudes.The correlation and RMSE is quite diverse, showing strong correlation for the Indian (r² = 0.77, RMSE = 0.19 ‰) and Pacific Ocean (r² = 0.74, RMSE = 0.32 ‰), medium correlation for the Atlantic Ocean (r² = 0.37, RMSE = 0.79 ‰) and no correlation for the Arctic Ocean (r² = 0.05, RMSE = 1.18 ‰).Overall, depleted δ 18 O w values are not very well simulated in the MITgcm, which is particularly recognizable in the Arctic, a region highly influenced by negative δ 18 O w values from river runoff (Yi et al., 2012).

Stable water isotopes and salinity relation
Similar physical processes determine the salinity and δ 18 O w distribution at the ocean surface.Thus locally a linear relationship between those two quantities can be expected.Restricting salinity range to 28 -38 psu in order to reflect open ocean conditions, the general features of the latter relationship are well captured in our model (Fig. 6).The modeled δ 18 O w -salinity relationship in the tropics agrees quite well with the observed one (Fig. 6a).Here, we find a simulated slope of 0.15 ‰ psu -1 , while the observed one is 0.22 ‰ psu -1 .A steeper slope is visible in the mid-latitudes for both the simulated and observed relationship (Fig. 6b).However, the agreement between those two slopes slightly decreases, with a simulated slope of 0.28 ‰ psu -1 and an observed slope of 0.49 ‰ psu -1 .Further, it underlines that we do not simulate salinity and δ 18 O w values as low as represented in the GISS data.

δ 18 O c distribution
The simulated δ 18 O c distribution at the surface increases from the tropical regions (~ 3‰) to high latitudes (~ 3.5 ‰), nicely reflecting the dependency on both δ 18 O w and temperature (Fig. 7).Even though most of the observational data is located in the North Atlantic, a latitudinal increase in δ 18 O c is recognizable.Thus the simulated large-scale pattern and latitudinal gradient match fairly well the measurements.Nevertheless, some model-data mismatch occurs.Simulated surface waters around Indonesia are too depleted as well as along coastal upwelling regions and in the Mediterranean Sea.Furthermore, regions in the North Atlantic, which are influenced by the Gulf Stream, are not sufficiently enriched compared to the observations.Overall, the simulated δ 18 O c distribution at the surface tends to be slightly too low.This becomes even clearer when performing a model-data comparison (Fig. 8a) by interpolating the δ 18 O c data to the nearest tracer grid point (analogously to the GISS data).Most of the data points plot below the 1:1 line.The cause for this bias is a matter of investigation and will be discussed in section 4.3.Nevertheless, we get a strong correlation with r² = 0.90 and RMSE = 1.01 ‰.

Sources of error
Results of the δ 18 O w distribution at the sea surface showed relatively large mismatches in the Arctic Ocean.As indicated by Eq. 11, there is no isotopic surface flux in areas that are covered by sea ice unless they are influenced by river runoff.Since parts of the Arctic Ocean are covered by sea ice all year round and others are seasonally influenced (not shown here), these areas do not experience any isotopic surface flux during most of the year.In this way the impact of precipitation and snow fall that is highly depleted is neglected, which could explain too high δ 18 O w values in the Arctic Ocean.Insufficient river discharge and neglecting the fractionation during sea ice formation could be further sources for the model-data deviations.
As part of the Pan-Arctic River Transport of Nutrients, Organic Matter and Suspended Sediments (PARTNERS) project, river discharges the most enriched freshwater (-14.9 ‰) while the water of the Kolyma river is most depleted in heavy isotopes (-22.2 ‰).This west-to-east trend is also recognizable in our model (Fig. 9), where the Ob' river contributes In the MITgcm these signals are slightly more enriched with δ 18 O w values of around 19.5 ‰ and 18.9 ‰ for the Yukon and Mackenzie rivers respectively.A consideration of the overall river discharge to the Arctic Ocean reveals a slight underestimate as the flow-weighted average for all six rivers is -18.8 ‰, while in the model it is only -17.9 ‰.Therefore, insufficient river runoff could add to the model-data mismatch in the Arctic Ocean.The general pattern of the simulated isotopic river discharge shows that river runoff is more enriched in low and mid-latitudes, which is in accord with observations (IAEA, 2012).Thus, simulating the isotopic composition of river runoff by taking the isotopic composition of the local precipitation is a reasonable first approximation, but should be overcome by implementing a bucket model in the MITgcm, which calculates the river discharge and its isotopic content for individual catchment areas over land.
During the formation of sea ice, the heavier isotopes are entrapped in the solid ice structure, while depleted sea ice brine is expelled (O'Neil, 1968).However, this fractionation process is relatively small.Lehmann and Siegenthaler (1991) reported an equilibrium fractionation constant of 2.91 • 10 -3 between pure water and ice under equilibrium laboratory conditions, while Melling and Moore (1995) estimated a fractionation constant of 2.09 • 10 -3 for ~1 m thick ice in the Canadian Beaufort Sea.So even though sea ice is highly dynamic, since it forms in one location and melts somewhere else, excluding not only the fractionation during the formation of sea ice but also the transportation of isotopes within the sea ice might lead to minor local changes but should not be one of the main sources of error.Indeed, Brennan et al. (2013) investigated the impact of a fractionation factor for sea ice on δ 18 O w within the University of Victoria Earth System Climate Model (UVic ESCM).They conclude that local changes in δ 18 O w due to the contribution of sea ice are smaller than 0.14 ‰ and therefore rather negligible.
Furthermore, the coarse resolution of our model may trigger some of the model-data discrepancies, since it is not able to resolve all of the physical processes.For instance, water that is transported towards the Nordic Seas as parts of the Gulf Stream System is displaced by water from the Labrador Sea due to the coarse grid system.Additionally, our isotopic forcing was not obtained from the same source as the atmospheric forcing for the freshwater, heat and momentum flux, whereby a maximum consistency cannot be ensured and an additional uncertainty to our sources of error is added.
Despite these sources of error, the simulated pattern of δ 18 O w both at the sea surface as well as in the deep ocean agrees fairly well with other recent studies such as the study by Xu et al. (2012) with an OGCM as well as the studies by Roche and Caley (2013) and Werner et al. (2016) with fully coupled models.

Planktonic foraminiferal δ 18 O c
To address questions regarding the evolution and history of the ocean and climate, oxygen isotopic records derived from measurements of foraminiferal shells have been used extensively.Particularly, the last glacial maximum (LGM) and last deglaciation are time periods, for which the evidence comes from proxy data recorded as oxygen isotopes in CaCO 3 .Hence, before using the model for paleostudies, an evaluation of modeled and measured δ 18 O c for the PI climate is necessary.
Apart from a few data points (red circle -Fig.8a) most of the modeled δ 18 O c values fall below the 1:1 line and thus are biased towards lower values (Fig. 8a).The outliers are all located in high latitudes in the Arctic Ocean, which is, as mentioned before, a critical area in our model.We checked whether the mismatch is caused by temperature or δ 18 O w .The temperature difference between modeled values and the WOA13 for the respective grid cells is rather small (differences not bigger than 0.15 °C) and can be ruled out as the main driving factor for our enriched δ 18 O c values.However, δ 18 O w is mainly increased by more than 1 ‰ and thus certainly leads to increased δ 18 O c values.The reason for the systematic offset of our modeled values is rather a matter of choice of the applied paleotemperature equation.The paleotemperature equation from Mulitza et al. ( 2004) is derived from plankton-tow specimens.Bé (1980), Duplessy et al. (1981) and Mulitza et al. (2004) as well as others argue that some planktonic foraminifera add an additional layer of calcite during reproduction (gametogenic calcification).This additional calcite layer is secreted in deeper and cooler water masses, introducing an 18 O enrichment in the shell.Duplessy et al. (1981) ascertained a δ 18 O mean enrichment of 0.78 ‰ and 0.92 ‰ in the shells of G. ruber and G.
sacculifer from core-top sediments, respectively.Mulitza et al. ( 2004) also showed that foraminiferal shells from the sediment are increased in δ 18 O by approximately 0.5-1 ‰.Our systematic offset of 0.74 ‰ falls into the range of these observed values and thus supports the notion that oxygen isotope values from core-top and sediment data are generally much higher than their plankton-tow counterparts.Hence, we assume that the core-top and Late Holocene data contain the signal formed by gametogenic calcification, while our δ 18 O c values calculated with the paleotemperature equation of Mulitza et al.
(2004) only represent surface water conditions.Furthermore, one has to keep in mind that we compare modeled annual mean values to measurements of foraminiferal shells that do not represent an annual average but likely contain a seasonal bias, which is not considered (Žarić et al., 2005).Additionally, the isotopic composition of foraminiferal shells can also be altered after deposition due to dissolution.This is especially the case, if the initial shell is dissolved rather than the crust formed during gametogenesis (gametogenic calcite is often more resistant to dissolution (Bé et al., 1975)), further shifting the δ 18 O towards higher values.
The observed offset vanishes, when applying the often used paleotemperature equation from Shackleton (1974) to our simulated values and comparing them to the dataset of Waelbroeck et al. (2005) (Fig. 8b).The RMSE slightly improves, while the r² and the slope almost remain unchanged (Table 2).However, this is not surprising, since the equation established by Shackleton (1974) is based on the shells of the benthic foraminifera Uvigerina spp, which, according to Zahn et al. (1986), are relatively enriched in 18 O.They assume that a low pH and decreased carbonate ion concentration, due to their habitat, might cause this latter enrichment.Based on these outcomes, we suggest, when modeling δ 18 O c of planktonic foraminifera and their respective surface water conditions, it is more appropriated to use paleotemperature equations derived from plankton-tow data (e.g.Mulitza et al., 2004) or more recent culture experiments (e.g.Spero et al., 2003), particularly since processes that affect the δ 18 O c of core-top sediments (e.g.gametogenic calcification) are not simulated and the commonly used paleotemperature equations like Shackleton (1974) or Kim and O'Neil (1997) tend to overestimate SST.
We further checked the agreement between modeled δ 18 O c and a data set derived from plankton tows (Duplessy et al., 1981, Kahn andWilliams (1981), Ganssen 1983, Bauch et al., 1997and Mulitza et al., 2003).Indeed, the MITgcm nicely reflects the δ 18 O c of living specimens (Fig. 8c) with a RMSE of 0.49 ‰ (Table 2).This further supports our assumption that δ 18 O c values from core-top sediments are a mixture of shells calcified at different water depths rather than just a pure surface water signal.
We conclude that the modeled δ 18 O c values can be compared to data successfully, although data from core-top data are generally much higher.Nevertheless, the latitudinal gradients of surface waters are well preserved and thus simulating temperature and δ 18 O w simultaneously in an ocean model provides a unique opportunity for questions regarding the climate changes.

Conclusions
Stable water isotopes have been successfully implemented in the MITgcm, using real freshwater and isotopic flux boundary conditions in conjunction with the non-linear free surface.The model captures well the broad pattern and magnitude of δ 18 O in annual mean seawater, reflecting accurately regions of net evaporation.The most enriched surface water occurs in the subtropical gyre of the Atlantic Ocean, while the surface water in the Arctic Ocean is isotopically most depleted.However, the latter ocean basin is the one with largest model-data discrepancies.They mostly result from a combination of sources of error such as the absence of highly depleted precipitation and snow fall in areas covered by sea ice as well as insufficient river runoff, which is simulated by using the isotopic composition of the local precipitation at the river mouth, but are also attributable to the coarse resolution of the model.The simulated δ 18 O w -salinity relationship is in good agreement with observations in tropical regions but less so in mid-latitudes, due to misrepresentation of δ 18 O w .But even though the δ 18 O w distribution at the sea surface reveals some deviations, the water mass structure of the deeper parts of the ocean and their characteristic δ 18 O w values are remarkably well captured in our model and show that δ 18 O w indeed can be used to characterize different water masses.Further, we tested simulated δ 18 O c against measurements of planktonic foraminiferal shells.Again, the latitudinal gradients and large-scale patterns are faithfully reproduced.However, we observed a systematic offset towards lower values in our modeled δ 18 O c .Gametogenic calcification probably causes this offset, since it is recorded in core-top and Late Holocene data but does not apply to living foraminiferal shells and their respective sea surface conditions.Hence, we suggest it is more appropriated to use paleotemperature equations derived from plankton-tow data or more recent culture experiments to investigate the δ 18 O c of planktonic foraminifera.Furthermore, for a better understanding about the factors that Dev. Discuss., doi:10.5194/gmd-2017-7,2017   Manuscript under review for journal Geosci.Model Dev. Discussion started: 21 February 2017 c Author(s) 2017.CC-BY 3.0 License.
the sea surface temperature [K].
Geosci.Model Dev.Discuss., doi:10.5194/gmd-2017-7,2017 Manuscript under review for journal Geosci.Model Dev. Discussion started: 21 February 2017 c Author(s) 2017.CC-BY 3.0 License.2.3.2 δ 18 O c data Waelbroeck et al. (2005) provide a global compilation of δ 18 O c that has been assembled within the Multiproxy Approach for the Reconstruction of the Glacial Ocean surface (MARGO) project.It contains over 2100 core-top and Late Holocene data from planktic foraminifera, mainly measured in Globigerinoides ruber pink and white, Globigerina bulloides, Neogloboquadrina pachyderma sinistral and dextral and Globigerinoides sacculifer and was used for comparison.The isotopic composition of calcite (δ 18 O c ) is linked to the temperature T during calcification and the δ 18 O w .We used the paleotemperature equation from Mulitza et al. (2004)  [°] = 14.32 − 4.28 • ( 18   −  18   ) + 0.07 • ( 18   −  18   ) 2(15) 3a, b).Highest values of about 1 ‰ occur in the subtropical gyre of the Atlantic Ocean, which are slightly higher than in the Pacific Ocean, reflecting the net freshwater transport by the trade winds.The Mediterranean Sea and Red Sea are regions of net evaporation and therefore contain δ 18 O w values of similar magnitude.The most depleted surface water is simulated in the high latitudes, showing values of -0.5 ‰ in the Southern Ocean and -1 ‰ in the Arctic Ocean.These depleted values result from negative δ 18 O w values in precipitation in combination with river/glacial runoff.Similar depleted values occur in surface waters around Indonesia.Geosci.Model Dev.Discuss., doi:10.5194/gmd-2017-7,2017 Manuscript under review for journal Geosci.Model Dev. Discussion started: 21 February 2017 c Author(s) 2017.CC-BY 3.0 License.The large scale pattern and latitudinal gradients of simulated δ 18 O w values match fairly well the observations.For example, the model captures the contrast between high and low latitudes and the Atlantic and Pacific Ocean.However, some notable discrepancies are recognizable when comparing the absolute range of δ 18 O w at the surface.The subtropical gyres are less enriched in the MITgcm (mean value of 0.6 ‰ as compared to 1.0 ‰, respectively), while δ 18 O w values in the Arctic Ocean are not as depleted as the observational data (mean value -0.6 ‰ as compared to -1.5 ‰, respectively).Especially near large river estuaries, the model-data mismatch is large.A clear distinction between different water masses based on the isotopic composition of sea water is recognizable in our simulation, both for the Atlantic Ocean and Pacific Ocean (Fig. 4a and 4b respectively).In our model, the North Atlantic Deep Water (NADW) in the Atlantic Ocean reaches down to approximately 3500 m depth and is rather enriched in H 2 18 O, resulting in a δ 18 O w content of around 0.11 ‰ (cf.Table 1).Most enriched δ 18 O w values (~ 0.6 ‰) occur in the upper water column of the tropics (20° -30 ° S and N).The NADW encounters the Antarctic Intermediate Water (AAIW) coming from the south at a water depth of approximately 1000 m.The latter is more depleted with a δ 18 O w value of around 0 ‰.The deepest parts of the Atlantic Ocean are characterized by negative δ 18 O w values of approximately -0.11 ‰ from the Antarctic Bottom Water (AABW) mixed with NADW.This water mass structure is also recognizable in the observational data and therefore in good agreement.However, the NADW is not enriched enough compared to the observational data (0.21 ‰), whereby the deepest parts of the Atlantic Ocean reveal too depleted δ 18 O w values.For the Pacific Ocean (Fig. 4b) the vertical structure is even more homogenously stratified.Enriched waters (~0.1 ‰) occur in the upper water column down to approximately 1000 m.Deeper parts of the Pacific are filled with depleted water of around -0.1 ‰.Compared to the observational data, the vertical and latitudinal gradients are in agreement.The large number of negative δ 18 O w measurements at 50° N is obtained from the Okhotsk Sea and thus is not representative for the North Pacific.
Cooper et al. (2008) published flow-weighted annual mean δ 18 O w data (collected between 2003 and 2006) from the six largest Arctic rivers.According to their estimates, δ 18 O w values of Eurasian rivers decrease from west-to-east, thus the Ob' Geosci.Model Dev.Discuss., doi:10.5194/gmd-2017-7,2017 Manuscript under review for journal Geosci.Model Dev. Discussion started: 21 February 2017 c Author(s) 2017.CC-BY 3.0 License.freshwater with a δ 18 O w value of around -15.6 ‰ and the Kolyma river of around -20.5 ‰.Even though it seems that the isotopic composition of the Ob' river is overestimated, all the other three rivers are not as depleted as seen in the PARTNERS data.Measurements of the Yukon and Mackenzie rivers reveal intermediate isotopic signals (-20.2 ‰ and -19.1 ‰ respectively).

Figure 1 :
Figure 1: Sea surface temperature anomaly (simulated in the MITgcm -WOA13) for (a) the global ocean and (b) the Arctic Ocean.

Figure 2 :
Figure 2: Sea surface salinity anomaly (simulated in the MITgcm -WOA13) for (a) the global ocean and (b) the Arctic Ocean.

Figure 3 :
Figure 3: Global surface δ 18 O w distribution simulated by the MITgcm in comparison to the observational GISS data (colored symbols (Schmidt et al., 1999) -averaged over the upper 50 m) for (a) the global ocean and (b) the Arctic Ocean.

Figure 4 :
Figure 4: Zonally-averaged cross section for the simulated δ 18 O w distribution in (a) the Atlantic and (b) the Pacific Ocean in comparison to the observational GISS data (colored symbols -Schmidt et al., 1999).Note that the GISS data does not represent a zonal mean, but rather values from specific locations.

Figure 5 :
Figure 5: Relationship between observed δ 18 O w from the GISS database (Schmidt et al., 1999) and simulated δ 18 O w from the MITgcm for the different ocean basins: (a) Atlantic Ocean, (b) Pacific Ocean, (c) Arctic Ocean and (d) Indian Ocean.Dashed lines represent the 1:1 line.

Figure 8 :
Figure 8: Relationship between measured δ 18 O c from various planktonic foraminiferas (core-top and Late Holocene (LH) -Waelbroeck et al., 2005) and modeled δ 18 O c from the MITgcm using the paleotemperature equation from Mulitza et al., 2004 (a) and Shackleton, 1974 (b).The red circle in (a) marks data points that do not fit to the observed systematic offset.The relationship between measured δ 18 O c from plankton tows (Duplessy et al., 1981, Kahn and Williams (1981), Ganssen 1983, Bauch et al., 1997 5

Figure 10 :
Figure 10: Combined T-S-δ 18 O w diagrams for the (a) GISS data and (b) simulated data in the Atlantic Ocean.The temperature and salinity ranges for the different water masses in the Atlantic Ocean are defined according to Emery and Meincke (1986).

Figure A1 :
Figure A1: Time series of the precipitation correction factor throughout the model integration.