δ 18 O water isotope in the i LOVECLIM model (version 1.0) – Part 2: Evaluation of model results against observed δ 18 O in water samples

. The H 218 O stable isotope was previously intro-duced in the three coupled components of the earth system model i LOVECLIM: atmosphere, ocean and vegetation. The results of a long (5000 yr) pre-industrial equilibrium simulation are presented and evaluated against measurement of H 218 O abundance in present-day water for the atmospheric and oceanic components. For the atmosphere, it is found that the model reproduces the observed spatial distribution and relationships to climate variables with some merit, though limitations following our approach are highlighted. Indeed, we obtain the main gradients with a robust representation of the Rayleigh distillation but caveats appear in Antarctica and around the Mediterranean region due to model limitation. For the oceanic component, the agreement between the modelled and observed distribution of water δ 18 O is found to be very good. Mean ocean surface latitudinal gradients are faithfully reproduced as well as the mark of the main intermediate and deep water masses. This opens large prospects for the applications in palaeoclimatic context.


Introduction
Water isotopes can be used as important tracers of the hydrological cycle. During phase transitions of water, such as evaporation or condensation processes, an isotopic fractionation occurs (Craig and Gordon, 1965, for example). This fractionation results from small chemical and physical differences between the main isotopic form of the water molecule (H 2 16 O, H 2 18 O). The isotopic composition of precipitation in the atmosphere has been observed to correlate with surface air temperature at mid-to high latitudes (Dansgaard, 1964) and could correlate to the amount of precipitation at low latitudes (Rozanski et al., 1993;Risi et al., 2010). In the ocean, the oxygen isotopic composition of seawater is a tracer for regional freshwater balance and water mass exchange (Östlund and Hut, 1984;Jacobs et al., 1985). As the fluxes of freshwater affect the concentration of both the oxygen isotopic composition of seawater and salinity, important regional correlation between these two parameters can be observed in most of the ocean (Craig and Gordon, 1965;LeGrande and Schmidt, 2006). Because oxygen isotope signals are preserved in an important range of records (marine and continental carbonates, ice) they are widely used as palaeoclimate proxies.
Within such a context, it is important and necessary to develop tools allowing the assessment of H 2 18 O variability under different climate conditions. The pioneering implementation of water isotopes used atmospheric general circulation models (AGCMs) (Joussaume et al., 1984;Jouzel et al., 1987) and paved the way to process based understanding of H 2 18 O in climate models. Today, most IPCC class AGCMs include the possibility for water isotopes tracing (Hoffmann et al., 1998;Noone and Simmonds, 2002;Mathieu et al., 2002;Lee et al., 2007;Yoshimura et al., 2008;Risi et al., 2010;Werner et al., 2011). The subsequent development of water isotopes modules in oceanic general circulation models (OGCM) (Schmidt, 1998;Delaygue et al., 2000;Xu et al., 2012) opens the prospect for coupled simulations of present and past climates, conserving water isotopes through the hydrosphere (Schmidt et al., 2007;Zhou et al., 2008;Tindall et al., 2009). In general, general circulation models (GCMs) have been used exclusively to simulate separately water isotopes in the atmosphere and ocean components. Indeed, given the computational constraints imposed by the use of coupled climate GCMs, it is still useful and insightful to use simpler models to simulate the evolution of water isotopes millennial timescales. Hence, the development of water isotope tracing in earth system models of intermediate complexity (EMIC) has some potential to fill this gap. The limitation in the latter class of models is the availability of enough physical mechanisms in the atmosphere modelled to provide a sufficiently realistic representation of water isotopes for our current climate in the first instance. The strongest constraint in that respect is consistent advection or transport of water since the largest spatial signals are due to along-path fractionation with the distance to the source of moisture (Craig and Gordon, 1965;Dansgaard, 1964).
This was proven possible in the CLIMBER-2 model that includes a 2.5 dimensional statistical-dynamical atmosphere with full computation for moisture advection (Roche et al., 2004), albeit under the constraint of a simplistic water closure assumption. In the present study, we evaluate the results obtained with the newly developed 18 O module for the iLOVECLIM coupled climate model, a derivative of the LOVECLIM model (Goosse et al., 2010). The difficulty does not reside presently in the advection of moisture which is computed explicitly but in the reduction of the atmosphere to three vertical layers, with a single moist layer (Opsteegh et al., 1998), complicating the vertical fractionations along the precipitation path (Roche, 2013).
The equations and verification of implementation was achieved in the first part (Roche, 2013). In the present manuscript, we present a comparison of fully coupled atmosphere-ocean-vegetation modelling results under preindustrial conditions with present-day δ 18 O data for the atmospheric and oceanic component. The performance of iLOVECLIM to sufficiently track the present-day water cycle and its potential for past climate studies will be enlightened. A third study and last part is evaluating the model from the perspective of a model-data integration with late Holocene carbonate data (Caley and Roche, 2013).

