Biogeosciences Mismatch between observed and modeled trends in dissolved upper-ocean oxygen over the last 50 yr

Observations and model runs indicate trends in dissolved oxygen (DO) associated with current and ongoing global warming. However, a large-scale observation-tomodel comparison has been missing and is presented here. This study presents a first global compilation of DO measurements covering the last 50 yr. It shows declining upperocean DO levels in many regions, especially the tropical oceans, whereas areas with increasing trends are found in the subtropics and in some subpolar regions. For the Atlantic Ocean south of 20 ◦ N, the DO history could even be extended back to about 70 yr, showing decreasing DO in the subtropical South Atlantic. The global mean DO trend between 50 ◦ S and 50 N at 300 dbar for the period 1960 to 2010 is−0.066 μmol kg−1 yr−1. Results of a numerical biogeochemical Earth system model reveal that the magnitude of the observed change is consistent with CO 2-induced climate change. However, the pattern correlation between simulated and observed patterns of past DO change is negative, indicating that the model does not correctly reproduce the processes responsible for observed regional oxygen changes in the past 50 yr. A negative pattern correlation is also obtained for model configurations with particularly low and particularly high diapycnal mixing, for a configuration that assumes a CO 2-induced enhancement of the C : N ratios of exported organic matter and irrespective of whether climatological or realistic winds from reanalysis products are used to force the model. Depending on the model configuration the 300 dbar DO trend between 50 ◦ S and 50 N is −0.027 to−0.047 μmol kg−1 yr−1 for climatological wind forcing, with a much larger range of −0.083 to +0.027 μmol kg−1 yr−1 for different initializations of sensitivity runs with reanalysis wind forcing. Although numerical models reproduce the overall sign and, to some extent, magnitude of observed ocean deoxygenation, this degree of realism does not necessarily apply to simulated regional patterns and the representation of processes involved in their generation. Further analysis of the processes that can explain the discrepancies between observed and modeled DO trends is required to better understand the climate sensitivity of oceanic oxygen fields and predict potential DO changes in the future.


Introduction
Anthropogenic carbon dioxide emissions alter our oceans in numerous ways, some of which have reached considerable attention, like ocean acidification and sea level rise.Recently, concern has emerged that CO 2 -driven climate change causes decreasing dissolved oxygen (DO) levels in the ocean (Keeling and Garcia, 2002;Bopp et al., 2002;Plattner et al., 2002;Frölicher et al., 2009) with potentially large impacts on marine habitats and ecosystems (Keeling et al., 2010).Volumes of the interior ocean with dissolved oxygen concentrations below the lethal threshold for many marine habitants are often called oxygen minimum zones (OMZs).As marine macro-organisms show different thresholds, different boundaries are used to define OMZs, varying from suboxic levels to more relaxed concentrations of 90 µmol kg −1 (Karstensen et al., 2008).Century-long climate model runs project an overall decline in oxygen concentration (Matear and Hirst, 2003;Schmittner et al., 2008).These simulated changes Published by Copernicus Publications on behalf of the European Geosciences Union.
L. Stramma et al.: Observed and modeled oceanic oxygen trends have, to a large extent, been explained by changes in ocean circulation and decreasing gas solubility under global warming conditions (Bopp et al., 2002).Declining oxygen levels may further result from direct impacts of elevated CO 2 concentrations on the marine biology.These include potential acidification-induced decline in calcium-carbonate ballast and associated change in export and remineralization profiles (Hofmann and Schellnhuber, 2009), or enhanced C : N ratios of postulated CO 2 -fertilized organic matter production (Oschlies et al., 2008).Particularly in the low-oxygen regions of the tropical oceans a local decline in DO can have global implications: whenever thresholds to suboxic conditions are surpassed, remineralization of organic matter continues by consuming fixed nitrogen instead of oxygen.The extent of suboxic conditions thereby controls the loss of fixed nitrogen from the ocean with direct impacts on global ocean productivity (Keeling et al., 2010).
OMZs exist in the eastern basins of the Pacific and Atlantic Ocean and in the northern Indian Ocean.These OMZs cover the upper ocean from about 100 to 900 m depth with a total lateral extent of the permanent OMZ's defined as DO < 20 µM (= µmol l −1 ) of about 8 % of the total oceanic surface area (Paulmier and Ruiz-Pino, 2009).Continuous time series in a few selected tropical areas with sufficient data (Stramma et al., 2008) as well as comparison between two time periods for the tropical oceans (Stramma et al., 2010a) indicate an ongoing oxygen decline and a vertical expansion of the OMZ in most tropical regions, contrasted by areas with DO increase, predominantly in the subtropical gyres.The probably longest continuous DO time series exists at the Ocean Station Papa (50 • N, 145 • W) in the subpolar North Pacific.It shows both multi-annual to bi-decadal oscillations as well as persistent trends of DO decline (Whitney et al., 2007).The DO trend might not be monotonous but will be influenced by variability on different time-scales e.g. in the Pacific by El Nino cycles or the Pacific Decadal Oscillation (PDO).
Along the eastern ocean margins, in the vicinity of the OMZs, productive upwelling systems support a large proportion of the world's fisheries.DO time series in the eastern North Pacific, e.g. in the California Current (Bograd et al., 2008) and off Oregon (Chan et al., 2008), show shoaling OMZs.The expansion of OMZs into shallower waters may interact with natural or eutrophication-induced hypoxic zones on the shelf and will intensify the processes increasing the count of dead zones (Diaz and Rosenberg, 2008).In the open ocean the shoaling of the tropical OMZ restricts the depth distribution of tropical pelagic fishes such as marlins, sailfish and tuna by compressing their habitat into a narrow surface layer (Prince and Goodyear, 2006) with a clear correlation with open ocean DO content (Prince et al., 2010;Stramma et al., 2012).
Interestingly, models have, so far and in contrast to the observations over the past decades, mostly shown an increase in oxygen concentrations in the tropical thermocline (Bopp et al., 2002;Matear and Hirst, 2003).These models show an associated reduction in the extent of suboxic areas under projected 21st century global warming, unless CO 2 -induced changes in C : N ratios (Oschlies et al., 2008) or CO 2 -induced ballast effects were assumed (Hofmann and Schellnhuber, 2009).The tendency towards higher oxygen concentrations in the models' tropical thermoclines under global warming has been attributed to enhanced stratification and reduced mixing of old and low-oxygen waters from below at the expense of a more immediate contact between the shallower thermocline and the sea surface (Gnanadesikan et al., 2007).Keeling et al. (2010) suggest that this increase in tropical thermocline DO levels might be an artifact resulting from an excessive diapycnal mixing in the thermocline of coarseresolution models.A recent modeling study, in contrast, finds that models with unrealistically high diapycnal mixing tend to simulate expanding suboxic water volumes, whereas models with low mixing intensity tend to simulate shrinking suboxic volumes over the 21st century (Duteil and Oschlies, 2011).With respect to direct observations of diapycnal mixing intensities in the subtropical thermocline (Ledwell et al., 2003), inferences about a reliable representation of such mixing levels in current ocean circulation models (Getzlaff et al., 2012), and metrics using nutrient and radio-carbon distributions (Duteil and Oschlies, 2011), models with lower mixing intensity appear more realistic than those with high mixing intensity.
The aim of the present work is a confrontation of biogeochemical model results with observed large-scale DO changes for which a new global compilation of oxygen measurements taken over the past decades is constructed.This allows, for the first time, a global overview of the DO changes over the last 50 yr, which is not confined to the difference of two temporal periods, and covers more than just the tropical and subtropical regions.In addition, observed differences covering about 70 yr in the tropical and South Atlantic are presented.A suite of model runs is employed to examine to what extent observed multi-decadal trends can be reproduced by current biogeochemistry-climate models.In particular, we investigate whether specific assumptions about the strength of diapycnal mixing, about direct CO 2 effects on the stoichiometry of the biological pump, and about historical wind patterns yield particularly good or bad representations of the observed changes in dissolved oxygen.

