Multiple stressors of ocean ecosystems in the 21st century: projections with CMIP5 models

. Ocean ecosystems are increasingly stressed by human-induced changes of their physical, chemical and biological environment. Among these changes, warming, acid-iﬁcation, deoxygenation and changes in

Abstract.Ocean ecosystems are increasingly stressed by human-induced changes of their physical, chemical and biological environment.Among these changes, warming, acidification, deoxygenation and changes in primary productivity by marine phytoplankton can be considered as four of the major stressors of open ocean ecosystems.Due to rising atmospheric CO 2 in the coming decades, these changes will be amplified.Here, we use the most recent simulations performed in the framework of the Coupled Model Intercomparison Project 5 to assess how these stressors may evolve over the course of the 21st century.The 10 Earth system models used here project similar trends in ocean warming, acidification, deoxygenation and reduced primary productivity for each of the IPCC's representative concentration pathways (RCPs) over the 21st century.For the "businessas-usual" scenario RCP8.5, the model-mean changes in the 2090s (compared to the 1990s) for sea surface temperature, sea surface pH, global O 2 content and integrated primary productivity amount to +2.73 (±0.72) • C, −0.33 (±0.003) pH unit, −3.45 (±0.44) % and −8.6 (±7.9) %, respectively.For the high mitigation scenario RCP2.6, corresponding changes are +0.71(±0.45) • C, −0.07 (±0.001) pH unit, −1.81 (±0.31) % and −2.0 (±4.1) %, respectively, illustrat-ing the effectiveness of extreme mitigation strategies.Although these stressors operate globally, they display distinct regional patterns and thus do not change coincidentally.Large decreases in O 2 and in pH are simulated in global ocean intermediate and mode waters, whereas large reductions in primary production are simulated in the tropics and in the North Atlantic.Although temperature and pH projections are robust across models, the same does not hold for projections of subsurface O 2 concentrations in the tropics and global and regional changes in net primary productivity.These high uncertainties in projections of primary productivity and subsurface oxygen prompt us to continue intermodel comparisons to understand these model differences, while calling for caution when using the CMIP5 models to force regional impact models.

Introduction
In recent decades, the ocean has undergone large physical and biogeochemical modifications in response to humaninduced global change, as revealed by a variety of in situ and remote sensing observations.These changes encompass ocean surface warming, changes in ocean salinity, modifications of density structure and stratification, as well as an increase in dissolved inorganic carbon concentrations and a decrease in seawater pH in response to ocean uptake of anthropogenic carbon (Doney, 2010).These physical and chemical modifications (Bindoff et al., 2007) have the potential to affect marine organisms and ecosystems (Doney et al., 2012).
Temperature for instance has a fundamental effect on biological processes, as most biological rates (e.g., enzyme reactions) are temperature-dependent (e.g., Eppley et al., 1972).Sea surface temperature has increased by about +0.7 • C over the last 100 yr (Bindoff et al., 2007).Some organisms may be able to acclimate to temperature changes that do not push them outside of their optimal range, but rising temperatures this century could cause poleward shifts in species thermal niches and a decline in diversity as shown by Thomas et al. (2012) for tropical phytoplankton.In addition, rising ocean temperatures lead to increased ocean stratification and changes in ocean mixing and ventilation that could induce indirect effects, which in turn reduce subsurface dissolved O 2 concentrations (Plattner et al., 2001) and supply of nutrients to the euphotic layer, increasing nutrient stress for phytoplankton (Bopp et al., 2001).
Sea surface pH has declined by about 0.1 pH unit since pre-industrial times (Bindoff et al., 2007).This decrease, along with other large-scale changes in seawater chemistry (increased dissolved CO 2 , decreased carbonate ion concentrations), is referred to as anthropogenic ocean acidification and is a direct consequence of the oceanic absorption of ∼ 30 % of the total anthropogenic emissions of carbon dioxide (CO 2 ) since 1750 (Sabine et al., 2004;Khatiwala et al., 2009).These changes may have adverse effects on marine organisms.For instance, lower carbonate ion concentrations reduce calcification for many calcifying organisms (Hofmann et al., 2010), while very high CO 2 levels impose increased physiological stress on heterotrophic organisms (Pörtner et al., 2008) and may lead to decreased stress on autotrophic organisms through CO 2 fertilization of photosynthesis, particularly for slow growers such as diazotrophs (Hutchins et al., 2007).
Reduced oceanic O 2 concentrations have been reported for coastal areas in response to human-induced eutrophication (Gilbert et al., 2010), and for the open ocean as a consequence of two major processes: a warming-induced reduction in solubility and increased stratification/reduced ventilation (Helm et al., 2011).This reduction in O 2 , referred to as ocean deoxygenation, between the 1970s and the 1990s amounts to a global average reduction of −0.93 mmol m −3 (Helm et al., 2011).Most organisms are not very sensitive to O 2 levels as long as concentrations are adequate.However, beyond a critical threshold referred to as hypoxia (below 50-80 mmol m −3 depending on taxa; see Vaquer-Sunyer and Duarte, 2008 for relevant references), many organisms start to suffer from several physiological stresses that could ultimately lead to death.
In addition, net production of organic matter by phytoplankton (net primary production, or NPP), which represents the first link in open-ocean marine food webs, may also be impacted by climate change and climate variability (Behrenfeld et al., 2006).Boyce et al. (2010) reported a global decrease of phytoplankton biomass of 1 % per year since 1900 based mostly on ocean transparency measurements.However, this result is debated and might arise from a temporal sampling bias (Rykaczewski and Dunne, 2011).Observations also suggest an extension of oligotrophic gyres by ∼ 15 % between 1998 and 2007 (Polovina et al., 2008).Conversely, Henson et al. (2011) emphasize that this trend cannot be unequivocally attributed to the impact of global climate change but may rather be a result of the strong La Niña in 1998, and that a much longer time series is needed to distinguish a global warming trend from natural variability.Whether or not these changes in primary productivity are attributable to global warming, they are likely to be felt further up the food chain, with potential implications on marine resources (Pauly and Christensen, 1995;Nixon and Thomas, 2001).
Marine ecosystems are adapted to particular environmental conditions.Temperature, pH, O 2 and NPP are fundamental controlling factors, which, when changed, affect how an ecosystem will be altered (Fischlin et al., 2007).Several previous modeling studies used coupled climate-marine biogeochemical models to investigate how these stressors could evolve under climate change scenarios.Most of these modeling studies will be discussed later when compared to our new results, but they all project large changes in these marine ecosystem stressors over the coming decades as described briefly below.
Using 13 ocean-only models, Orr et al. (2005) showed that global-mean surface pH could drop by 0.3-0.4pH unit in 2100 under the IS92a scenario (a "business-as-usual" scenario less severe than the RCP8.5 scenario used in this study).In this scenario, Southern Ocean surface waters start to become undersaturated with respect to aragonite as soon as year 2030 (McNeil and Matear, 2008).In 2100, aragonite undersaturation extends throughout the entire Southern Ocean and into the subarctic Pacific.
Coupled climate-marine biogeochemical models used over the past 15 yr all project a long-term decrease in the oceanic O 2 inventory in response to anthropogenic global warming (e.g., Sarmiento et al., 1998;Plattner et al., 2001;Bopp et al., 2002).In a recent intercomparison of seven Earth system models, Cocco et al. (2013) found that the oceanic O 2 inventory would decline by 2-4 % in 2100 under SRES-A2, an earlier IPCC scenario similar to the RCP8.5 scenario used here.
Finally, most coupled climate-marine biogeochemical models also project a decline in NPP in the coming decades as a response to anthropogenic global warming (e.g., Bopp et al., 2001).In a recent model intercomparison study of four coupled models, Steinacher et al. (2010) reported a decrease in global mean NPP of 2-20 % by 2100 relative to preindustrial conditions in the SRES A2 emission scenario.
These physical and biogeochemical changes in temperature, pH, O 2 and NPP are known to interact with each other, potentially leading to synergistic effects, both on biogeochemical cycles and on marine ecosystems (Gruber, 2011).Some of these synergistic effects on biogeochemical cycles occur at large spatial (regional, global) scales, one example being the still-debated impact of ocean acidification on ocean deoxygenation.Using ocean biogeochemical models and a simple parameterization based on mesocosm experiments and relating C / N ratios of organic matter to CO 2 levels (Riebesell et al., 2007), Oschlies et al. (2008) and Tagliabue et al. (2011) showed that increasing dissolved CO 2 could induce large increases in subsurface O 2 utilization, hence expanding the volume of suboxic waters.Synergistic effects also take place at the physiological level.Temperature and dissolved CO 2 may affect levels of tolerance to low-O 2 concentrations (Pörtner et al., 2004;Pörtner et al., 2007), whereas elevated CO 2 and lower O 2 levels may reduce thermal tolerance of some organisms (Pörtner, 2010, Metzger et al., 2008).These studies emphasize that multiple stressors should be studied together simultaneously in order to be able to evaluate synergistic effects.
Here, we use the most recent simulations performed in the framework of the Coupled Model Intercomparison Project 5 (CMIP5, Taylor et al., 2012) to assess how these stressors may evolve over the course of the 21st century.We focus on four stressors: temperature, pH, O 2 and NPP, while comparing four different representative concentration pathways (RCPs) scenarios across 10 different Earth system models (ESMs), all including a marine biogeochemical component.The choice of pH as one of the main stressors instead of calcite/aragonite undersaturation states or carbonate ion concentrations has been dictated by the fact that pH is a more generic variable of ocean acidification and thus concerns many more processes than solely calcification (e.g., see Hutchins et al., 2009, for the impact of pH changes on marine nutrient cycles).Note that some of the interactions between these stressors are not included in the models used here and could still well cause further changes.
The remainder of the paper is organized as follows.Section 2 describes the models and the set of simulations used in this study.Section 3 presents the main results and a discussion, describing major evolution of the different stressors at the global scale, at the regional scale and relations between the stressors.Conclusions and perspectives follow in Sect. 4.