Simulation set-up
In the following, we present results of a 5000 yr equilibrium run under fixed pre-industrial boundary conditions that was used already in the first part of this study to verify the implementation of the water isotopic scheme. The atmospheric pCO 2 is chosen to be 280 ppm, methane concentration is 760 ppb and nitrous oxide concentration is 270 ppb. The orbital configuration is calculated from Berger (1978) with constant year 1950. We use present-day land-sea mask, freshwater routing and interactive vegetation.
With regards to the water isotopes, the atmospheric moisture is initialised at −12 ‰ and the ocean at 0 ‰ for δ 18 O. The consistency of our integration is checked by ensuring that the water isotopes are fully conserved in our coupled system.

Annual mean results
Starting from the annual mean distribution of δ 18 O in precipitation (cf. Fig. 1), we obtain a qualitatively good agreement with the Global Network for Isotopes in Precipitation (GNIP) data (IAEA, 2006;Masson-Delmotte et al., 2008). The main gradients (depletion towards drier and colder areas) as well as the land-sea contrasts are globally faithfully reproduced. Our results also capture the lower δ 18 O value in precipitation in the equatorial regions with respect to the tropical regions (the so-called "equatorial trough", Craig and Gordon, 1965) that is predicted from the lower E/P ratio that exists in these regions. Two exceptions may be noted: first, the too low gradient towards low δ 18 O content in precipitation in northern North America toward the Arctic ocean; second, the overall δ 18 O content in precipitation for continental eastern Africa that is too low with respect to what is expected from the measurements. The latter is due to a displacement of the zone of high tropical precipitations simulated by ECBilt towards the east in comparison to climatology, hence the higher continental fractionation over the continent. Finally, the largest bias that may be noted is the high δ 18 O content of precipitation over the Mediterranean Sea and the adjacent regions. Since this is also the case for the δ 18 O of seawater (see below), we interpret this mismatch as a consequence as the low resolution of the ocean model and thus inadequate exchange of waters between the Mediterranean and the North Atlantic Ocean. Indeed, as can be seen in Goosse et al. (2010) (their Fig. 4), running CLIO at about 3 • resolution implies the need for a 6 • latitudinal extent for the Gilbraltar Strait, hence the absence of Spain. This caveat in the model together with the low resolution of the atmospheric model displace the continentality effect over Europe and could explain the late apparition of the eastern δ 18 O trends together with much enriched values on the western side of Europe whereas data indicate more depleted values (Fig. 1). Indeed, water evaporated from the Mediterranean Sea is greatly enriched and contributes to higher δ 18 O across Europe and into Asia (LeGrande and Schmidt, 2006). Nonetheless, although it appears further inland in the model, the gradient over Europe for the δ 18 O (reflecting gradual rainout from air masses with decreasing mean annual air temperature and increasing continentality) is represented (Fig. 1), indicating a robust representation of the Rayleigh distillation, as was noted in Roche (2013).