Observational data
The basic data set used for our analysis is the freely available HydroBase-2 data set (Curry, 2008).The quality-controlled HydroBase-2 data available as of 10 October 2008 are augmented with additional recently collected data sets (Stramma et al., 2010a).Oxygen data available for the time period 1960 to 2010 is used.Spatial data coverage changes over time, following the changing foci of research (Fig. 1).For the earlier years, the data distribution shows a focus on densely sampled station data near the continents with few trans-ocean sections.In contrast, DO data coverage is dominated by the trans-oceanic sampling pattern of the World Ocean Circulation Experiment (WOCE) since the 1990s and generally sparse since 2000.This leads to a temporal trend covering slightly different periods in different location.An additional oxygen data set is extracted from the historical data for the tropical and South Atlantic for the measurements of the Meteor 1925-1927expedition (Wattenberg, 1938).
There are not enough oxygen measurements available to derive past variability and trends in great detail in the ocean.More or less continuous oxygen time series exist only at a few locations (e.g.Whitney et al., 2007).To obtain trends for the open ocean with limited data availability, one can bin data from larger areas to construct a single "average" time series valid for this area (e.g.Stramma et al., 2008).Alternatively, to gather a reasonable amount of data, time periods of the order of decades can be compared to other time periods to study changes that occurred between the two periods (Stramma et al., 2010a).While quality control and data preparation are similar to that of Stramma et al. (2010a), the approach taken here combines the advantages of both procedures: It does not restrict the analysis to two periods, but takes into account anomalies on all time scales in a global mapping of multi-decadal changes.The drawback is that the trend is computed over the period of locally available data and not over a uniform period of time.On the other hand biases introduced in trend analysis using other methods, like for instance the difference between two periods, are likely larger since that approach does not make use of all available data, and spatial extrapolation (necessary to create the distribution for each period) has to be applied over larger distances than in the procedure chosen here.Despite the fact that our method is less likely to overestimate the significance of trends, decadal variability certainly can bias our maps depending on temporal data distribution.
We map all available DO data since 1960 onto a 1 • × 1 • grid between 70 • N and 70 • S. All DO profiles with measurements between 275 dbar and 325 dbar are used.Data points with implausible physical parameters are rejected.Data within this depth interval spanning 50 dbar are interpolated onto the 300 dbar pressure surface.The local DO anomaly on this surface is then derived by the difference of the measurement to the DO climatology of the World Ocean Atlas (WOA) 2009 (Garcia et al., 2006).The data set underlying the interpolated WOA DO climatology is essentially identical to the HydroBase-2 data set.Anomaly data are further quality controlled by an interquartile range (IQR) filter, rejecting data exceeding three times the IQR above the third or below the first quartile.A first quartile is larger than 25 % of the data, the second quartile is the median (larger than 50 % of the data) and the third quartile represents 75 % of the data.The IQR represents the center 50 % of the data.An IQR filter does not require a normal Gaussian distribution and is thus a more objective way to quality control data in comparison to a standard deviation threshold (e.g.Smith and Sandwell, 1997).All anomaly data is then binned in a global 1 • × 1 • grid for each year since 1960.This minimizes distortions during the mapping process, which can arise from sparse irregular temporal/spatial sampled areas next to intensively sampled areas like oceanographic stations with an abundance of oxygen data.Thus any measurements within one year and a 1 • × 1 • box are combined to one data point by taking the mean of each parameter.The DO anomaly data is than mapped with a least squares linear model approach, using locally weighted scatter plot smoothed with 1500 km standard tri-cube distance weighting (see Cleveland, 1979).The local trend and standard error is computed over the period of locally available data.The local trend derived is a weighted least square solution to a linear system on time and observed DO values.The solution is found by minimizing the system with distance weights applied.The mean squared error and the estimated standard errors can be derived in this process, as detailed in the literature (e.g.Strang, 1986).The validity of the method was tested by subsampling results of our numerical model (see below) at the times and locations of the real observations, and then comparing the DO trend by the interpolation method from the subsampled "synthetic" data set against the "true" model trend computed from the complete model output (see Fig. 4

below).
Between the 1925-1927 Meteor Expedition and the 1960s only few DO data with sometimes questionable accuracies exist.Hence, instead of estimating the trend since 1925, the DO differences between the 1925-1927 expedition and a mean oxygen distribution for the period 1990 to 2008 derived as in Stramma et al. (2010a) was computed.These differences cover a time period of about 70 yr.