Models and simulations
The latest generation of Earth system models (ESMs) were used to carry out simulations within the framework of CMIP5 (Taylor et al., 2012).These simulations include four future scenarios referred to as RCPs (Moss et al., 2010;van Vuuren et al., 2011): RCP8.5, RCP6.0,RCP4.5 and RCP2.6.The RCPs are labeled according to the additional radiative forcing level in 2100 with CO 2 concentrations reaching 936, 670, 538 and 421 ppm, respectively.RCP2.6 is also referred to as RCP3PD for "peak and decline": the atmospheric CO 2 peaks at a concentration of 443 ppm in 2050 before declining in the second half of the 21st century.
The selection of the 10 models used for this study was based on the availability of all variables necessary to discuss the four stressors we focus on: temperature, pH, O 2 and NPP.They were also selected on the requirement that at least one representative concentration pathway (RCP) was performed up to 2099.We used historical simulations from 1870 to 2005, all climate change scenarios (RCPs) from 2006 to 2099, and pre-industrial control simulations.Historical and climate change scenarios are forced not only by prescribed atmospheric CO 2 but also by other greenhouse gas and aerosols concentrations, anthropogenic land use evolution, as well as by natural forcings such as solar and volcanic forcings.In case several ensemble members were run for the same scenario and with the same model, only the first member was used.Note that the simulations used here are simulations in which atmospheric CO 2 is prescribed (as in Jones et al., 2013), and not the ones also performed with several of these Earth system models, in which atmospheric CO 2 is explicitly computed from prescribed anthropogenic emissions and the simulated ocean and land carbon sources and sinks (as in Friedlingstein et al., 2013).We did not use these latter simulations, as different atmospheric CO 2 levels between the different models for the same scenario would have added another degree of complexity in our model intercomparison study.
Table 1a and b present the 10 models used for this study.Each includes representations of the general circulation and physics of the atmosphere and the ocean (Table 1a), as well as biogeochemical components (Table 1b), including a representation of the ocean carbon cycle and the lowest trophic level of marine ecosystems (i.e., phytoplankton).However, the models differ in many respects.Thus attributing the causes of differences between models to particular processes or parameterizations is difficult.
Atmospheric and ocean resolutions vary widely across the different models (Table 1a).Typical atmospheric horizontal grid resolution is ∼ 2 • , but it ranges from 0.94 to 3.8 • .Typical oceanic horizontal resolution is ∼ 1 • , but it ranges from 0.3 to 2 • .The number of vertical levels varies from 24 to 95 in the atmosphere and 31 to 63 in the ocean.All marine biogeochemical components are typical nutrient-phytoplankton-zooplankton-detritus (NPZD) models, but with varying degrees of complexity, illustrated for instance by the number of phytoplankton functional groups (from 1 to 3) or limiting nutrients (from 3 to 5) explicitly represented (Table 1b).Note that while these models may be www.biogeosciences.net/10/6225/2013/Biogeosciences, 10, 6225-6245, 2013 well suited to study the evolution of the four considered stressors (temperature, pH, dissolved O 2 and NPP), they are not designed to study the impact of changes in these stressors on marine ecosystems per se.For instance, most of these models do not include any effect of pH, dissolved CO 2 or carbonate ion concentration on phytoplankton photosynthesis or calcification.