Large-scale climatic to δ 18 O relationship
To further test the applicability of the isotopic atmospheric component, we present in Figs. 2 and 4 two classical δ 18 O to climatic variable relationship: the first compares annual mean δ 18 O to annual mean temperature (the well-known Dansgaard relationship, Dansgaard, 1964) and the second the δ 18 O to annual mean precipitation relationship.
Comparing the values obtained in the model (coloured points) with those from the GNIP dataset (grey dots), we see that the Dansgaard relationship does exist in our model (as already noted in Roche, 2013) albeit with a slightly lower slope. We interpret the lower slope as an underestimation of the fractionation towards lower temperatures/latitudes, originating from the use of the very simplified approximation of a single moist layer for the atmosphere in ECBilt (Roche, 2013). Another factor is the altitudinal effect which characterises the decrease in δ 18 O with height. This production of isotopically depleted precipitation is physically related to the Rayleigh distillation that takes place when an air parcel, lifted uphill, condenses. Taking into account the low resolution of the model that flattens all high elevation topography, we can explain why the δ 18 O values modelled in Greenland (−20 to −25 ‰) are not sufficiently depleted in comparison with data (−25 to −30 ‰) (Fig. 3). We note as well that the relationship breaks at low temperature over Antarctica, as indicated by their latitude. We have investigated this aspect in the development and verification step (Roche, 2013) and found that this mismatch is probably due to a numerical issue in the advection-diffusion scheme at very low humidity content. We have so far not been able to find a satisfactory solution to deal with this issue and thus keep it as such for the time being.
We have previously computed the linear relationship between temperature and precipitation (Dansgaard, 1964) as the regression between the mean annual temperature and the mean precipitation-weighted δ 18 O. Assuming that δ 18 O in precipitation can be used as a proxy for local temperature, only temperature at times when precipitation occurs can be recorded in the δ 18 O signal. A good example of such an issue is given by the use of the Dansgaard relationship during the Last Glacial Maximum in Greenland: during this cold period the precipitation was only falling in summer in Greenland, biasing the δ 18 O-temperature relationship (Werner et al., 2000). Therefore, it appears physically more relevant to calculate the regression between precipitationweighted temperature and precipitation-weighted δ 18 O. The precipitation-weighted temperature is computed as for a given cell i, where the sum over the j is the sum over the twelve months of the year in our case. The relationship between the obtained T * and δ 18 O in precipitation is slightly more linear (cf. Fig. 3) than the classical Dansgaard relationship (Fig. 2). Indeed, a linear fit to the modelled weighted relationship obtains a R 2 value of 0.87 against 0.83 for the standard relationship. In particular, tropical regions are more aligned with the mid-latitudes ones (as shown by the quasi   identical R 2 value obtained by a second-order polynomial fit, 0.88); same goes for the high latitude regions, apart from Antarctica. For Antarctica (Masson-Delmotte et al., 2008), the spread of anomalous points is still visible, with the possibility that the few clustered points around T * = −40 • C are more or less consistent with the linear model while in general those between −30 and −60 • C are clearly biased.
As tropical regions have a fairly constant temperature over the year and are more driven by the amount of precipitation, it is usual to plot the relation between that amount and the δ 18 O of precipitation. Model results in the δ 18 O-precipitation amount space (Fig. 4) is in very good agreement with data showing no relationship at low mean annual precipitation and a positive relationship at high precipitation amount. Our model only fails to reproduce the very few regions with very high precipitation amount, above 350 cm per year. The results are nonetheless very good given the simplicity of our atmospheric model.

Results from monthly data
To further evaluate our iLOVECLIM simulation results on a global scale, it is tempting to compare directly the seasonal evolution at specific stations representative of various climate conditions: Reykjavik (northern Atlantic), Vienna (central Europe), Ankara (eastern Mediterranean) and Belem (South America). Since our results are at relatively coarse resolution and since we have some regional biases as described in the annual mean results, we cannot expect to reproduce faithfully the mean. Hence, results are presented centred around the mean for GNIP stations and for the closest model data point (Fig. 5). While we model the correct seasonal amplitude over the year at each chosen GNIP station, we do not obtain an evolution of δ 18 O similar to the data at Belem and Reykjavik, but we do in Vienna and Ankara. As noted before, the two latter regions are areas where the mean annual discrepancy is largest between the data and the model (around 10 ‰) whereas the bias is much smaller for Belem and Reykjavik. This reinforces the notion that for a large part of Europe, the annual cycle is simply shifted to heavier values in the model with respect to the data.
The observed discrepancy between model and data for climate variables (temperature and precipitation) affects our data model comparison for oxygen isotopes. For the tropical region of Belem, we expect a more important control of precipitation on the δ 18 O signal than temperature that stays relatively constant through the year. The abnormal high precipitation rate in the model from September to December (Fig. 5) explains probably the more δ 18 O depleted values that we observe for the same period in our model whereas the data are more enriched. To evaluate whether iLOVECLIM could be used to study the interannual variability of the precipitation in the tropics, we evaluate the correlation of the interannual relationship between monthly anomalies (seasonal cycle subtracted) of δ 18 O in precipitation and temperature and precipitation rate (Figs. 6 and 7). A negative and significant correlation is found in the tropics (−30 • , 30 • ) between δ 18 O and precipitation rate, mainly over the oceanic regions. On the contrary, for the relationship between δ 18 O in precipitation and temperature, a stronger and significant correlation is observed at higher latitudes whereas the correlation is insignificant at lower latitudes. Exceptions occur in the northern high latitudes where significant negative correlations are observed in Siberia and the North America regions. We conclude that oxygen isotopes mainly record interannual variability of the precipitation in the tropics. However, in our model, results are not significant for a large part of tropical continents. The validity of these observed relationships at the interannual timescale would be tested by simulating long time periods and past climates such as the Last Glacial Maximum. By comparison, δ 18 O in precipitation also mainly record interannual variability of the precipitation in the tropics in LMD-z v4 GCM (Risi et al., 2010) and previous studies Ramirez et al., 2003).