The model
The model used is the University of Victoria (UVic) Earth System Climate Model (Weaver et al., 2001;version 2.8) in the configuration described and validated against observed global distributions of nutrients, oxygen, dissolved inorganic carbon and alkalinity by Schmittner et al. (2008) and Oschlies et al. (2008).The oceanic component is a fully threedimensional primitive-equation model with 19 levels in the vertical ranging from 50 m thickness near the surface to 500 m in the deep ocean.It contains a simple marine ecosystem model with the two major nutrients nitrate and phosphate and two phytoplankton classes, nitrogen fixers and other phytoplankton, with the former being limited only by phosphate.Detritus sinks with a sinking velocity increasing linearly with depth from 7 m day −1 at the surface to 40 m day −1 at 1000 m depth and below.This ocean component is coupled to a single-level energy-moisture balance model of the atmosphere, a dynamic-thermodynamic sea-ice component, and × 3.6 • .Diapycnal mixing is parameterized as the sum of tidally induced mixing (Simmons et al., 2004) and background mixing described by a uniform diffusivity k bg and an additional Southern Ocean enhancement of the vertical diffusivity by 1.0 cm 2 s −1 .The tidally induced diffusivity rapidly decays in the water column above the seafloor with an exponential vertical scale of 500 m, so that in most parts of the thermocline diapycnal mixing is determined by the background diffusivity.
This mixing background is set to 0.15 cm 2 s −1 in the standard simulation (Oschlies et al., 2008), to 0.01 cm 2 s −1 in a low-mixing sensitivity experiment, and to 0.3 cm 2 s −1 in a high-mixing sensitivity experiment.Duteil and Oschlies (2011) show that all three model configurations yield a relatively good agreement with the observed large-scale distribution of pre-bomb δ 14 C. Another sensitivity experiment uses the standard model's 0.15 cm 2 s −1 background diffusivity, but instead of the constant Redfield C : N stoichiometry of or-ganic matter, assumes that C : N ratios increase with increasing atmospheric CO 2 levels according to a linear response inferred from mesocosm experiments run under different atmospheric CO 2 levels (Riebesell et al., 2007;Oschlies et al., 2008).
All model configurations are spun up for more than 10 000 yr under pre-industrial atmospheric and astronomical boundary conditions and subsequently run under historical conditions from year 1765 to 2000 using fossil-fuel and land-use carbon emissions.Following the IPCC AR4 protocol, we switch to emissions according to the SRES A2 (Nakicenovic et al., 2000) scenario in year 2000 and continue the model to year 2010.Emission data now available for the first years of the 21st century are very close to the SRES A2 predictions (Manning et al., 2010) and we do not expect our model results to be sensitive to small emission changes in a few years.The model results shown here correspond to the simulations of Oschlies et al. (2008) and Duteil and Oschlies (2011).In contrast to the earlier studies that concentrated on potential future changes in global oxygen budgets, the current paper investigates the spatial patterns of the simulated changes and confronts these with observational estimates.An important feature of this Earth system model of intermediate complexity is the considerably simplified atmospheric model component, which requires prescribed wind fields.The standard configuration of the model employs seasonally cycling climatological wind fields taken from the NCAR/NCEP monthly climatology (Kalnay et al., 1996).In this configuration, the model cannot simulate natural or anthropogenically induced interannual and longer-term variability in the wind field and instead focuses on CO 2 -driven changes in the radiative energy balance.The role of interannually varying wind forcing is explored in two sensitivity simulations that employ realistic monthly-mean wind forcing of the Co-ordinated Ocean-ice Reference Experiments (CORE, Large and Yeager, 2004) for the period 1958-2007.