Model output, model mean and robustness
In this study, we make use of the standard CMIP5 output from the Program for Climate Model Diagnosis and Intercomparison (http://pcmdi3.llnl.gov/esgcet/home.htm).We use output provided by the different models listed in Table 1a and available at the time of writing.Not all models have performed all RCP simulations.The variables selected for this study are annual-mean 3-D temperature, salinity, meridional velocity, pH, carbonate ion, O 2 as well as 2-D vertically integrated NPP and export production of organic particles at 100 m.For models that did not provide pH or carbonate ion concentration fields to the database at the time of this writing (MPI-ESM-LR and HadGEM2-ES for pH, MPI-ESM-LR for carbonate ion), we computed pH and carbonate ion values using output of 3-D Dissolved Inorganic Carbon (DIC), alkalinity, temperature, and salinity and and routines based on the standard OCMIP carbonate chemistry adapted from earlier studies (Orr et al., 2005).
To be able to compute model means and inter-model standard deviations, all variables were interpolated onto a common 1 • × 1 • regular grid and to standard ocean depths, using a Gaussian weighted average.Preindustrial control simulations, for which atmospheric CO 2 and other greenhouse gases and aerosols are set to preindustrial levels, were used to remove potential century-scale model drifts.This correction is applied for all regional-scale analyses and for global-average time series of NPP, export production, O 2 and heat content.It is not applied for time series of surface pH and surface sea temperature for which no apparent long-term drift was detectable.
In the following analysis, we use both individual model results (for RCP8.5) and model means (for all scenarios).In the case of model means, no weighting functions are applied: all models contribute the same to the model mean.
Our core analysis is based on sea surface temperature, sea surface pH, vertically integrated NPP, and O 2 averaged over 200-600 m.For O 2 , we focus first on this subsurface depth interval to be able to describe changes in and around the oxygen minimum zones as well as changes in mid-latitude ventilated waters.Later, we complete the core analysis by extending it to make use of the full 3-D distribution of temperature, pH and dissolved O 2 , as well as by using additional variables such as export production and carbonate ion concentrations.
The inter-model difference or model spread is used as an estimate of uncertainty around the projections.To show model agreement on global-mean time series, we use one inter-model standard deviation.To show model agreement at the regional scale, we use stippling on the maps based on a simple robustness (or agreement) measure.Similar to the approach used for surface temperature changes in Meehl et al. (2007), high robustness for sea surface temperature and for pH is defined when the model-mean-simulated change exceeds the inter-model standard deviation (the robustness index, defined as the ratio between model-mean-simulated Table 1b.Main characteristics of the 7 marine biogeochemical components of the 10 ESMs used in this study: list of nutrients limiting phytoplankton growth, the number of explicit phytoplankton groups represented, the number of explicit zooplankton groups represented, explicit or implicit representation of heterotrophic bacteria, indications of the use of fixed Redfield (R) or variable (V) ratios for organic matter production, and Q 10 for temperature dependency of biogeochemical processes (autotrophic/heterotrophic when 2 numbers are given).change and inter-model standard deviation, is then larger than 1).Similar to the approach used for precipitation changes in Meehl et al. (2007), high robustness for oxygen and NPP is defined when at least 80 % of models agree on the sign of the mean change (robustness index larger than 1).Note that the estimate of uncertainty based on model spread may be biased by an arbitrary distribution of CMIP5 model output for a specific variable and by similarities between models (e.g., IPSL-CM5A-LR and IPSL-CM5A-MR share the same components but differ in atmospheric resolution; see Table 1a).

Water mass analysis
We used a global framework to group together water masses of similar behavior.Four classes were defined: tropical water (TW) masses, mode and intermediate water (MIW) masses, deep water (DW) masses and bottom water (BW) masses.For example, the class MIW aims at gathering mode and intermediate waters of all basins, which share common features but are not distributed in the same range of density (Hanawa and Talley, 2001).Limits between classes were defined using several criteria (salinity, stratification, meridional velocities) resulting in different density thresholds in the five different basins (North Atlantic, South Atlantic, North Pacific, South Pacific and Indian Ocean) and for the different models.These density thresholds were computed for each model using the first 10 yr (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) of scenario RCP4.5 (tem-perature, salinity, etc.) and the density referenced to 2000 m (σ 2 ).
The limit between well-stratified TW and homogeneous MIW was defined using a stratification criterion (∂σ 2 /∂z > 0.02 kg m −4 , with z the depth).The lower boundary of MIW was defined as the depth where the salinity is equal to its deep minimum value +0.05.DWs were distinguished from BW using the deepest change in sign of meridional velocities, orientated northward in BW and southward in DW.Note that in the North Atlantic MIWs are not associated with a deep salinity minimum.Instead the limit between MIW and DW was defined as a minimum in meridional velocities, which are oriented southward in both water masses.

Results and discussion
Most models, including their marine biogeochemical components, have been evaluated individually elsewhere (see references in Table 1a).Here, we briefly discuss global mean values of present-day sea surface temperature (SST), surface pH, oceanic O 2 content and integrated NPP, as listed in Table 2 for each individual model.We also show a regional comparison to data-based products of present-day modelmean SST, surface pH, subsurface O 2 and NPP (Fig. 1).Finally the skill of the different models in representing spatial variability of data-based fields is quantified while relying on Taylor diagrams (Taylor, 2001, Fig. 2), for which we extend the model-data comparison to 3-D temperature, 3-D pH, 3-D dissolved O 2 as well as surface and 3-D carbonate ion concentrations fields.
Whereas large-scale patterns are well represented for SST and subsurface O 2 , the comparison to data-based products is much less satisfying for surface pH and NPP (Fig. 1).This is also reflected in the Taylor diagrams, with correlation coefficient (R) ranging from 0.98 to 1.0 for SST, from 0.7 to 0.95 for O 2 , but from 0.1 to 0.6 for surface pH and from 0.2 to 0.6 for NPP (Fig. 2).For NPP, the models score poorly because we do not include the seasonal cycle in the correla-   tion, and because they under-represent high productivity associated with coastal regions, as discussed for a previous set of Earth system models in Schneider et al. (2008).For surface pH, the models score poorly because of the very uniform nature of this field (Fig. 1).When turning to 3-D pH, or to carbonate ion concentrations (surface and 3-D), the modeldata correlation coefficients increase to 0.6-0.85,0.93-0.98 and 0.85-0.96,respectively.Additionally, we can note that models that do perform well for one variable may perform Table 2. Observed and modeled present-day global average of sea surface temperature ( • C), surface pH (pH unit), dissolved oxygen (mmol m −3 ), volumes of waters (10 15 m 3 ) in which O 2 is less than 80, 50 and 5 mmol m −3 , net primary productivity (PgC yr −1 ) and export production of organic particles at 100m (PgC yr −1 ).Interannual standard deviations, when meaningful, are also indicated in parenthesis.Observed estimates are from Reynolds et al. (2008) for SST, computed using DIC and alkalinity from Key et al. (2004) for surface pH, from WOA (2009) for O 2 content, from Bianchi et al. (2012) for the volumes of low oxygen waters, from Behrenfeld et al. (1997)  poorly for others (e.g., IPSL-CM5A-MR performs well for NPP (R = 0.6), but poorly for surface pH (R = 0.2).
Global mean values of present-day SST, surface pH, oceanic O 2 content, and integrated NPP also show some striking differences between models, and when compared to observations.This is especially true for the globally averaged O 2 concentration, with some models clearly underoxygenated (e.g., IPSL-CM5As) and other models overoxygenated (e.g., NorESM1-ME).For NPP also, some models simulate global integrated values as low as 30.9Pg C yr −1 (IPSL-CMA-MR), whereas others simulate NPP as high as 78.7 Pg C yr −1 (GFDL-ESM2M).While these model differences in reproducing present-day patterns and values may explain some of the differences in the model projections we detail below, they also lead us to use relative quantities when comparing NPP changes or O 2 changes, as done in the rest of the manuscript.

Changes of multiple stressors at the global scale
The ocean warms because it takes up much of the additional heat that accumulates in the Earth system due to increas-ing greenhouse gas concentrations.The intensity of simulated sea surface warming in the coming decades is mostly dictated by the RCP scenario, i.e., by the amount of greenhouse gases emitted to the atmosphere, with an inter-model range depending on the strength of the simulated climate feedbacks.From the 1990s to the 2090s, model-mean global average SST increases by +2.73 (±0.72), +1.58 (±0.48), +1.28 (±0.56) and +0.71 (±0.45) • C for RCP8.5, RCP6.0,RCP4.5 and RCP2.6, respectively (Fig. 3, Table 3).Note that because of our model selection process, these model-mean values would differ from the standard CMIP5 analyses that include a wider selection of models.This simulated increase is, as expected, lower than for global-mean air surface temperature, which amounts to +4.2, +2.5, +1.9 and +1.0 • C for the same four RCP scenarios and from 1960-1990 to the end of the 21st century (Knutti and Sedlacek, 2012).
The model spread for each scenario is used as an estimate of uncertainty around the model-mean projection.Part of this model spread is due to internal variability simulated by the climate models.However, most of it arises from model differences in (1) the way RCP scenarios are set up (i.e., aerosols  and some greenhouse gases concentrations other than CO 2 may differ between models for the same scenario; Szopa et al., 2012), and in (2) climate sensitivities (Knutti and Hegerl, 2008).For example, the SST warming for the RCP8.5 scenario reaches +3.5 • C in three of the models (MPI-ESM-LR, IPSL-CM5A-LR and IPSL-CM5A-MR) and only 2.25 • C in two others (GFDL-ESMs) (Fig. 4).This has been explained by the differences in climate sensitivity: IPSL-CMs and MPI-ESM-LR have high 2 × CO 2 equilibrium climate sensitivities, whereas the GFDL-ESMs are on the low range of climate sensitivities as demonstrated by Andrews et al. (2012).
Sea surface pH decreases as a consequence of the ocean taking up a significant fraction of anthropogenic carbon accumulated in the atmosphere.Even more than for SST, the magnitude of pH decrease is entirely dictated by the scenario (for a given atmospheric CO 2 concentrations).In the 2090s, the drop in global-average surface pH compared to 1990s values amounts to −0.33 (±0.003), −0.22 (±0.002), −0.15 (±0.001) and −0.07 (±0.001) pH unit, for RCP8.5, RCP6.0,RCP4.5 and RCP2.6, respectively (Fig. 3, Table 3).The model-mean projection for RCP8.5 is slightly larger than that by Orr et al. (2005) for the IS92a scenario, in which atmospheric CO 2 reaches 712 ppmv in 2100 as compared to 935 ppm in RCP8.5.Projections of surface carbonate ion concentrations, a better variable than pH to discuss potential impacts on calcification and calcifiers, are detailed in Sect.3.2.2.In contrast to SST projections, the model spread for global surface pH projections (estimated as the inter-model standard deviation) is very low (less than 0.003 pH unit).This is explained by (1) the weak interannual variability in global mean surface pH (Fig. 4), (2) a weak climate-pH feedback, as demonstrated in Orr et al. (2005) for earlier Earth system models, (3) the similar carbonate chemistry equations and well-defined constants based on the OCMIP-2 protocol used by most, if not all, models (Orr et al., 2000) and (4) the uniqueness of the ocean acidification forcing (i.e., of the atmospheric CO 2 trajectory of each RCP scenarios; Moss et al., 2010).Indeed, changes in surface ocean pCO 2 , and hence corresponding changes in carbonate chemistry, closely track changes in atmospheric CO 2 because the equilibration time for CO 2 between the atmosphere and mixed layer is fast enough in most areas (global average of ∼ 8 months) to allow near equilibration.
All models lose O 2 from the ocean in response to climate change under every RCP scenario (Fig. 3).The modelmean reduction in global ocean oxygen content reaches −3.45 (±0.44), −2.57(±0.39), −2.37 (±0.30), and −1.81 (±0.31) % in the 2090s relative to the 1990s, for RCP8.5, RCP6.0,RCP4.5 and RCP2.6, respectively (Fig. 3 and Table 3).For RCP8.5, this translates into a −9.03 (±1.15)Tmol O 2 or −6.13 (±0.78) mmol m −3 decrease in the 2090s relative to the 1990s, based on the reference global O 2 inventory from WOA 2009.This long-term decline in O 2 inventory is a consistent trend simulated in many coupled climate-marine biogeochemical models (e.g., Sarmiento et al., 1998).Our simulated rates of deoxygenation in the CMIP5 models under RCP8.5 are quite similar to those reported in a previous model intercomparison study (−2 to −4 %) of seven models under the SRES-A2 scenario (Cocco et al., 2013).The model spread is limited (moderate uncertainty) largely because the global O 2 inventory is an integrated quantity with low interannual variability (Fig. 4).The models with higher rates of deoxygenation, such as GFDL-ESMs and CESM-BGC1, are not always those that simulate more intense sea surface warming.This suggests that the warming-induced reduction in solubility may not be the first driver of deoxygenation, at least in explaining model differences.We further detail the contribution of solubility changes to deoxygenation in Sect.3.3.1.
All models simulate a decrease in integrated global primary productivity during the 21st century, but there is a large spread.Relative changes in NPP amount to −8.6 (±7.9), −3.9 (±5.7), −3.6 (±5.7), and −2.0 (±4.1) % in the 2090s relative to the 1990s, for RCP8.5, RCP6.0,RCP4.5 and RCP2.6, respectively (Fig. 3 and Table 3).This reduction is consistent with a recent inter-model comparison in which four models all simulate a decline in integrated NPP, ranging from −2 % to −13 %, from 1860 to 2100 under the SRES A2 scenario (Steinacher et al., 2010).Those simulated decreases were primarily attributed to reduced nutrient supply to the surface in the tropics and temperate ocean due to enhanced stratification and reduced ocean ventilation.However other modeling studies find that integrated NPP increases during the 21st century, even though simulated new production and export production still decrease (Schmittner et al., 2008;Taucher and Oschlies, 2011).The decoupling between export production and net primary production in these latter studies has been shown to be driven by the direct effects of temperature on biological processes such as remineralization rates.These links are explored further in Sect.3.2.3.
The large inter-model spread in NPP is obvious in Fig. 3 for all scenarios and illustrated for RCP8.5 in Fig. 4.While this large spread is partly linked to the large interannual variability of NPP, as simulated by some models (e.g., HadGEM2-ES, NorESM1-ME and MPI-ESMs), it is mainly due to large model differences.Some models (GFDL-ESMs) simulate a slight decrease of globally integrated NPP (between 0 and −2 % in 2100), whereas others (MPI-ESMs, and HadGEM2-ES) simulate a large decrease (between −12 and −18 % in 2100).The fact that models of the same family (IPSL-CM5A-LR and IPSL-CM5A-MR, MPI-ESM-LR and MPI-ESM-MR, GFDL-ESM2G and GFDL-ESM2M, respectively), which share the same biogeochemical component, project similar changes in globally integrated NPP (Fig. 4) suggests that differences in marine biogeochemical model parameterizations explain model differences, at least from a global perspective.

General description
Unlike global mean projections, simulated changes in temperature, pH, O 2 and NPP have contrasting regional patterns.Figures 5 and 6 show the changes of SST, surface pH, averaged O 2 concentrations over 200-600 m and vertically integrated NPP, between 1990s and 2090s, and for the two extreme scenarios, RCP8.5 and RCP2.6, respectively.Local model agreement based on a robustness index detailed in Sect.2.2 is indicated by stippling.
Changes in SST are not spatially uniform.Stronger warming occurs in the tropics, in the North Pacific and in the Arctic Ocean, with SST increases larger than 4 • C in the RCP8.5 scenario (Fig. 5).On the contrary, much weaker warming, even cooling, is simulated in the North Atlantic and in some parts of the Southern Ocean, where deep convection is strongly reduced or where sea ice remains unchanged.Robustness of these regional projections is high (Knutti and Sedlacek 2012), even for the low-emission scenario RCP2.6 (Fig. 6).Only the regions with a weak signal show low robustness, which is a consequence of the signal-to-noise metric used here to estimate robustness (the ratio of model-mean change over inter-model standard deviation).
Changes in surface pH are smoother and more uniform than SST changes, and very robust across models.Surface pH changes range from −0.25 to −0.45 pH unit in RCP8.5 and from −0.05 to −0.15 in RCP2.6 (Figs. 5 and 6).Larger changes occur in the surface Arctic Ocean for both scenarios.This is consistent with recent results obtained with the MIROC ESM under the RCP8.5 scenario (Yamamoto et al., 2012).The authors showed that Arctic sea-ice melting amplifies the decrease of surface pH due to the uptake of anthropogenic carbon, consistent with Steinacher et al. (2009).
Changes in subsurface (200-600 m) O 2 are not spatially uniform, and there is less agreement among models.But despite a strong difference in magnitude, the complex patterns of spatial changes are very similar across the two scenarios and reflect the influence of changes in several processes (ventilation, vertical mixing, remineralization) on O 2 levels (Figs. 5 and 6).The North Pacific, the North Atlantic, the Southern Ocean, the subtropical South Pacific and South Indian oceans all undergo deoxygenation, with O 2 decreases of as much as −50 mmol m −3 in the North Pacific for the RCP8.5 scenario.In contrast, the tropical Atlantic and the tropical Indian show increasing O 2 concentrations in response to climate change, in both RCP8.5 and RCP2.6 scenarios.The equatorial Pacific displays a weak east-west dipole, with increasing O 2 in the east and decreasing O 2 in the west.Apart from changes in the equatorial Pacific, these regional changes in subsurface O 2 are consistent across models under the RCP8.5 scenario (stippling in Fig. 5), and they are quite similar to those from a recent inter-model com-parison of the previous generation of Earth system models (Cocco et al., 2013).
Over the mid-latitudes, patterns of projected changes in subsurface O 2 are broadly consistent with observations collected over the past several decades (Helm et al., 2011;Stendardo and Gruber, 2012;Takatani et al., 2012).Yet there is no such model-data agreement over most of the tropical oceans.Observed time series suggest a vertical expansion of the lowoxygen zones in the eastern tropical Atlantic and the equatorial Pacific during the past 50 years (Stramma et al., 2008), conversely with models that simulate increasing O 2 levels with global warming over the historical period (Andrews et al., 2013).A more detailed analysis of the simulated evolution of volumes of low-oxygen waters is given in Sect.3.2.4.
Similar to subsurface O 2 and in line with previous modeling studies (Bopp et al., 2001;Steinacher et al., 2010), projected changes in NPP are spatially heterogeneous.A decrease in NPP is consistently simulated across models and scenarios in the tropical Indian Ocean, in the west tropical Pacific, in the tropical Atlantic and in the North Atlantic (Figs. 5 and 6).This decrease reaches as much as −150 g C m −2 yr −1 regionally in the 2090s for the RCP8.5 scenario, more than a 50 % decrease in historical levels of NPP in the North Atlantic, while at the same time there is a 30 % decrease in the tropical Indian and west tropical Pacific.The main mechanisms responsible for NPP decreases in the tropics and in the North Atlantic have been identified in a previous model inter-comparison study (Steinacher et al., 2010) and are linked to a reduced supply of nutrients to the euphotic zone in response to enhanced stratification and slowed circulation.In the eastern equatorial Pacific, the model mean also indicates a large decrease of NPP, but this response is not consistent across models, with three models (GFDL-ESMs and CESM1-BGC) simulating an increase in NPP in response to climate change in that region (eastern equatorial Pacific).It is worth pointing out that in some regions of the largest changes in the model-mean NPP (e.g., equatorial Pacific), the models do not agree even in the sign of changes.
Finally, the mean of the CMIP5 models exhibits increased NPP in the western North Pacific, the Arctic Ocean, and in parts of the Southern Ocean (Figs. 5 and 6).These changes are consistent across models for the RCP8.5 scenario but not for the RCP2.6 scenario.Mechanisms that could explain such an increase in NPP are related to an alleviation of light limitation and/or temperature limitation as shown by Steinacher et al. (2010), and to increased nutrient supply in upwelling regions under a shoaling nutricline (Rykaczewski and Dunne, 2010).
We now turn to four specific questions related to the changes of surface carbonate ion concentrations, the changes of pH at depth, the links between projected changes of NPP and export production, and the evolution of oxygen minimum zones.

Projected changes of carbonate ion concentrations
Whereas we used pH as one of our four major stressors (with temperature, dissolved O 2 and NPP), carbonate ion concentrations or saturation states are more relevant variables to focus on when discussing impacts of ocean acidification on calcification and calcifiers.undersaturation in the surface waters of the Southern Ocean and the North Pacific in 2100 as already shown by Orr et al. (2005) for the IS92a scenario.Even earlier aragonite undersaturation could be reached in the Arctic waters as demonstrated by Steinacher et al. (2009).

Projected pH changes over depth
Direct observations of pH from available time series show consistent trends of decreasing pH of about 0.02 pH unit per decade for the surface ocean (Bindoff et al., 2007).However, several recent studies have reported more rapid reductions of pH at subsurface (between 200 and 300 m) in the subtropical oceans (Bates et al., 2012;Dore et al., 2009;Ishii et al., 2011), which can be explained solely by a lower carbonate buffering capacity of these subsurface waters (Orr, 2011).These larger subsurface changes in pH could occur at depth relative to the surface, consistent across all models used here.In the 2090s and for the RCP8.5 scenario, maximum pH changes are located at the surface in the high latitude oceans and in upwelling regions, but at 100 to 400 m depth in the subtropics (Fig. 8).These changes of subsurface pH are 0.05 pH unit larger than the comparable surface changes (Fig. 8).The buildup of this subsurface low-pH water may have implications on the timing of pH changes at the surface when these water masses are brought/mixed back to the surface (Resplandy et al., 2013).

Projected changes of NPP and export production
Model-mean projections of changes in export production are shown for RCP8.5 in Fig. 9.These changes display very similar patterns to those of NPP (Fig. 5), with increasing export production in the Arctic Ocean and in parts of the Southern Ocean, and decreasing export production in the tropical oceans, and in the North Atlantic.Interestingly, robustness of the model-mean projection is larger for export production than for NPP projections: 38 % of 1 • × 1 • ocean pixels have a robustness index of more than 1 for projected NPP, whereas 45 % of the same pixels have a robustness index greater than 1 for projected export production.This is especially the case for the eastern equatorial Pacific, where eight out of nine models project a decrease of export production with climate change, whereas the response for NPP is split between an increase (three models) and a decrease (six models).For instance, GFDL-ESM2G and GFDL-ESM2M simulate a decoupling of export and primary production in that region (i.e., decreased export production and increased NPP).Some potential mechanisms explaining such decoupling are discussed below.
For the global average, this increased model consistency for simulated changes in export production versus those for NPP is also apparent (Fig. 9).For the RCP8.5 scenario, projected changes in the 2090s for export production reach −7 to −18 % (inter-model standard deviation is 5.9 %) , whereas for NPP they range between −2 and −16 % (inter-model standard deviation is 7.9 %).Taucher and Oschlies (2011) have explored the parallel responses of NPP and export production to climate change using a simple NPZD model in which they vary the dependence of production and remineralization rates on temperature.They show that when rates that are dependent on temperature are taken into account, elevated temperatures lead to enhanced NPP and accelerated carbon recycling, which is opposite to what is shown here with the CMIP5 models.However, whether or not temperature-dependent rates are used, simulated export production decreases with climate change in their simulations, as with previous model studies (Bopp et al., 2001).Note that almost all of the biogeochemical components of the CMIP5 models include temperaturedependent production and remineralization rates (Table 1b).It is not yet clear if different parameter values for temperature dependency or different routing functions of living organic material to particulate or dissolved pools may explain the large discrepancy in NPP projections.Whereas global export production seems to be mostly controlled by the balance between reduction in nutrients supplied to the euphotic layer and alleviation of light and temperature limitations (but also by changes in plankton community as shown by Bopp et al., 2005), larger uncertainties in projections of NPP emphasize the need to improve model representations of more direct biological effects, such as the temperature dependency of physiological rates studied by Taucher and Oschlies (2011).

Projected changes in the extension of oxygen minimum zones
Oxygen minimum zones (OMZs) are key oceanic regions because of their role in the marine nitrogen cycle (watercolumn denitrification occurs almost exclusively in O 2deficient waters) and because of the unusual ecosystems associated with low-O 2 regions (OMZs represent a respiratory barrier for many organisms).O 2 <5 mmol m −3 (as in Cocco et al., 2013).Only models that simulate volumes "close" ( * ) to the ones estimated with WOA-corrected oxygen climatology (Bianchi et al., 2012) are colored.* in between +100 and −50 % of the data-based estimated volume.
Following Cocco et al. (2013), we used three different thresholds (5, 50 and 80 mmol m −3 ) to characterize water volumes of low-oxygen waters.Suboxic waters are defined with a threshold of 5 mmol m −3 , whereas hypoxic waters are defined here with a threshold of 50 mmol m −3 .Figure 10 presents the relative evolution of these three volumes as simulated by the CMIP5 models over 1870 to 2100 and for the RCP8.5 scenario.Although each model is plotted for each of the three volumes in Fig. 10, we highlight the models for which the simulated volumes over 1990-1999 fall within +100 and −50 % of the observed volumes (126, 60 and 2.6 millions of km 3 , respectively), as estimated from the revised WOA 2005 database (Bianchi et al., 2012).
By 2100, all models project an increase in the volume of waters below 80 mmol m −3 , ranging from +1 % (MPI-ESMs) to +9 % (CESM1-BGC), relative to 1990-1999.This response is more consistent than that of the previous generation of Earth system models, i.e., changes varying from −26 to +16 % over 1870 to 2099 under the SRES-A2 scenario (Cocco et al., 2013).Conversely, for lower oxygen levels, there is less agreement among the CMIP5 models.Simulated volumes do not agree with observations, and changes due to climate change may be either negative or positive.For the volume of waters below 50 mmol m −3 , four models project an expansion of 2 to 16 % (both GFDL-ESMs, HadGEM2-ES and CESM1-BGC), whereas two other models project a slight contraction of ∼ 2 % (NorESM1-ME and MPI-ESM-MR).For the volume of waters below 5 mmol m −3 , only one model (IPSL-CM5A-MR) is close to the volume estimated from observations and simulates a large expansion of this volume (+30 % in the 2090s), These results for low-O 2 waters (5 and 50 mmol m −3 ) agree with those of Cocco et al. (2013), with large model-data and model-model discrepancies and simulated responses varying in sign for the evolution of these volumes under climate change.
The ability of climate models to represent O 2 observations has been questioned in recent studies.Stramma et al. (2012) used an ESM of intermediate complexity to perform simulations over the historical period and compared the simulated subsurface O 2 trends with observations.They showed that the model was unable to reproduce the spatial patterns of observed changes.Similarly, Andrews et al. (2013) compared output of MPI-ESM-LR and HadGEM2-ES to observations over the historical period.They reported that both models fail to reproduce the pattern of O 2 loss recorded by observations in low-latitude OMZs.
A more thorough analysis of the mechanisms responsible for the model-data discrepancies as well as the mechanisms driving the simulated future changes is necessary.Gnanadesikan et al. (2012) performed such an analysis with simulations carried out with a previous version of the GFDL Earth system model (GFDL-ESM2.1)under the SRES-A2 scenario.They show that the volume of suboxic waters does not increase under global warming in the tropical Pacific.A detailed analysis of the different terms contributing to the O 2 budget showed that an increase in O 2 in very low oxygen waters is associated with an enhanced supply of O 2 through lateral diffusion and increased ventilation along the Chilean coast.These results cast doubt on the ability of the present generation of models to project changes in O 2 accurately at the regional level, especially for low-O 2 waters, and stress the need for more model-data comparisons over the historical period alongside a better understanding of reasons for model biases.

A global-scale perspective
The existence of potential synergistic effects between the different stressors discussed here emphasizes the need to study them together (Boyd et al., 2008).Figure 11 shows the temporal model-mean evolution of global surface pH, global oxygen content and global NPP vs. global sea surface warming for each of the RCPs over 2006-2099.For RCP8.5, all these relationships appear linear, implying a constant fraction of acidification, deoxygenation, and NPP reduction per degree of warming: mean surface pH decreases by 0.127 pH unit • C −1 , global oxygen content decreases by 1.3 %/ • C and global NPP loses 3.3 % • C −1 (with R 2 of 0.99, 0.99 and 0.95, respectively).
For the other RCP scenarios, relationships are similar for surface pH vs. SST and for NPP vs. SST, but the relationship breaks down for oxygen content versus sea surface warming (Fig. 11).That is, deoxygenation continues long after sea surface temperatures have stabilized.In particular for the RCP2.6 scenario, the total content of oxygen in the ocean loses an additional 1 % in the second half of the 21st century after sea surface warming has been stabilized at +0.7 • C in 2050 (Fig. 3).When plotted against heat content change however, a single linear relationship for O 2 content emerges across the different scenarios with a slope of ∼ −0.149 %/10 22 J or 3.9 nmol J −1 (R 2 > 0.99) in the RCP8.5 scenario (Fig. 11).This slope of 3.9 nmol O 2 J −1 is slightly lower than that found in early deoxygenation studies (e.g., 6 nmol/J in Bopp et al., 2002), but more consistent with the recent study of Frölicher et al., (2009), thus indicating a larger contribution (here around 40 %) of warminginduced solubility reduction in global deoxygenation (whose contribution is estimated at 1.5 nmol J −1 , Bopp et al., 2002).

Analysis within a water-mass framework
For more insight into regional relationships between the different stressors, we computed trends of temperature, pH and O 2 in distinct water masses as described in the methodology section.Because coupled climate models have strong biases in the way they simulate the distribution of the main oceanic water masses (see Sallee et al., 2013, for an evaluation of Southern Ocean water masses in the CMIP5 models), this water mass framework is a natural approach to analyze model behavior (Iudicone et al., 2011) or to compare models.Furthermore, it avoids averaging biogeochemical properties between different water masses.
Figure 12 details the global relationships between temperature and pH and between temperature and O 2 , respectively, across the different RCP scenarios and for the four distinct water masses (WMs).
Not surprisingly, this analysis reveals large differences between the different WMs, and for instance between tropical waters (TWs) and mode/intermediate waters (MIWs).TWs are characterized by relatively low acidification-to-warming and low deoxygenation-to-warming ratios (Fig. 12).In contrast, MIWs are characterized by much higher acidificationand deoxygenation-to-warming ratios, demonstrating the importance of these water masses in the propagation of the acidification and deoxygenation signals to the ocean's interior.In addition, MIWs have a low buffer capacity of the carbonate system thus amplifying the response of pH to the uptake of anthropogenic carbon (see Sect. 3.2.3).The deep and bottom water masses (DWs and BWs) show very little acidification and warming trends, but some significant deoxygenation trend for the high-emission scenarios (RCP8.5).The climate change patterns differ when the model results are distinguished by an oceanic basin.Figure 13 shows the relationships between pH and temperature, and between O 2 and temperature for three of the four WMs (omitting TW), averaged over large oceanic basins (North Atlantic, North Pacific and Southern Ocean) for RCP8.5 (2090s minus 2000s).Remarkably, the responses of different stressors in the water masses of different basins are distinct and relatively robust across the range of models used here, demonstrating the promise of such an approach.The North Pacific and the Southern Ocean display similar behavior with similar ratios of acidification-to-warming and deoxygenation-to- warming for the different water masses with larger changes in both stressors for MIW and less in deep and bottom waters.The Southern Ocean shows, however, stronger signals than the North Pacific for MIW and BW, whereas the North Pacific shows stronger signals for DW (Fig. 13).WMs of the North Atlantic are characterized by much larger uncertainties (model spread) and by acidification-and deoxygenationto-warming ratios different from the ratios of other basins (Fig. 13).The North Atlantic MIW shows very small changes in temperature and oxygen, but a large and robust signal in pH (−0.26 pH unit).In the North Atlantic DW, the models project a very strong acidification and deoxygenation signals, with −0.16 pH unit and −13 mmol m −3 for a warming of only 0.1-0.2• C.
Although a more thorough analysis is needed to assess the mechanisms responsible for these changes and their relationship with ocean circulation and water mass formation, this first analysis illustrates the unique character of MIW, across the ocean, in terms of its intensity of deoxygenation and acidification.This study also reveals strong differences between the North Atlantic relative to the Southern Ocean and the North Pacific.The large model spread for warming and deoxygenation in the North Atlantic, for all water masses, is likely to be linked to differences between models for the weakening of the Atlantic Meridional Overturning Circulation (Cheng et al., 2013) and for changes in the distribution of North Atlantic WMs.

Conclusions
Using here 10 ESMs participating to the recent Coupled Model Intercomparison Project 5, we assess how the major stressors of marine ecosystems, namely ocean temperature, pH, dissolved oxygen concentrations and primary productiv- ity, may evolve over the course of the 21st century and under several atmospheric CO 2 pathways.All models project consistent warming, acidification, deoxygenation and lowering of primary productivity, whose intensities strongly depend on the scenario.Sea surface warming, sea surface pH reduction, decrease in global oxygen content and decrease in integrated primary productivity range from +0.71 to +2.73 • C, from −0.07 to −0.33 pH unit, from −1.81 to −3.45 %, and from −2.0 to −8.6 %, respectively, in the 2090s compared to the 1990s and for all RCP scenarios covered here.
Turning to the regional scale, the simulated evolutions of the different stressors differ (Fig. 5 and 6). Figure 14 summarizes the regional responses of the different stressors at the end of the 21st century under RCP8.5 scenario.To do so, we highlight the regions where each stressor change exceeds a threshold defined using the global-scale mean change.In addition, we also highlight the simulated presentday low-oxygen regions in which a decrease of dissolved O 2 would have much more impact on ecosystems than in an O 2enriched water mass.In brief, the tropical oceans are characterized by high warming rates, low acidification rates, inconsistent changes in subsurface oxygen, and high rates of NPP decrease.The North Pacific is characterized by high warming rates, high acidification rates, large decrease in subsurface oxygen and a mixed and inconsistent response in NPP.The North Atlantic is characterized by low warming rates, high acidification rates, medium-to-large decrease of subsurface oxygen and large decrease of NPP.Finally, the Southern Ocean is characterized by very low warming rates, high acidification rates, medium-to-large decrease of subsurface oxygen and a mixed response in NPP.
The mechanisms explaining these modifications are not analyzed in detail here, but consistency with previous model  intercomparison studies suggests continued robustness in the general drivers of these changes.Ocean acidification is driven by the atmospheric evolution of CO 2 with little feedback from climate change.Regional modulations of acidification rates are explained by differences in uptake rates of anthropogenic carbon and in carbonate buffering capacities of the different water masses.Deoxygenation is due to not only a warming-induced reduction in solubility but also to increased stratification and reduced ventilation.Reduction of NPP in the tropics and in the North Atlantic is also due to increased stratification and reduced ventilation, but increasing NPP at higher latitudes is explained by an alleviation of light limitation and/or temperature limitation.
Lack of robustness for these projections is especially critical for low-latitude subsurface O 2 concentrations (and OMZ evolution) and for NPP in some high-productivity regions such as the east equatorial Pacific.These high uncertainties in some of the projections urge the continuation of inter-model comparisons to understand these model differences.In addition, these uncertainties call for much more caution when using the CMIP5 models to force regional impact models, especially when relying on a single model and projection.
One critical improvement in future ESMs will be amelioration of biases in the representation of OMZs.A related long-term goal will be vastly enhanced resolution to represent the scales of coastal upwelling and other mesoscale phenomena such as eddies.Finally, representation of ecosystems in models such as these is an evolving science.They represent only a small set of the processes controlling the ecosystem and biogeochemical function.While the models are each constructed in mathematically defensible forms, they are all different in the underlying assumptions.Rather than representing discrete biological forms, they represent ecosystems as a biological continuum with infinite biodiversity in some ways (i.e., the role of temperature), and an artificial rigidity in others (i.e., fixed half-saturation constants).

Fig. 2 .
Fig. 2. Taylor diagrams showing the correspondence between annual-mean model results and observations for (a) SST, 3-D temperature, subsurface O 2 averaged over 200-600 m, 3-D dissolved O 2 , vertically integrated NPP; and (b) surface pH, 3-D pH, surface carbonate, and 3-D carbonate ion.Data-based products are from Reynolds et al. (2008) for SST, computed using DIC and alkalinity from Key et al. (2004) for pH and CO 2− 3 , from WOA 2009 for dissolved O 2 concentrations and 3-D temperature, from Behrenfeld et al. (1997) for NPP.For NPP, all coastal points (identified with a bathymetry of less than 500 m) are left out.All data are from 1990-1999 except observed NPP from 1997-2006, and O 2 concentration and 3-D temperature from WOA 2009 climatologies.The angular coordinate indicates the correlation coefficient (R), the radial coordinate shows the normalized standard deviation (stdmodel/stdobs).A model perfectly matching the observations would reside in point [1,1].

Fig. 3 .
Fig. 3. Model-mean time series of global sea surface warming ( • C), surface pH change (pH unit), ocean O 2 content change (%), and global NPP change (%) over 1870-2100 using historical simulations as well as all RCP simulations.Shading indicates one inter-model standard deviation.All variables are plotted relative to 1990-1999.

Fig. 5 .
Fig. 5. Change in stressor intensity (defined as the change in the magnitude of the considered variable) in 2090-2099 relative to 1990-1999 under RCP8.5.Multi-model mean of (a) sea surface warming ( • C), (b) surface pH change (pH unit), (c) subsurface dissolved O 2 concentration change (averaged between 200 and 600 m, mmol m −3 ), and (d) vertically integrated NPP (gC m −2 yr −1 ).Stippling marks high robustness.Robustness is estimated from inter-model standard deviation for SST and pH, from agreement on sign of changes for O 2 and NPP.Dark red color shading is used to mark the change in stressor that is detrimental for the marine environment.

Fig. 6 .
Fig. 6.Change in stressor intensity (defined as the change in the magnitude of the considered variable) in 2090-2099 relative to 1990-1999 under RCP2.6.Multi-model mean of (a) sea surface warming ( • C), (b) surface pH change (pH unit), (c) subsurface dissolved O 2 concentration change (averaged between 200 and 600 m, mmol m −3 ), and (d) vertically integrated NPP (gC m −2 yr −1 ).Stippling marks high robustness.Robustness is estimated from inter-model standard deviation for SST and pH, from agreement on sign of changes for O 2 and NPP.

Fig. 7 .
Fig. 7.Model-mean time series of mean surface carbonate ion concentrations (µmol L −1 ) for the global ocean, the Arctic Ocean (north of 70 • N), and the Southern Ocean (south of 60 • S) over 1870-2100 using historical simulations as well as all RCP simulations.Also indicated are estimates of aragonite and calcite saturation levels for the Arctic and the Southern oceans (dashed lines).

Fig. 8 .
Fig. 8. Change in pH in 2090-2099 relative to 1990-1999 under RCP8.5.(a) Model-mean depth of the maximum changes in pH.Robustness estimated from inter-model standard deviation and (b) model-mean changes of pH at that depth.

Fig. 9 .
Fig. 9. Change in export production of organic particles at 100 m under RCP8.5.(a) Multi-model mean change of export production in 2090-2999 relative to 1990-1999 (%).Stippling marks robustness from agreement on sign of changes.(b) Time series for all models of the change of export production relative to 1990-1999 (%) from 1870-2099 using historical and RCP8.5 simulations.

Fig. 11 .
Fig. 11.Relations between model-mean changes in surface pH (pH unit), global O 2 content (%) and global NPP and model-mean sea surface warming ( • C) for all scenarios.For global O 2 content, changes are also plotted against changes in total heat content (10 22 J).All changes are relative to1990-1999 and plotted over  2006-2100.

Fig. 12 .
Fig. 12. Relations between model-mean changes in pH (pH unit), dissolved O 2 (mmol m −3 ) and model-mean temperature change ( • C), in four distinct global water masses (tropical water mass, modal and intermediate water mass, deep water mass and bottom water mass) and for the four different RCPs.For a definition of the different water masses, please see Sect.2.3.Inter-model standard deviations are also indicated.

Fig. 13 .
Fig. 13.Relation between model-mean changes in pH (pH unit), dissolved O 2 (mmol m −3 ) and model-mean temperature change ( • C), in three distinct global water masses (modal and intermediate water mass, deep water mass and bottom water mass), for three different basins (North Atlantic, North Pacific and Southern Ocean) and for RCP8.5.For a definition of the different water masses, please see Sect.2.3.Inter-model standard deviations are also indicated.

Fig. 14 .
Fig. 14.Change in multiple stressor intensity (defined as the change in the magnitude of the considered variable) in 2090-2099 relative to 1990-1999 under the RCP8.5 scenario.The thresholds (T SST , T pH , T O 2 , T NPP ) for SST, surface pH, subsurface O 2 and integrated NPP have been defined as the sum of the global mean change in 2090-2999 and of the model-mean spatial standard deviation, such that ∼ 1/3 of the global ocean with the largest changes is depicted for each of the stressor.Red indicates where sea surface warming exceeds its threshold T SST (+3.64 • C), hatched green indicates where surface pH decreases by more than T pH (−0.36 pH unit), hatched yellow indicates where subsurface (200-600 m) oxygen concentrations decrease by more than T O 2 (−19.9 mmol m −3 ), and hatched blue indicates where vertically integrated NPP decreases by more than T NPP (−79.8 gC m −2 yr −1 ).When the simulated modelmean change is not robust (see text for the definition of the robustness index), no color is shown.In addition, hatched orange indicates present-day low-oxygen (< 50 mmol m −3 ) subsurface (200-600 m) waters.

Table 1a .
A brief description of the models used in this study, indicating atmospheric and oceanic resolution (lev is the number of levels on the vertical, and then the horizontal resolution is indicated in degrees), the marine biogeochemical component included in the ESM, and the different RCP scenarios performed and used in this study.

Table 3 .
Evolution of SST, pH, global O 2, integrated NPP for all models and all scenarios.Model means in2090-2099 (compared to  1990-1999)and inter-model standard deviations in 2090-2099.