Annual mean near-surface δ 18 O
Surface mean ocean δ 18 O results obtained (Fig. 8) are in very good agreement with data from the GISS database (Schmidt et al., 1999). The latitudinal gradients are faithfully reproduced with lower δ 18 O in high latitude regions and higher δ 18 O in tropical regions. We obtain as expected from data a lower δ 18 O around the equator the latter being a less evaporative region than the tropics, as seen already in the atmospheric part. The contrast between the evaporative zones and the mid-latitudes is well represented in the model. We observe nonetheless some notable discrepancies between the modelled distribution and data. One clearly apparent mismatch is the western Indian Ocean where we simulate much lower δ 18 O in near-surface ocean than observed in reality. This bias is due to a shift of the African precipitation regions from the west to the east of the continent, leading to much less saline waters (and unrealistically depleted δ 18 O) in the western Indian Ocean. Another region where model and data do not agree is offshore California: the isotopic signal of the North Pacific depleted values does not penetrate as far south as observed. A similar pattern is observed in the GISS model (Schmidt et al., 2007); we do not have an explanation yet for such a disagreement. Figure 9 shows details of the North Atlantic and the Arctic Ocean region. We chose the region because of the strong gradients in δ 18 O occurring due to the different water masses. Our results show an excellent match between the model results and the near-surface ocean data from the GISS database. In particular, we can clearly follow the δ 18 O enriched water masses entering the Arctic ocean, mixing with the δ 18 O depleted waters there. The fronts are generally reproduced in the right place. From these results we may infer a too little influence of Arctic waters in Baffin Bay in our model in  comparison to the data and a slightly too small western extent for the North Atlantic drift in the Nordic seas, showing the potential for regional scale evaluation of our coupled climate model. The match between the modelled and observed near-surface ocean δ 18 O is rather good in this critical region where deep water masses are formed that fill the whole North Atlantic.
Since surface water δ 18 O and salinity are affected by very similar processes (evaporation, precipitation, water masses advection and mixing), it is important to ascertain whether the observed δ 18 O to salinity relationship can be reproduced in our control simulation. Figure 10 shows that for the Atlantic Ocean, simulated δ 18 O-salinity is in very good agreement with observed data for the same ocean basin for all latitudes. We restricted the lower salinity range to 33 ‰ since the GISS database also contains data points with very nonsalty conditions, reflecting rivers mouths and not open ocean. Since we are using a relatively coarse grid model, we cannot expect to reproduce such conditions. The obtained modelled gradient for the δ 18 O-salinity relationship is 0.43 (Roche, 2013), in close accordance with the observed one (0.52),  while still bearing the too low fractionation towards high latitude, as observed in the atmospheric model.
Regarding the Pacific Ocean near-surface waters (Fig. 11), while we still obtain a fair agreement, modelled data points from the northern tropical latitudes to the north are shifted towards higher δ 18 O values. This reflects the overestimation of δ 18 O already seen in Fig. 8 especially in the southern tropics. Since the salinities are in good accordance with the data in that particular region it implies that there is an overestimation of fractionation at evaporation in this region. One likely source is an underestimation of near-surface ocean humidity since this term is the most effective in the governing equation (Roche, 2013). The notable difference between the northern Pacific and the southern Pacific yield two parallel lines, the northern one being offset half a per mil in δ 18 O and one per mil in salinity. Though the spread observed is still within the observed data range, further investigation would be needed to fully understand that difference.