Observed oxygen changes
The investigation of DO changes over the past 50 yr was carried out for the 300 dbar layer, which is at about 300 m depth.A pressure layer is often used, since a pressure layer is easier to conceive for a reader than an isopycnal layer varying in depth in space and also in time.The 300 dbar pressure surface covers the upper part of the OMZ of the tropical oceans.
Adding the mapped DO anomalies to the coarse WOA fields gives a higher resolution DO distribution than using WOA DO fields alone.Figure 2a shows that despite the large horizontal influence radius employed in the mapping routine, the derived mean dissolved oxygen distribution fits very well the known distributions of the OMZs in the tropical eastern Pacific, tropical eastern Atlantic and northern Indian Ocean.Similar distributions are described in earlier large-scale surveys (e.g.Karstensen et al., 2008;Paulmier and Ruiz-Pino, 2009).Well visible is the poleward spreading of the OMZ in the eastern basins along the continental slopes, especially in the Pacific as well as in the North Atlantic.
The DO changes on the isobaric 300 dbar surface for the time period 1960 to 2010 show some large variability and trends of opposite signs in all oceans, with stronger trends often above the 95 % confidence standard error interval (Fig. 2b).The tropical open oceans between 20 • N and 20 • S are dominated by negative trends of up to 0.83 µmol kg −1 yr −1 .The subtropical and subpolar oceans show a more equal distribution between positive as well as negative trends.Positive trends reach up to 0.68 µmol kg −1 yr −1 in the extratropics.In particular, some regions of the subtropical gyres indicate a DO increase.
The circulation of the subtropical gyres has intensified in recent years and is likely to have a variability on decadal time scales, e.g. an intensification of the circulation of the South Pacific Ocean subtropical gyre, extending from the sea surface to mid-depth, was observed in the 1990s (Roemmich et al., 2007).A comparison of WOCE and Argo float trajectories at 1000 dbar confirms the gyre spin-up during the 1990s (Roemmich et al., 2007).Such an enhanced circulation could accelerate the transport of freshly ventilated oxygen-rich wa-ters of the subtropical gyres into the subtropical thermocline, consistent with the observed DO increase.
The DO decrease in the eastern North Pacific, interrupted only by a band of positive DO change at about 20 • N is in agreement with results by Deutsch et al. (2005) showing that the DO decrease in the subtropical North Pacific is due to long-term reduction in ventilation, while the DO increase along the southern flank of the subtropical gyre is caused by a shift of the gyre boundary.
The subpolar regions in the northern hemisphere show mainly decreasing oxygen trends (Fig. 2b).DO changes in the subpolar and subtropical North Atlantic south of Iceland can be related in part to the North Atlantic Oscillation as was reported for the period 1988 to 2003 (Johnson and Gruber, 2007) although model runs assessed that the NAO accounts for only 30 % of the total variability in the subpolar and subtropical North Atlantic (Frölicher et al., 2009).In the southern hemisphere positive as well as negative oxygen trends appear (Fig. 2b).
Large areas of DO increase are found in the subtropical gyre of the South Indian Ocean.A similar DO increase was observed over most parts of hydrographic sections in the South Indian Ocean at about 32 • S from 1987 and 2002 (Mc-Donagh et al., 2005).A survey of changes in the Southern Ocean (Aoki et al., 2005) for three regions in the western Atlantic and south of Australia between 1960-1980and 1980-1996 showed mainly an oxygen decrease, however, an oxygen increase was reported for a part of the Atlantic area.
The mapped observation-based trends can be compared to local time series data.The best record of a long term trend is the DO time series at 50 • N, 145 • W in the eastern North Pacific showing a DO decrease of 0.39 to 0.70 µmol kg −1 yr −1 at 100 to 400 m depth (Whitney et al., 2007).The negative DO trend derived by our method of 0.30 µmol kg −1 yr −1 at 300 dbar agrees reasonably well with these local observations (Fig. 2b), the differences might be related to the interpolation used here and the restriction to just one pressure layer.For several areas in all three tropical oceans Stramma et al. (2008) identified negative oxygen trends for the 300-700 m layer over the last 50 yr.These areas coherently appear as regions of DO decrease in our new analysis as well.Summaries of observed local changes (Keeling and Garcia, 2002;Keeling et al., 2010) list decreasing oxygen levels for most observations similar to the large-scale global patterns derived here.
The mean trend for the 300 dbar layer between 50 • N and 50 • S is −0.066 ± 0.051 µmol kg −1 yr −1 hence a mean decrease in oxygen.Helm et al. (2011) reported a global oxygen decrease of 0.93 µmol l −1 for the comparison of a period centered at 1970 to the one centered at 1992 for 5 • × 10 • grid boxes between 100 and 1000 m, which leads to an annual trend of −0.042 µmol l −1 and which is of similar order as our trend for 300 m depth.In the tropical ocean Helm et al. ( 2011 Off California the observed oxygen anomalies were high in the mid-1980's due to the PDO, but were lower in the 1960's and the 2000's, related to anomalies in suboxic volume (Deutsch et al., 2011).Although our observations cover the entire period between the two low oxygen periods, a larger negative trend might be related to the higher PDOrelated oxygen in the 1980's and 1990's as many data were collected in this time period during WOCE.
A comparison with the DO changes at 200 dbar and 200 to 700 dbar between the periods 1960-1974 and 1990-2008 derived between 40 • N and 40 • S by Stramma et al. (2010a) shows similarities but also some differences with respect to our results obtained here for 300 dbar, e.g. the positive trend in the eastern tropical Pacific off Peru at 300 dbar.The differences can partly be explained by the different depth layers used and the different time periods investigated.The changes related to different time periods are illustrated here by exploiting one of the earliest basin-scale data set of oxygen measurements available, namely from the 1925-1927 Meteor Expedition.
Despite difficulties at earlier times to correctly determine the geographical position as well as the depth of the samples from bottles attached to a wire, we present the differences of the DO measurements from the Meteor Expedition to the period 1990 to 2008.Identical quality control thresholds were applied to this data as to the DO analysis above.The differences at 300 m depth show the largest DO changes in the subtropical South Atlantic (Fig. 3) even in the area off Argentina where the DO trend for the last 50 yr was positive.In the tropical Atlantic there is not such a clear DO decrease from 1925-1927 to present as for the trend since 1960 to present.This indicates lower oxygen concentrations in some regions of the tropical Atlantic in the late 1920s than in the 1960s.The differences in the trends since 1925-1927 compared to the trends since 1960 suggest that multi-decadal differences in observed oxygen concentrations may well be dominated by temporal variations other than linear trends.This result might be caused in part by the different methods used and the reduced data accuracy for the 1925-1927 measurements.

Modelled oxygen changes
We first investigate changes in dissolved oxygen concentrations simulated by the UVic Earth System model configuration that is forced by anthropogenic CO 2 emissions and seasonally cycling climatological wind fields and thereby effectively filters out wind-driven interannual and longer-term variability.Results of this model configuration show a decline of DO concentrations over large parts of the ocean over the period 1960 to 2010 at 300 dbar (Fig. 4a).Typical rates of the simulated DO decline (300 dbar average is −0.027 µmol kg −1 yr −1 in the standard configuration, see Table 1) appear consistent with the observed trend and, in agreement with the observations, the model simulates a substantial DO decline in the northwestern North Pacific, in the North Atlantic and in parts of the Southern Ocean.However, general patterns of simulated DO changes differ considerably from those observed.In fact, the global pattern correlation between the 300 dbar oxygen trend estimated from observations (Fig. 2b) and simulated by the standard model (Fig. 4a) is negative (−0.21 when considering only those areas on the 300 dbar surface with significance levels of the observational estimate exceeding 95 %, see Fig. 2b).The difference between simulated and observed DO changes is most pronounced in the tropical oceans where simulated DO shows a slight increase from 1960 to 2010 in the standard model configuration while the observations show the largest DO decrease in the zonally averaged view (Fig. 5).A global-warming induced increase in the tropical thermocline DO was also found in the earlier model studies by Bopp et al. (2002) and Matear and Hirst (2003).
As pointed out by Gnanadesikan et al. (2007), a local oxygen increase in the tropical oceans can, in the models, be related to enhanced stratification and a shoaling of the OMZ, resulting in enhanced supply of younger and oxygen-rich waters from above into the shallower OMZ.This argument relies on the presence of a substantial diapycnal tracer flux, which may be systematically overestimated in typical coarseresolution climate models (Keeling et al., 2010).To investigate the impact of the strength of diapycnal mixing on the simulated evolution of tropical oxygen levels, we performed a sensitivity experiment with the vertical background diffusivity reduced from 0.15 cm 2 s −1 to 0.01 cm 2 s −1 , and a second sensitivity experiment with an elevated vertical background diffusivity of 0.3 cm 2 s −1 .Although the results of the low-diffusivity experiment show some reduction in the areal extent of increasing tropical DO levels (Fig. 4c) and a somewhat larger average DO decline (−0.034 µmol kg −1 yr −1 , the overall correlation with the observed pattern of DO changes is still negative (−0.16 for the observational > 95 confidence area).A reduction in the extent of tropical regions with a simulated increase in DO is also found for the high-diffusivity sensitivity experiment (Fig. 4d), in which DO declines on average by −0.041 µmol kg −1 yr −1 , again with a slightly negative pattern correlation (−0.09) with observed DO trends.
Employing even higher vertical diffusivities of 0.5 cm 2 s −1 yields less and less agreement with observed δ 14 C distributions (Duteil and Oschlies, 2011), while displaying a mean DO decline of −0.037 µmol kg −1 yr −1 with a slightly positive pattern correlation (+0.01) with the observational estimate of DO change.Pattern correlations are less negative for the low and high mixing sensitivity runs mainly because of the reduced positive simulated DO trends in the tropics (Fig. 4).The reduced positive trend at very low diapycnal mixing rates can, in our model, be explained by a relatively large poleward movement of isopycnal outcrop areas under global warming and at low vertical mixing, and the associated increase in the distance waters have to travel before reaching the tropical thermocline.In the high-mixing case, the reduced positive trend can be related to the enhanced downward transport of warming surface waters with lower oxygen solubility into the thermocline.
The sensitivity experiment with C : N ratios increasing proportionally to atmospheric CO 2 levels yields results which are very similar to those of our constant C : N standard configuration shown Fig. 4a.Although the two model configurations were found to simulate different signs of DO changes in the tropical thermocline over the 21st century under a business-as-usual emission scenario (Oschlies et al., 2008), assumed CO 2 -dependent variations in C : N ratios are yet too small over the period 1960 to 2010 to have a substantial impact on simulated DO fields.
Thus, in the zonal mean all model configurations with climatological wind forcing show an increase in tropical DO at 300 dbar, whereas the zonally averaged trend is most negative in the tropics (Fig. 5).Zonally averaged observed and simulated DO trends agree better at mid and high latitudes, while the simulations cannot reproduce the observed zonally averaged oxygen increase in the subtropics.Given the overall poor agreement of both standard and low-mixing simulations with observed patterns of oxygen change in the tropical ocean, we postulate that processes other than excessive diapycnal diffusion in the model, as had been suggested by Keeling et al. (2010), must be responsible for the systematic discrepancy.
To better understand why the models simulate increasing DO in the tropical thermocline, we isolated the impact of the individual processes responsible for changes in dissolved oxygen concentrations.Changes in solubility are mostly driven by temperature and lead to a decrease in saturation concentrations when temperatures increase.At 300 dbar, solubility decreases almost everywhere in the model, with an average value of about 0.03 µmol kg −1 yr −1 (Fig. 6).Local changes in biotic oxygen consumption simulated by the model on the 300 dbar surface turn out to be about an order of magnitude smaller, with a mean value of about −0.002 µmol kg −1 yr −1 .Most of the regional signal in simulated DO change is attributed to changes in transport pathways and processes, which are here diagnosed as the residual between total simulated DO changes and the sum of shows the trend estimated for the model results subsampled according to the dates and locations of the read DO measurements and using the same routine that was applied to the in situ data in Fig. 2. Units are µmol kg −1 yr −1 .Fig. 5. Zonally averaged change in observed oxygen at 300 dbar (dashed black line) and simulated oxygen at 300 dbar for the standard configuration of Duteil and Oschlies (2011) with a diapycnal background diffusivity of k bg = 0.15 cm 2 s −1 and constant (Redfield) C : N stoichiometry (black curve, Fig. 4a), for a much reduced background diffusivity of k bg = 0.01 cm 2 s −1 and constant C : N stoichiometry (red curve, Fig. 4c), for increased background diffusivity of k bg = 0.3 cm 2 s −1 and constant C : N stoichiometry (green curve, Fig. 4d), and for k bg = 0.15 cm 2 s −1 and pCO 2 -sensitive C : N ratios as used by Oschlies et al. (2008) (blue curve).solubility and biotically induced DO changes (Fig. 6).The relative contribution of the different processes and the resulting total simulated oxygen change (Fig. 6) are consistent with results from earlier model studies (Bopp et al., 2002;Matear and Hirst, 2003).
While our model results of increasing oxygen levels in the tropical thermocline under global warming can also be found in other models and do not seem to be a specific feature of the model used here, we still have to explain why simulated and observed patterns of oxygen change are anticorrelated.Given the relatively uniform contribution of solubility changes and the very small changes in biotic oxygen consumption, we expect that the model-data differences must to a large extent be attributed to errors in the modelled transport pathways and processes.A number of processes of potential importance for the control of oxygen levels in the tropical ocean's thermocline are not properly represented in the current model.These include equatorial jets, which cannot be resolved by the coarse-resolution model and which also are notoriously difficult to simulate even in high-resolution models (Brandt et al., 2008;Eden and Dengler, 2008).Sensitivity studies with a somewhat refined lateral grid resolution in the tropics (from 1.8 to 0.9 degrees meridionally), however, showed essentially no change in the simulated DO trends (Somes, C., personal communication, 2011).Other sensitivity studies, that employed enhanced zonal diffusion in the tropical Fig. 6.Top panel: Zonally averaged change in simulated oxygen at 300 dbar (black curve as in Fig. 5), the contribution due to changes in solubility diagnosed from simulated changes in temperature and salinity on the 300 dbar surface (blue), changes of simulated biotic oxygen consumption on the 300 dbar surface (green) and transport processes, which are computed as the residual between total oxygen changes and the sum of solubility and biotically induced oxygen changes (red).Bottom panel: Zonally averaged change in simulated dissolved oxygen between 1960 and 2010.Units are µmol kg −1 yr −1 .thermocline and sub-thermocline waters in order to mimick some effects of the equatorial jets, found some improvement in DO mean fields and also in trends, without changing the overall results presented here (Getzlaff, J. and Dietze, H., personal communication, 2012).Significantly finer grid resolutions may still lead to systematic changes in the simulated DO fields and, possibly, in their sensitivity to climate change.Another process not included in the standard configuration of the model using climatological winds, is the impact of interannual to decadal changes in the wind forcing.
To test the potential effect of variations in the wind forcing on oxygen concentrations, we performed additional experiments employing interannually varying monthly wind stresses of the Co-ordinated Ocean-ice Reference Experiments (CORE, Large and Yeager, 2004) for the period 1958-2007.The CORE and NCEP/NCAR mean wind stress fields show some difference with somewhat stronger westerlies over the Southern Ocean in the NCEP/NCAR climatology and somewhat stronger westerlies in the northern hemisphere mid latitude in the CORE wind stresses (Fig. 7).It cannot be assumed that either of these data sets gives a good represen-Fig.7. Zonally averaged annual-mean zonal wind stress (in dyn cm −1 ) for the standard model employing the NCEP/NCAR climatology (black), the CORE forcing (Large and Yeager, 2004) averaged over the period 1958-2007 (red) and the CORE forcing averaged over 1958-1967 (green).
tation of wind patterns prior to 1958, a time period that will, however, still impact on the oxygen fields simulated for at least part of the 1960-2010 time period under investigation here.Therefore, two sensitivity experiments are performed with different wind forcing from year 1765 to the end of 1957: experiment CORE-1 uses the NCEP/NCAR climatology until the end of year 1957 and then abruptly switches to the monthly core winds, whereas experiment CORE-2 employs, for the period 1765-1957, the monthly-mean CORE winds averaged over the decade 1958-1967 representative of conditions near the beginning of the time period under investigation (green curve in Fig. 7).
Both CORE-1 and CORE-2 simulations reveal a much larger spatial variance in simulated DO trends (Fig. 8a and b) as well as a larger uncertainty of the estimated trends because of the wind-induced interannual variability (Table 1).Compared to the standard simulation with climatological wind forcing, the use of interannually varying winds leads to a negative trend in simulated DO in the thermocline of the eastern tropical Pacific and also the eastern tropical Atlantic in both CORE-1 and CORE-2 experiments.However, both model simulations display an essentially opposite behavior in the western part of the tropical Pacific and in large parts of the Indian Ocean.Overall, patterns of DO trends among CORE-1 and CORE-2 simulations vary considerably, although both models differ only in the wind forcing applied during the years 1765 to 1957, i.e. prior to the period investigated.Simulated DO trends thus depend crucially on the poorly-known wind forcing prior to the analysis period.None of the two simulations using interannually varying CORE winds yields a good pattern correlation with the observational trend estimate (Table 1).While the overall performance with respect to DO trends of the model configurations forced by interannually varying winds does not seem to be better than that of the standard configuration using climatological wind  8c and d).
The above results indicate that neither parameter changes reflecting typical uncertainties in diapycnal mixing, the consideration of potential CO 2 effects on C : N stoichiometry, nor interannual variations in the wind forcing can bring the modelled patterns of 300 dbar DO changes into reasonable agreement with the observational estimate.Still, the simulated decline of DO at 300 dbar is, on average, smaller but, given the observational uncertainties, relatively similar to the observed decline (0.066 ± 0.051 µmol kg −1 yr −1 of the observational estimate; 0.034/0.027/0.41µmol kg −1 yr −1 for constant C : N ratios and low/standard/high vertical diffusivities) during the period 1960-2010.By confronting the simulated oxygen changes with observed multi-decadal trends, we can here demonstrate that none of the model configurations can satisfactorily reproduce the observed patterns of DO changes in the thermocline of the tropical ocean.However, given the uneven distribution of oxygen data in space and time (Fig. 1) the data interpolation method might have led to deviations from the real trends.The test by subsampling the results of our numerical standard model configuration at the times and locations of the real observations (Fig. 4b) resulted in some deviations from the complete model output (Fig. 4a), especially in the Southern Ocean, where data are sparse (Fig. 1).Future work will investigate the role of different wind patterns prior and during the analysis period in more detail.

Conclusions
Open-ocean DO trends were quantified from observations covering the past 50 yr as well as an about 70-year period in the Atlantic Ocean.In summary, a DO decrease is the dominant oxygen change within the tropical regions during last 50 yr in observations.The subtropical gyres, in contrast, show areas with increasing DO levels as well as regions with a DO decline.The zonal mean DO trend at 300 dbar is positive in the subtropics between 55 • S and 20 • S and 20 • N to 35 • N (Fig. 5).The DO increase in the subtropical gyres might be related to changes in the basin's largescale circulation (Deutsch et al., 2005) or an acceleration of the subtropical gyre in the South Pacific (e.g.Roemmich et al., 2007).The subpolar regions in the Northern Hemisphere show mainly decreasing oxygen (Fig. 2b), while in the Southern Ocean positive as well as negative trends appear.Historical sea-level records reveal that a strengthening of the Pacific subtropical cells (STCs) since the early 1990s has reversed a multi-decadal weakening tendency.It appears that neither the trends prior to the early 1990s nor those after the early 1990s can be fully attributed to greenhouse gas forcing (Feng et al., 2010).The mean DO trend between 50 • S and 50 • N at 300 m depth for the period 1960 to 2010 is −0.066 µmol kg −1 yr −1 , as the strong negative DO trends in the tropics counterbalance the areas with positive DO trends.Given the non-uniform data coverage in space and time (Fig. 1), this mean trend is generally consistent with, though somewhat larger than, the 100-1000 m trend of −0.042 µmol kg −1 yr −1 estimated by Helm et al., (2011) for the period 1970 to 1992.
Changes on a pressure layer could be influenced by density changes.However, we performed a computation also for the trends on a density surface located at about 300 dbar in the tropical regions, and the results for the observations as well as for the model run were similar to the results on the pressure layer.The simulated density change is almost completely explained by warming.All sensitivity runs using different mixing intensities show decreasing densities at 300 dbar (Fig. 9a).In the tropical Pacific, the simulated density change is almost completely explained by warming.We analyzed this point further by computing the density trend from the same bottle data used to in the observational estimate of 300 dbar oxygen trends (Fig. 9b).It turns out that observed density shows increasing trends particularly in the western part of the tropical Pacific.This is consistent with observed temperature decreases in this region (Harrison and Carson, 2007).Hence density trends at 300 dbar differ between model and observations, but as they are primarily related to warming and with solubility changes contribution little to the total DO trends (Fig. 6), these cannot explain the difference between observed and simulated oxygen trends.
DO changes in model runs are evaluated for the same period and compared against observations in order to investigate the potential reliability of such models for projecting future change scenarios.Depending on the model configuration the DO trend at 300 dbar varies between −0.047 and −0.027 µmol kg −1 yr −1 , which is smaller -albeit consistent within the error bars -than the observed mean changes.An unexpected result of our study is that simulated and observed patterns of thermocline DO change are essentially anticorrelated.The model results for the last 50 yr differ from the observations primarily in the tropical oceans, where the model shows predominantly increasing oxygen levels in contrast to the decreasing trend estimated from the observations.The weaker oxygen decrease or even increase in the tropics is a common feature in coarse-resolution biogeochemicalclimate models (Bopp et al., 2002;Matear and Hirst, 2003).
Here, we show that this discrepancy is not significantly reduced for substantially lower levels (or higher) of sub-grid scale mixing that was previously suggested as a main cause for the model-data discrepancy (Keeling et al., 2010).These differences may instead be caused by the coarse models' inability to resolve the complicated set of equatorial current bands, which are believed to be a major component in the oxygen balance of the OMZs (Stramma et al., 2010b).Also mesoscale eddies are a mechanism for large scale variability of chlorophyll (Chelton et al., 2011) as well as for oxygen, and eddy transport is not resolved in the model runs.
Interannual to decadal changes in the wind forcing, which is not included in our standard models but is accounted for in two sensitivity simulations, has been demonstrated to affect our estimates of long-term trends and make its attribution to global warming alone problematic.Our sensitivity simulation with interannually varying reanalysis winds indicates that relatively small changes in wind patterns may have a large impact on the regional patterns of DO changes, and that the poorly known wind history prior to the start of the analysis period (1960-2010) can substantially impact the DO trends simulated by the models.On the observational side the data include multi-decadal variability and future global scale measurements will have to be undertaken to distinguish (multi-) decadal variability from long-term trends.Further sensitivity studies using observed changes in wind forcing, as well as simulations using higher spatial resolution models are required to investigate these model-data discrepancies further.

Fig. 2 .
Fig. 2. Global mean dissolved oxygen (DO) distribution in µmol kg −1 (a) at 300 dbar and (b) DO changes at 300 dbar in µmol kg −1 yr −1 for the period 1960 to 2010.In (b) areas with DO increase are red and decreases blue and areas with DO changes below the 95 % confidence standard error interval are white.

Fig. 3 .
Fig. 3. Dissolved oxygen (DO) differences at 300 m depth in µmol kg −1 yr −1 in the tropical and South Atlantic Ocean between the stations (circles) collected during the 1925 to 1927 Meteor expedition and the period 1990-2008.Positive (red) values represent a DO increase in time.The background color field is the 1960 to 2010 DO trend at 300 dbar from Fig. 2b.

Fig. 4 .
Fig. 4. Simulated linear trend of dissolved oxygen at 300 dbar for the period 1960 to 2010 (a) for the standard configuration of Duteil and Oschlies (2011) with a diapycnal background diffusivity of k bg = 0.15 cm 2 s −1 and constant (Redfield) C : N stoichiometry, (c) for a much reduced background diffusivity of k bg = 0.01 cm 2 s −1 and constant C : N stoichiometry, and (d) for increased background diffusivity of k bg = 0.3 cm 2 s −1 .(b) shows the trend estimated for the model results subsampled according to the dates and locations of the read DO measurements and using the same routine that was applied to the in situ data in Fig. 2. Units are µmol kg −1 yr −1 .

Fig. 8 .
Fig. 8. Simulated linear trend of dissolved oxygen at 300 dbar for the period 1960 to 2010 (a) for the experiment CORE-1 switching from NCEP/NCAR climatological to interanually varying CORE winds on 1. January 1958, (b) running the model with 1958-1967 mean CORE forcing from 1765 until the end of 1957 and then switching to interannually varying CORE wind stress fields.DO changes at 300 dbar simulated with the climatological NCEP/NCAR wind stress are shown in (c), and simulated with the CORE 1958-1967 monthly mean wind stress in (d).Units are in µmol kg −1 yr −1 .

Table 1 .
Mean trend and standard deviation in dissolved oxygen and pattern correlation with observational estimate (wherever exceeding the 95 % threshold) simulated by the various model configurations.Standard model configurations (i) to (iii) employ climatological seasonally cycling monthly winds and Redfield C : N stoichiometry.Only configuration (iv) employs variable C : N stoichiometry.Model experiments (v) and (vi) use Co-ordinated Ocean-ice Reference Experiments (CORE) monthly wind stress fields for the period 1958-2007 with two different spin-up wind fields prior to 1958 (see text).