Annual mean deep ocean δ 18 O
Since we obtain satisfactory results at the surface of the oceans it is worth comparing our results to available data for the interior of the oceans. It shall be noted that representing the deep ocean δ 18 O is complicated since it depends strongly on the δ 18 O content at the location of deep water formation, a rather restricted area. A relatively small shift in δ 18 O may lead to a substantial drift in deep ocean δ 18 O with an opposite drift in surface waters δ 18 O. The cross sections presented are averaged over all the considered ocean basin, using a mask to define the ocean basins. The data points from the GISS database are all collapsed on the same cross section but not averaged. Thus, certain data points may represent values from specific locations that are not comparable to the mean.
In the Atlantic Ocean (Fig. 13), the simulated δ 18 O distribution clearly show the mark of the main water masses. From the north, the North Atlantic deep waters (NADW) are marked by higher δ 18 O content, around 0.2 ‰. NADW extends in our model almost to the bottom of the ocean (as shown in Goosse et al., 2010) where it mixes with the Antarctic bottom water (AABW) coming from the south. The latter water mass is marked by low δ 18 O content. In the Southern Ocean, the Antarctic intermediate waters are also marked  (Schmidt et al., 1999) for the Atlantic Ocean. Colour scale indicates the latitude for the iLOVECLIM circles.   (Schmidt et al., 1999) (coloured scale) presented as a latitudinal-depth (m) transect. field is otherwise in very good agreement with data, we infer that this discrepancy may arise from two distinct causes. First, the fact that we compare a zonal average field in the model versus point-based measurements from data tends to smooth the east-west contrasts with a western deep Atlantic having higher δ 18 O content (by 0.1 to 0.2 ‰) than the eastern deep Atlantic, as shown on Fig. 12. The intrusion of Labrador deep water is particularly prominent in that respect, with values of 0.3 ‰. Second, the fact that even when taking this aspect into account there is still an underestimation of the δ 18 O in the model (see in the tropical deep Atlantic on Fig. 12 for example) shows that the influence of NADW with positive δ 18 O content might not be strong enough in the deep  (Schmidt et al., 1999) (coloured scale) presented as a latitudinal-depth (m) transect.
North Atlantic. Interestingly, Fig. 12 clearly shows that the southern deep Atlantic is reproduced faithfully in the model. This indicates that the entrance of southern source water in the southern deep Atlantic is correctly represented and is not the cause of the problem mentioned. Whether the cause is not enough northern sourced deep water export to the deep ocean or inadequate mixing of that northern sourced waters in the ocean interior is a matter for investigation.
For the Pacific Ocean, presented in Fig. 14, the results are similar though less clearly visible. We obtain a good agreement for the δ 18 O distribution until 2000 m deep in all the Pacific and reasonable low δ 18 O content around the Antarctic at all depth. As in the Atlantic, the deep ocean (from 30 • S) has a too low δ 18 O content in our simulation by about 0.2 ‰. Again, we can only state that it may be due to incorrect mixing in the deep ocean or not enough influence of δ 18 O rich waters. Closer to the surface, our simulation reproduce faithfully the vertical and latitudinal gradients. In the North Pacific, the data is dominated by the large number of points from the Okhotsk Sea that does not represent the mean of the ocean basin. In addition, we have seen previously that models (iLOVECLIM and GISS) failed to reproduce the North Pacific depleted values in offshore California. Together, this can explain the disagreement between model and data for the North Pacific zonal mean comparison.

Conclusions
In the present manuscript, we have evaluated a pre-industrial control run against present-day observation of δ 18 O in precipitation and in the ocean waters. For the atmospheric part, we found that apart from central Antarctica, our model produces results in good accordance with what is known from the GNIP dataset, with some caveats due to the simplicity of the approach taken. If the general gradients (latitudinal, oceanic to continental) are well reproduced, the trends towards colder and drier are generally underestimated, indicating a too low fractionation in our single moist layer atmosphere.
From the ocean perspective, we obtain a very good agreement between the simulated δ 18 O values in the oceans and the observed values as recorded in the GISS database. The slightly too low gradient found in the atmosphere is naturally also present in the oceanic part of our coupled model since the system is coupled and closed with respect to water and isotopic water content.
Overall, with the caveats mentioned above, we obtain a water isotope enabled coupled climate model well suited for long-term simulation of the changing climate. Since our aim is to use the model in palaeoclimatic context, the next logical step is to evaluate the model against available δ 18 O proxy for the late Holocene. This is the subject of the third part of our study (Caley and Roche, 2013) and will enable us to identify the advantage and caveats of our model in such comparisons under the well-known conditions of the pre-industrial climate.