Twenty-first century ocean warming, acidification, deoxygenation, and upper-ocean nutrient and primary production decline from CMIP6 model projections

Anthropogenic climate change is projected to lead to ocean warming, acidification, deoxygenation, reductions in near-surface nutrients, and changes to primary production, all of which are expected to affect marine ecosystems. Here we assess projections of these drivers of environmental change over the twenty-first century from Earth system models (ESMs) participating in the Coupled Model Intercomparison Project Phase 6 (CMIP6) that were forced under the CMIP6 Shared Socioeconomic Pathways (SSPs). Projections are compared to those from the previous generation (CMIP5) forced under the Representative Concentration Pathways (RCPs). A total of 10 CMIP5 and 13 CMIP6 models are used in the two multi-model ensembles. Under the high-emission scenario SSP5-8.5, the multi-model global mean change (2080–2099 mean values relative to 1870–1899) ± the inter-model SD in sea surface temperature, surface pH, subsurface (100–600 m) oxygen concentration, euphotic (0–100 m) nitrate concentration, and depth-integrated primary production is $+\mathrm{3.47}\pm \mathrm{0.78}$ ^∘C, $-\mathrm{0.44}\pm \mathrm{0.005}$, $-\mathrm{13.27}\pm \mathrm{5.28}$, $-\mathrm{1.06}\pm \mathrm{0.45}$ mmol m⁻³ and $-\mathrm{2.99}\pm \mathrm{9.11}$ %, respectively. Under the low-emission, high-mitigation scenario SSP1-2.6, the corresponding global changes are $+\mathrm{1.42}\pm \mathrm{0.32}$ ^∘C, $-\mathrm{0.16}\pm \mathrm{0.002}$, $-\mathrm{6.36}\pm \mathrm{2.92}$, $-\mathrm{0.52}\pm \mathrm{0.23}$ mmol m⁻³, and $-\mathrm{0.56}\pm \mathrm{4.12}$ %. Projected exposure of the marine ecosystem to these drivers of ocean change depends largely on the extent of future emissions, consistent with previous studies. The ESMs in CMIP6 generally project greater warming, acidification, deoxygenation, and nitrate reductions but lesser primary production declines than those from CMIP5 under comparable radiative forcing. The increased projected ocean warming results from a general increase in the climate sensitivity of CMIP6 models relative to those of CMIP5. This enhanced warming increases upper-ocean stratification in CMIP6 projections, which contributes to greater reductions in upper-ocean nitrate and subsurface oxygen ventilation. The greater surface acidification in CMIP6 is primarily a consequence of the SSPs having higher associated atmospheric CO₂ concentrations than their RCP analogues for the same radiative forcing. We find no consistent reduction in inter-model uncertainties, and even an increase in net primary production inter-model uncertainties in CMIP6, as compared to CMIP5.

1 Introduction

1.1 Ocean warming, acidification, deoxygenation, nutrient stress, and reduced primary production

Since the pre-industrial period the global oceans have experienced fundamental changes in physical and biogeochemical conditions as a result of anthropogenic climate change. Although these changes reflect the climate services that the oceans provide through heat and carbon storage, they also have major implications for the health of marine ecosystems. Ocean ecosystems are affected by the direct and indirect consequences of climate change. Atmospheric warming and rising CO₂ concentrations drives ocean warming and acidification, while these direct factors cause changes that modulate other important components of the ocean system, such as oxygenation, nutrient levels, and net primary production.

Temperature is a principal determinant of biological metabolism in the ocean (e.g. Eppley, 1972) and plays a major role in shaping the global distribution of marine species (e.g. Thomas et al., 2012; Sunagawa et al., 2015). The radiative forcing associated with greenhouse gas emissions results in an accumulation of heat in the Earth system, most of which is taken up by the oceans (Frölicher et al., 2014). Globally averaged sea surface temperature (SST) has increased by +0.7 ^∘C over the last 100 years (Bindoff et al., 2007), with observations indicating that the heat content trend in the upper 2000 m of the ocean has increased from 0.55 to 0.68 $\mathrm{J}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ since 1991 (Cheng et al., 2019).

Earth system models project twenty-first century increases in SST under all of the Representative Concentration Pathways (RCPs; Bopp et al., 2013). While certain marine organisms may have the potential to acclimate to rising ocean temperatures, poleward range shifts of many species have already been observed (Gregory et al., 2009; Sorte et al., 2010), with associated declines in tropical diversity projected (Thomas et al., 2012). Concurrently, the frequency, intensity, and duration of ocean heat waves has increased in the observational record and is projected to substantially increase in the future (Frölicher et al., 2018). This has already had serious impacts on marine foundation taxa such as corals, seagrasses, and kelps (Garrabou et al., 2009; Hobday et al., 2016; Smale et al., 2019).

A consequence of ocean warming is an increase in vertical density gradients and enhanced stratification. This results in a reduction in the supply of nutrients to the euphotic zone, with enhanced nutrient limitation generally leading to observed declines in net primary production (Behrenfeld et al., 2001; Behrenfeld et al., 2006). Earth system model projections consistently show enhanced stratification and associated reductions in euphotic zone nutrient concentrations under scenarios of climate change (Bopp et al., 2001; Sarmiento et al., 2004; Cabré et al., 2014; Fu et al., 2016). This generally results in projected global reductions in net primary production that are driven by enhanced phytoplankton nutrient limitation in the low-latitude oceans (Steinacher et al., 2010; Bopp et al., 2013; Krumhardt et al. 2017; Kwiatkowski et al., 2017; Moore et al., 2018). The projected magnitude of net primary production declines is highly uncertain across model ensembles (Bopp et al., 2013; Krumhardt et al., 2017), in part due to concurrent changes in phytoplankton light and temperature limitation, as well as altered top-down grazing, all of which can compensate for nutrient-driven production declines (Taucher and Oschlies, 2011; Laufkötter et al., 2015). However, declines in phytoplankton primary production are consistently amplified in higher trophic levels such as zooplankton (Chust et al., 2014; Stock et al., 2014; Kwiatkowski et al., 2018) and fish (Lotze et al., 2019).

Dissolved oxygen in the ocean exerts a strong control on marine ecosystems. At low O₂ levels, marine animals are unable to sustain aerobic metabolism, which can lead to mortality (Vaquer-Sunyer and Duarte, 2008). Oxygen levels also affect many oceanic biogeochemical cycles through an impact on redox reactions and microbial metabolism (e.g. on the nitrogen cycle, Gruber, 2004). Global warming is driving a global decline of dissolved oxygen in the ocean, referred to as ocean deoxygenation, because of a warming-induced reduction in O₂ solubility and increased stratification and reduced ventilation (Keeling et al., 2010; Oschlies et al., 2018). A recent assessment, based on three different analyses (Helm et al., 2011; Schmidtko et al., 2017; Ito et al., 2017), concluded that the oxygen content over the first 1000 m of the ocean has decreased by 0.5 to 3.3 % (or 0.7 to 3.5 % between 100 and 600 m) over 1970–2010 (Bindoff, et al., in press). In coastal systems, this warming effect is exacerbated by the effects of increased loading of nutrients and organic matter, which also lead to oxygen decline and an increase in coastal ocean dead zones (Breitburg et al., 2018). Earth System model projections consistently show continuing declines in oxygen over the twenty-first century as a function of the employed scenario (Bopp et al., 2013; Cocco et al., 2013), with large uncertainties in the tropics and for the evolution of oxygen minimum zones (Cabré et al., 2015).

The uptake of carbon by the oceans affects marine chemistry via ocean acidification (Gattuso and Buddemeier, 2000; Orr et al., 2005; Doney et al., 2009), a process that increases seawater concentrations of CO₂, H⁺, and ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and reduces pH and ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ ion concentrations. The oceans have absorbed approximately 30 % of anthropogenic carbon emissions since the pre-industrial era (Sabine et al., 2004; Khatiwala et al., 2009, 2013; Gruber et al., 2019), resulting in global surface pH declines of approximately 0.1 units (Bindoff et al., 2007). Declines in global open ocean surface pH are 0.018 units per decade over 1991–2011 (Lauvset et al. 2015), with individual time series stations exhibiting declines of 0.017 to 0.027 units per decade (Bindoff, et al., in press). Earth system models have projected twenty-first century global surface ocean pH declines of up to 0.33 units under previous high-emissions scenarios (Bopp et al., 2013), with associated changes in the seasonal cycles of seawater carbonate chemistry (McNeil and Sasse, 2016; Kwiatkowski and Orr, 2018; Landschützer et al., 2018).

The impact of ocean acidification on marine species is extensive and diverse. Calcifying species, such as echinoderms, bryozoans, and cnidarians, exhibit depressed calcification, growth, and survival under acidification (Kroeker et al., 2010; Albright et al., 2016; Kwiatkowski et al., 2016), altering the competitive balance in ecosystems (Kroeker et al., 2013). In teleost fish and marine invertebrates, ion exchange is reduced under acidification, depressing protein synthesis and metabolic rates (Langenbuch et al., 2006; Pörtner, 2008). Physiological and behavioural functioning is also sensitive to acidification, with olfactory discrimination (Munday et al., 2009) and predator–prey responses (Watson et al., 2014, 2017) shown to be impaired under more acidified conditions.

Marine organisms typically experience changes in multiple physical and geochemical conditions simultaneously, with impacts determined by the interactions between potential stressors. For example, the combined effect of warming and deoxygenation is projected to force poleward and vertical contractions of metabolically viable habitat for marine ectotherms (Deutsch et al., 2015). At the physiological level, experimental studies indicate that synergistic effects between potential marine stressors are common (Gunderson et al., 2016). Compound warming and acidification has been shown to exacerbate negative impacts on photosynthesis, calcification, reproduction, and survival of marine organisms (Harvey et al., 2013), while compound exposure to acidification and low oxygen can also have synergistic effects (McBryan et al., 2013) and may reduce the thermal tolerance of certain species (Pörtner, 2010).

Here we assess future projections of climate-related drivers of marine impacts within the Coupled Model Intercomparison Project Phase 6 (CMIP6; Eyring et al., 2016; O'Neill et al., 2016) simulations, evaluating how these differ from previous CMIP5 (Taylor et al., 2011) simulations. We focus on projected changes in ocean temperature, pH, dissolved O₂ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration, and net primary production across 13 CMIP6 and 10 CMIP5 Earth system models. Expanding on previous CMIP intercomparison studies (e.g. Steinacher et al., 2010; Bopp et al., 2013; Cocco et al., 2013), which typically concentrate on upper-ocean mean state changes in biogeochemistry, we also assess benthic impacts and changes in seasonal cycles.

1.2 Ocean biogeochemical model development since CMIP5

A comprehensive assessment of changes between CMIP5 and CMIP6 in the ocean biogeochemical components of Earth system models (ESMs) and their associated skill is provided in Séférian et al. (2020). Since CMIP5, CMIP6 has seen a general increase in the horizontal grid resolution of physical ocean models and a limited increase in vertical resolution. The latter may be particularly important for ecosystem projections as it directly affects simulated stratification, a key factor influencing changes in ocean impact drivers (Capotondi et al., 2012; Bopp et al., 2013; Laufkötter et al., 2015; Kwiatkowski et al., 2017) and their impact on higher trophic levels (Stock et al., 2014; Chust et al., 2014; Kwiatkowski et al., 2018; Lotze et al., 2019). Updates in the representation of ocean biogeochemical processes between CMIP5 and CMIP6 have generally included increases in model complexity (Séférian et al., 2020). Specifically, CMIP6 models provide more widespread inclusion of dissolved oxygen, micronutrients, such as iron, variable stoichiometric ratios, and improved representation of lower trophic levels including heterotrophic bacteria and the cycling and sinking of organic matter (Séférian et al., 2020).

Relative to CMIP5, the CMIP6 Earth system models display an improved ability to reproduce the modern mean-state distribution of a number of key biogeochemical tracers (Séférian et al., 2020). Although global-scale totals of ocean carbon flux and net primary production estimates have been improved in CMIP6 with respect to CMIP5, the simulated geographical distribution of present-day mean state air–sea CO₂ fluxes and surface chlorophyll concentrations show only moderate improvements between CMIP5 and CMIP6. There are also moderate improvements in the representation of subsurface dissolved oxygen concentrations in most ocean basins. Model skill in the representation of surface macronutrient concentrations in CMIP6 has improved for dissolved silicate but declined slightly for nitrate.

2 Methodology

The analysis of projected multiple ocean impact drivers presented here focuses on three key depth levels: the upper ocean, the thermocline, and the benthic zone. The surface zone is where most biological activity is concentrated in the oceans and where impacts from climate change are typically the greatest. Specifically, we assess projections of surface ocean temperature, surface ocean pH, subsurface dissolved O₂ concentration (averaged between 100 and 600 m), upper-ocean ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration (averaged between 0 and 100 m) and net primary production (depth integrated over the full water column). The choice of vertical range for O₂ reflects the potential importance of the expansion of oxygen minimum zones, which are more prominent at such depths. The choice of vertical range for ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reflects its importance as a critical macronutrient supporting primary production in the euphotic zone. Both vertical ranges are chosen to be compatible with the recent assessment of marine drivers in the IPCC Special Report on the Ocean and Cryosphere (Bindoff, et al., in press). Additionally, for the CMIP6 models we assess benthic ecosystem drivers, focussing on projections of bottom temperature, pH, and O₂ concentration. The benthic level is defined as the bottom-ocean model layer at each grid point. As such, its exact depth depends on vertical discretisation and bathymetry, which differs across the CMIP6 ensemble. All benthic model outputs were corrected for potential drift at the grid cell level (e.g. Gehlen et al., 2014; Séférian et al., 2016) using coincident pre-industrial control simulations.

2.1 Processing and analysis of model outputs

All ESMs assessed in the CMIP5 and CMIP6 ensembles (Tables 1 and 2) include physical ocean models and coupled ocean biogeochemistry schemes that account for some or all of the potential ocean impact drivers: temperature, pH, O₂, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, and net primary production. A total of 10 CMIP5 and 13 CMIP6 models are assessed with the model ensemble size differing among scenarios depending on contributions from each model group. The CMIP5 ensemble is the same as that used in the comprehensive assessment of projected ocean drivers provided by Bopp et al. (2013). Only one ensemble member per model is used for a given scenario. That is, in CMIP terminology we typically use ensemble member “r1i1p1” from each CMIP5 model and “r1i1p1fx” from each CMIP6 model (where “fx” is the recommended set of external forcings employed by the various modelling groups). Consequently, we do not assess the role of internal variability in the emergence of climate-related changes in marine ecosystems drivers (e.g. Frölicher et al., 2016; Lovenduski et al., 2016; Krumhardt et al. 2017; Freeman et al., 2018). Two of the CMIP6 models included in our analysis (GFDL-CM4 and ACCESS-ESM1.5) do not include ${\mathrm{NO}}_{\mathrm{3}}^{-}$ as a prognostic tracer. Hence, their ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations were calculated from modelled total dissolved inorganic phosphorus assuming a constant Redfield ratio of 16:1.

Table 1 The CMIP6 Earth system models used in this study; their individual components used to represent ocean, sea ice, and marine biogeochemistry; and the ocean impact drivers and simulations that were assessed.

Table 2 The CMIP5 Earth system models used in this study; their individual components used to represent ocean, sea ice, and marine biogeochemistry; and the simulations that were assessed. All models provided temperature, pH, oxygen, nitrate, and NPP outputs, as in Bopp et al. (2013).

To facilitate intercomparison, model output on each native grid was regridded to the same regular $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ horizontal grid using distance weighted average remapping (climate data operators; remapdis). Model outputs were kept on their native vertical grids, with vertical discretisation ranging from 40 (MPI-ESM1.2) to 75 (IPSL-CM6A-LR, CNRM-ESM2-1 and UKESM1-0-LL) levels, except for models using hybrid or isopycnic vertical coordinates for which model outputs were vertically regridded to 35 (GFDL-ESM4, GFDL-CM4) and 70 (NorESM2-LM) levels. Following generally adopted practice (e.g. Bopp et al., 2013), all models were given equal weighting in the respective CMIP6 and CMIP5 ensemble mean. However, within the CMIP6 ensemble two modelling groups contributed two ESMs and within the CMIP5 ensemble three modelling groups contributed two ESMs, which is likely to influence the extent of model independence (Masson and Knutti, 2011; Knutti et al., 2015; Sanderson et al., 2015; Lovenduski and Bonan, 2017).

The CMIP5 historical simulations had variable start dates between 1850 and 1861, all of which finished in 2005; the subsequent RCP simulations started in 2006 and were run until at least 2099. In CMIP6, there is greater temporal consistency. All CMIP6 historical simulations were made over 1850–2014, while the subsequent SSP scenarios started in 2015 and ran until at least 2100. To facilitate comparison between CMIP5 and CMIP6, the historical and future projections of ocean impact drivers in both phases of CMIP are presented as anomalies relative to 1870–1899 mean values of their respective historical simulations. When solely evaluating twenty-first century projections in the SSPs, however, the last 20 years of the CMIP6 historical simulations (1995–2014) are used as a baseline period. Throughout the analysis, the uncertainty associated with global mean projections is assessed using the inter-model SD (given ± uncertainties). At regional scales, projection robustness is evaluated using previously adopted approaches (e.g. Bopp et al., 2013), including whether the magnitude of the multi-model mean anomaly exceeds the inter-model SD or if there is at least 80 % model sign agreement. The interquartile range of regional projections is given in the annexes.

2.2 From Representative Concentration Pathways to Shared Socioeconomic Pathways

Aside from changes in ESMs, a fundamental difference between CMIP5 and CMIP6 is that they differ in the future scenarios used for anthropogenic emissions and land use change. Those scenarios are derived from integrated assessment models and based on plausible future pathways of societal development. In CMIP6, the Shared Socioeconomic Pathways (SSPs) provided via the Scenario Model Intercomparison Project (ScenarioMIP) are used instead of the RCPs that were used in CMIP5 (O'Neill et al., 2016). The SSPs provide revised emission and land use scenarios relative to the RCPs (Riahi et al., 2017).

In this study, we confine our assessment of ocean impact drivers to concentration-driven simulations, focussing on Tier 1 SSPs of ScenarioMIP (O'Neill et al., 2016), that is SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5, which result in end-of-century approximate total radiative forcing levels of 2.6, 4.5, 7.0, and 8.5 W m⁻², respectively. The SSPs have generally higher associated concentrations of atmospheric CO₂ and lower associated atmospheric concentrations of CH₄ and N₂O relative to their RCP counterparts (O'Neill et al., 2016; Meinshausen et al., 2019). This is particularly the case for SSP5-8.5, which in comparison to RCP8.5 assumes that coal constitutes a greater proportion of the primary energy mix in the second half of the twenty-first century (Kriegler et al., 2017). Given that differences among projections of surface ocean acidification are dominated by scenario uncertainty, with relatively little inter-model uncertainty and internal variability (e.g. Bopp et al., 2013; Frölicher et al., 2016), such changes in atmospheric concentrations of CO₂ are expected to have a large impact on projections of ocean pH and related carbonate system variables.

Alongside the assessment of the SSP concentration-driven model outputs, outputs from models forced under the four RCPs (RCP2.6, RCP4.5, RCP6.0 and RCP8.5) are also assessed in parallel. This allows some comparison with past CMIP5 assessments (e.g. Bopp et al., 2013). However, RCP6.0 has no direct SSP analogue, while SSP3-7.0 has no direct RCP analogue.

3 Results and discussion

3.1 Comparison with historical global trends

Observed historical trends in global mean surface ocean temperature, pH, and subsurface oxygen were compared with the multi-model mean of the CMIP6 ensemble over the corresponding years of historical simulations (Table 3). Global observations of historical trends in euphotic-zone nitrate concentrations and integrated primary production were deemed insufficiently robust, given the associated interannual–decadal variability, to be assessed in the models (Elsworth et al. 2020). The observed 1901–2012 SST warming of +0.06 ^∘C per decade is well reproduced in the CMIP6 ensemble, although SST warming between 1979 and 2012 is warm-biased in the multi-model mean. Global surface ocean acidification of −0.018 units per decade between 1991 and 2011 is also well reproduced by the CMIP6 models, particularly considering that the observed trend is a reconstruction based on discrete surface ocean fCO₂ measurements and alkalinity estimates that are derived from temperature and salinity (Lauvset et al., 2015). Finally, with respect to subsurface deoxygenation, the observed dissolved oxygen trend of −0.30 to −1.52 mmol m⁻³ per decade from 1970 to 2010 (90 % confidence range; Bindoff, et al., in press) encompasses the CMIP6 multi-model mean response over the corresponding years. Given the performance of the CMIP6 models at reproducing ocean biogeochemical mean conditions (Séférian et al., 2020) and trends, they are deemed appropriate to project future trends in biogeochemistry under the SSPs.

Table 3 Historical trends in ocean impact drivers in observations and CMIP6 models. Observed multi-decadal historical trends in global mean sea surface temperature, surface pH, and subsurface O₂ (averaged between 100 and 600 m) and the corresponding trend in the CMIP6 ensemble. Uncertainty estimates are the inter-model SD.

3.2 Global upper-ocean projections

Under all SSPs, global multi-model mean sea surface temperature is projected to increase, while surface pH, subsurface dissolved oxygen concentration, euphotic-zone nitrate concentration and net primary production are projected to decline during the twenty-first century (Fig. 1). The projected change in the five ocean impact drivers increases with associated radiative forcing across the four SSPs. Under the high-mitigation SSP1-2.6 scenario, the end-of-century model mean changes (2080–2099 mean values relative to 1870–1899) in sea surface temperature, surface pH, subsurface oxygen concentration, euphotic nitrate concentration, and net primary production are $+\mathrm{1.42}\pm \mathrm{0.32}$ ^∘C, $-\mathrm{0.16}\pm \mathrm{0.002}$, $-\mathrm{6.36}\pm \mathrm{2.92}$, $-\mathrm{0.52}\pm \mathrm{0.23}$ mmol m⁻³, and $-\mathrm{0.56}\pm \mathrm{4.12}$ %, respectively. Under the high-emissions scenario SSP5-8.5, the corresponding changes are 3.47± 0.78 ^∘C, $-\mathrm{0.44}\pm \mathrm{0.005}$, $-\mathrm{13.27}\pm \mathrm{5.28}$, $-\mathrm{1.06}\pm \mathrm{0.45}$ mmol m⁻³, and $-\mathrm{2.99}\pm \mathrm{9.11}$ % (Table 4), respectively. Across these two scenarios, the separation between CMIP6 projections of sea surface temperature and pH and to a lesser extent oxygen and nitrate, further demonstrate the effectiveness of intense mitigation strategies in limiting twenty-first century marine ecosystem exposure to potential stress. This is in agreement with assessments of previous multi-model projections (e.g. CMIP5; Bopp et al., 2013).

Table 4 Global mean changes in multiple ocean impact drivers across the CMIP6 and CMIP5 ensembles. Global mean anomalies of sea surface temperature, surface ocean pH, subsurface dissolved O₂ concentration (averaged between 100 and 600 m), upper-ocean ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (averaged between 0–100 m), and depth-integrated net primary production for the CMIP6 SSPs and CMIP5 RCPs. Anomalies are given as 2080–2099 mean values relative to the 1870–1899 mean. Uncertainty estimates are the inter-model standard deviation.

Figure 1 Global mean projections of upper-ocean impact drivers. Global mean projections of (a) sea surface temperature (^∘C), (b) surface ocean pH, (c) subsurface dissolved O₂ concentration (averaged between 100–600 $\mathrm{m};\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}$), (d) euphotic-zone ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (averaged between 0–100 $\mathrm{m};\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}$), and (e) depth-integrated net primary production ( %). Values are anomalies relative to the 1870–1899 reference period. CMIP6 mean anomalies for the historical and SSP simulations are shown as solid lines with shading representing the inter-model SD. CMIP5 projections only show the multi-model mean. The model ensemble size for each scenario is given in parentheses.

Following previous assessments (Bopp et al., 2013), inter-model uncertainty is estimated as the inter-model SD. Although some of this model spread is due to internal variability, this contribution is limited for global averages and its relative contribution to inter-model uncertainty is expected to decline throughout the twenty-first century (Frölicher et al., 2016). Relative to scenario uncertainty, which is estimated as the maximum difference between mean SSP projections, inter-model uncertainty is extremely low for surface pH projections, which show distinct separation between the SSPs prior to 2050. The low inter-model uncertainty associated with projections of surface ocean pH is well characterised and associated with the identical CO₂ forcing used by all ESMs in concentration-driven SSP and RCP projections (Lovenduski et al., 2016), a weak climate–pH feedback (Orr et al., 2005; McNeil and Matear, 2007), limited interannual variability, and consistently adopted standards for ESM ocean carbonate chemistry equations (Orr et al., 2017). Surface ocean pCO₂ and corresponding carbonate chemistry generally follow changes in atmospheric CO₂ with a global mean equilibration time of approximately 8 months (Gattuso and Hansson, 2011). The differences between projected surface pH across the SSPs therefore reflect the divergence of prescribed atmospheric CO₂ concentrations, i.e. the different scenarios.

In comparison to pH, projections of SST exhibit greater inter-model uncertainty (Fig. 1). This uncertainty is likely to result from differences in climate sensitivity between models. Historically, such differences have been attributed to diversity in cloud feedbacks and to a lesser extent water vapour and lapse-rate feedbacks (Andrews et al., 2012; Vial et al., 2013). For projections of subsurface oxygen and euphotic-zone nitrate concentrations, inter-model uncertainty is greater still and can exceed scenario uncertainty. This greater inter-model uncertainty is a result of oxygen and nitrate concentrations being strongly influenced by both physical changes (e.g. changes in solubility, circulation, and mixing) and changes in biological sources and sinks (Stramma et al., 2012; Fu et al., 2016; Bopp et al., 2017; Oschlies et al., 2018).

The inter-model uncertainty associated with CMIP6 net primary production projections is consistently larger than the scenario uncertainty. Indeed, for each SSP, individual models project both increases and decreases in global primary production, with the inter-model SD encompassing positive and negative anomalies (Fig. 1). This is a consequence of net primary production changes reflecting a diverse and delicately balanced suite of bottom-up and top-down ecological processes, which are variously parameterised across models. Bottom-up changes in phytoplankton growth rates are typically driven by a combination of enhanced nutrient limitation and reduced temperature and light limitation (Doney, 2006; Steinacher et al., 2010; Taucher and Oschlies, 2011), while top-down changes in zooplankton grazing rates can simultaneously influence the stock of phytoplankton biomass (Laufkötter et al., 2015). The accurate simulation of many of the biogeochemical tracers upon which net primary production (NPP) depends (e.g. the distribution of iron; Tagliabue et al., 2016) represents a significant and ongoing challenge to ESMs (Séférian et al., 2020).

3.3 Regional patterns of upper-ocean change

Global-scale projections of end-of-century upper-ocean impact drivers (2080–2099 anomalies relative to 1995–2014 mean values) exhibit spatial variability that is both ocean impact driver and SSP dependent (Fig. 2). CMIP6 projections of SST show near global increases under both SSP1-2.6 and SSP5-8.5, with generally high regional robustness across the model ensemble. The greatest warming is evident in the Northern Hemisphere, particularly the Arctic Ocean and high-latitude North Pacific, where multi-model mean warming can exceed 2 ^∘C in SSP1-2.6 and 5 ^∘C in SSP5-8.5. This Arctic amplification is well established in both observations (Bekryaev et al., 2010) and models, and thought to be primarily driven by temperature and surface albedo feedbacks (Screen and Simmonds, 2010; Pithan and Mauritsen, 2014). The notable exception to warming is in the subpolar North Atlantic where there is minor cooling in SSP1-2.6 and limited warming in SSP5-8.5. Although it is associated with high inter-model uncertainty (Figs. 2b,c and A1), this “warming hole” is well documented in both observations and models and typically related to a slow down in the Atlantic meridional overturning circulation (Drijfhout et al., 2012; Menary and Wood, 2018). Spatial patterns of SST anomalies are broadly consistent with those of the CMIP5 ensemble (Bopp et al., 2013).

Figure 2 Projections of multiple upper-ocean impact drivers under SSP1-2.6 and SSP5-8.5. CMIP6 mean historical climatologies and twenty-first century anomalies in (a–c) sea surface temperature (^∘C), (d–f) surface ocean pH, (g–i) subsurface dissolved O₂ concentration (averaged between 100 and 600 $\mathrm{m};\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}$), (j–l) euphotic-zone ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (averaged between 0 and 100 $\mathrm{m};\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}$), and (m–o) depth-integrated net primary production ($\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{y}}^{-\mathrm{1}}$). Anomalies are 2080–2099 mean values relative to the 1995–2014 baseline period. Stippling designates areas of projection robustness. For temperature and pH this is defined as the magnitude of the mean anomaly exceeding the inter-model SD. For O₂, NO₃, and NPP this is defined as at least 80 % model sign agreement.

Anomalies in surface ocean pH are ubiquitously negative under both SSP1-2.6 and SSP5-8.5, with very low associated inter-model uncertainty (Figs. 2e and f and A1). Consistent with past model projections (McNeil and Matear, 2007; Steinacher et al., 2009; Bopp et al., 2013), the greatest declines in pH are projected in the higher latitudes and the Arctic Ocean in particular, where model mean declines can exceed 0.45 in SSP5-8.5 (2080–2099 anomalies relative to 1995–2014; Fig. 2e, f). The enhanced Arctic Ocean acidification reflects both the large surface ocean warming in the Arctic, which acts to decrease pH, as well as the related loss of permanent and semi-permanent sea ice (McNeil and Matear, 2007; Steinacher et al., 2009). Sea ice melt increases anthropogenic carbon uptake and decreases pH by both providing a greater surface area for air–sea gas exchange and simultaneously enhancing air–sea CO₂ fluxes by dilution of dissolved inorganic carbon concentrations with freshwater (Yamamoto-Kawai et al., 2009; Yamamoto et al., 2012).

Although global mean subsurface (100–600 m) O₂ concentration is projected to decline under all SSPs, there is a high degree of variability in projections at regional scales (Fig. 2h,i). The largest declines in subsurface O₂ generally occur at higher latitudes and in particular in the North Pacific, where declines in the multi-model mean can exceed 45 mmol m⁻³ in SSP5-8.5. In equatorial regions of the Atlantic and Indian Ocean and upwelling regions of the Pacific, increases in subsurface O₂ concentration are projected under both SSP1-2.6 and SSP5-8.5. However, these increases have high inter-model uncertainty and are at odds with historical observations of OMZ expansion (Stramma et al., 2008; Helm et al., 2011; Schmidtko et al., 2017; Ito et al., 2017), with an assessed range of 3.0 to 8.3 % between 1970 and 2010 (Bindoff et al., in press), though it has been suggested that such observations are the result of climate variability (Deutsch et al., 2011; Bindoff et al., in press). These subsurface O₂ increases are, however, consistent with previous projections, including those from CMIP5, which have highlighted that coarse-resolution models struggle to reproduce subsurface ventilation pathways in these regions (Stramma et al., 2012; Andrews et al., 2013; Busecke et al., 2019).

Figure 3 Change in subsurface oxygen saturation and apparent oxygen utilisation under SSP5-8.5. CMIP6 multi-model mean changes in (a) subsurface dissolved O₂ concentration (averaged between 100 and 600 $\mathrm{m};\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}$), (b) subsurface O₂ saturation (O_2 sat), and (c) subsurface apparent oxygen utilisation (AOU) in 2080–2099 of SSP5-8.5 relative to 1995–2014. Stippling designates robustness, as defined by at least 80 % model sign agreement (only a subset of CMIP6 models were considered).

For a subset of the CMIP6 models (CanESM5, CNRM-ESM2-1, GFDL-CM4, IPSL-CM6A-LR, MIROC-ES2L, MPI-ESM1.2-HR and UKESM1-0-LL), projected changes in subsurface O₂ concentration under SSP5-8.5 were decomposed into changes in O₂ saturation (O_2 sat) and apparent oxygen utilisation (AOU), where $\mathrm{\Delta }{\mathrm{O}}_{\mathrm{2}}=\mathrm{\Delta }{\mathrm{O}}_{\text{2\hspace{0.17em}sat}}-\mathrm{\Delta }\text{AOU}$ (Fig. 3). O_2 sat was computed from model temperature and salinity outputs and represents the effect of oxygen solubility changes on dissolved O₂ concentration, while AOU was calculated as $\mathrm{\Delta }\text{AOU}=\mathrm{\Delta }{\mathrm{O}}_{\text{2\hspace{0.17em}sat}}-\mathrm{\Delta }{\mathrm{O}}_{\mathrm{2}}$, and is affected by both changes in biological consumption of O₂ and in ventilation and stratification. The heightened reductions in subsurface O₂ in the North Pacific and North Atlantic are shown to be the result of consistent reductions in O_2 sat and increases in AOU, which act to reinforce O₂ concentration declines. In contrast, the projected increases in O₂ in the tropical Indian and Atlantic Oceans are shown to be the result of reductions in AOU that lower oxygen demand more than the concurrent reductions in O_2 sat. These tropical reductions in AOU are generally robust across the model ensemble despite this not being the case for the coincident increases in O₂. The spatial patterns of CMIP6 projected changes in subsurface O_2 sat and AOU under SSP5-8.5 are similar to that of the CMIP5 models under RCP8.5 (Bopp et al., 2017). The general reduction in O_2 sat has been shown to be predominantly due to warming-driven reductions in solubility, while the heightened AOU declines in the North Pacific and North Atlantic have been primarily attributed to reductions in ventilation and an increase in the age of these waters (Bopp et al., 2017; Tjiputra et al., 2018).

CMIP6 multi-model mean projections of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations in the euphotic zone (0–100 m) show variable regional declines under SSP1-2.6 and SSP5-8.5 (Fig. 2k,l). These projected declines are robust and largest in the Arctic Ocean, equatorial eastern Pacific, North Atlantic, and North Pacific, where they can exceed 3 mmol m⁻³ in SSP5-8.5. ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations show limited anomalies in the subtropical gyres where concentrations are already very low. The CMIP6 spatial pattern of euphotic-zone ${\mathrm{NO}}_{\mathrm{3}}^{-}$ anomalies is in broad agreement with CMIP5 projections (Fu et al., 2016).

Projections of primary production anomalies are highly diverse across regions (Fig. 2n,o). The global CMIP6 multi-model mean decline in primary production is shown to be primarily driven by declines in the North Atlantic and the western equatorial Pacific, which can exceed 40 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{y}}^{-\mathrm{1}}$ under SSP5-8.5. In the high latitudes, primary production generally increases, with anomalies approaching 20 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{y}}^{-\mathrm{1}}$ in parts of the Arctic and Southern Oceans under SSP5-8.5. Such changes have historically been associated with enhanced stratification as the upper ocean warms (Doney, 2006). In tropical and mid-latitude waters, where phytoplankton are nutrient limited, this tends to reduce the vertical nutrient supply and exacerbate nutrient stress. In contrast, in high-latitude waters, where phytoplankton are typically light limited, enhanced stratification can reduce mixing below the euphotic depth and therefore result in reduced light stress. However, the simplicity of this paradigm has been challenged by observations, which show regionally variable coupling between changes in stratification and primary production on interannual timescales (Lozier et al., 2011; Dave and Lozier, 2013), with recent studies demonstrating the additional importance of changes to the horizontal advection of nutrients (Whitt, 2019) and zooplankton grazing (Laufkötter et al., 2015) in shaping regional primary production responses.

Although the general pattern of NPP changes is similar in CMIP6 compared to CMIP5, regional declines are reduced in magnitude, less spatially extensive, and are typically less robust. This is particularly notable in the Indian Ocean and subtropical North Pacific, which were regions of consistent NPP decline in CMIP5 projections (Bopp et al., 2013) but exhibit limited robust trends in CMIP6. The projected increase in NPP in the high latitudes is, however, broadly consistent with previous model intercomparisons (Steinacher et al., 2010; Bopp et al., 2013), which have attributed these increases to the retreat of sea ice, reduced deep mixing and corresponding reductions in light limitation.

3.4 Stratification and mixed-layer depth

The difference between densities at the surface, 200, and 1000 m are used as stratification indices in the CMIP6 projections and given alongside changes in the maximum mixed-layer depth (Fig. 4). Global multi-model mean stratification is consistently projected to increase over the twenty-first century, with greater increases under SSPs that have higher radiative forcing. Under SSP1-2.6, the global change in stratification index is $+\mathrm{0.13}\pm \mathrm{0.05}$ kg m⁻³ between the surface and 200 m and $+\mathrm{0.26}\pm \mathrm{0.08}$ kg m⁻³ between the surface and 1000 m, while under SSP5-8.5 this increase to $+\mathrm{0.58}\pm \mathrm{0.11}$ and $+\mathrm{0.90}\pm \mathrm{0.20}$ kg m⁻³, respectively. Over the same period, the global mean maximum annual mixed-layer depth shoals by 7.0±3.3 m in SSP1-2.6 and 19.5±2.6 m in SSP5-8.5. With the exception of the Arctic Ocean, multi-model projections of increasing stratification are typically consistent across most regions, with robustness increasing under SSP5-8.5 and when the upper 1000 m of the water column is considered. Projected shoaling of the maximum mixed-layer depth is also generally robust across the multi-model ensemble, albeit with slightly less model consistency, as would be expected given the coincident climatic changes in wind stress and surface heat fluxes (Fig. 4h,i).

Figure 4 Change in stratification and mixed-layer depth in SSP1-2.6 and SSP5-8.5. CMIP6 mean historical climatologies and anomalies in (a–c) stratification index between 200 m and the surface (kg m⁻³), (d–f) stratification index between 1000 m and the surface (kg m⁻³), and (g–i) maximum annual mixed-layer depth (m). Anomalies are 2080–2099 mean values relative to 1995–2014. The stratification index is defined as the difference in density between given depths and the surface. Stippling designates robustness, as defined by the magnitude of the mean anomaly exceeding the inter-model SD.

The regions of projected primary production decline in the North Atlantic and western equatorial Pacific are typically associated with heightened increases in stratification, notably in the upper 200 m of the water column. In these regions, particularly the North Atlantic, the maximum mixed-layer depth also shoals. As such, there is strong evidence that reduced vertical mixing and entrainment of nutrients into the upper ocean is, at least partially, responsible for these regional declines in primary production. However, similar increases in stratification and reductions in mixed-layer depth occur in regions such as the North Pacific and Indian Ocean, where declines in primary production are largely absent. Therefore, further assessment of simultaneous changes in processes such as nutrient advection (e.g. Whitt, 2019), nitrogen fixation (Riche and Christian, 2018), the microbial loop (e.g. Schmittner et al., 2008; Taucher and Oschlies, 2011), and top-down grazing pressure (e.g. Laufkötter et al., 2015) are required to fully understand the regional primary production response in CMIP6.

3.5 Compound stressors

The projected occurrence of multiple potential ecosystem stressors in the upper ocean was determined across the SSPs using prescribed thresholds of surface warming ($>+\mathrm{2}$ ^∘C), surface acidification ($<-\mathrm{0.2}$ units), subsurface deoxygenation ($<-\mathrm{30}$ mmol m⁻³) and euphotic-zone ${\mathrm{NO}}_{\mathrm{3}}^{-}$ decline ($<-\mathrm{1}$ mmol m⁻³), with anomalies calculated as 2080–2099 mean values relative to 1995–2014 (Fig. 5). It should be noted that our choice of stressor thresholds, based on the magnitude of biogeochemical anomalies, are somewhat arbitrary. Indeed, it could be argued that absolute biogeochemical thresholds, for example as defined for hypoxia or oligotrophy, may better reflect potential ecosystem stress. Moreover, the thresholds take no account of regional differences in natural temperature, pH, O₂, and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ variability, which may mediate ecosystem responses to changes in mean conditions (e.g. Kroeker et al., 2020). That being said, a single threshold that encompasses the variety of ecosystem responses to a particular stressor likely does not exist.

Figure 5 Compound upper-ocean impact drivers. Regions where projected CMIP6 sea surface warming exceeds 2 ^∘C, euphotic-zone (0–100 m) ${\mathrm{NO}}_{\mathrm{3}}^{-}$ decline exceeds 1 mmol m⁻³, surface ocean pH decline exceeds 0.2, and subsurface (100–600 m) dissolved O₂ concentration decline exceeds 30 mmol m⁻³ in (a) SSP1-2.6, (b) SSP2-4.5, (c) SSP3-7.0, and (d) SSP5-8.5. The exceedance of driver thresholds is determined from 2080 to 2099 anomalies relative to 1995–2014 values.

The concurrent exceedance of multiple thresholds increases with associated radiative forcing across the SSPs, indicative of greater potential compound ecosystem stressors. The tropical and subtropical oceans are generally characterised by projected compound warming and acidification under SSP3-7.0 and SSP5-8.5, with additional nutrient thresholds exceeded in regions of equatorial upwelling. The North Pacific is characterised by high sensitivity to potential compound stressors, with all thresholds of warming, acidification, deoxygenation, and nutrient decline exceeded under SSP5-8.5. In contrast, the projected occurrence of compound stressors is limited in the Southern Ocean, where only the acidification threshold is consistently exceeded. The North Atlantic is characterised by sensitivity to combined acidification and nutrient stress, while the Arctic Ocean is sensitive to compound warming, acidification, and nutrient stress.

3.6 CMIP6 vs. CMIP5 projections

While the temporal behaviour of changes in ocean impact drivers is similar across the CMIP5 and CMIP6 model ensembles (Fig. 1), the CMIP6 Earth system models generally project greater global surface ocean warming, acidification, subsurface deoxygenation, and euphotic zone NO₃ reductions than the CMIP5 projections performed with comparable radiative forcing (Fig. 6, Table 4). The CMIP6 models, however, project reduced global primary production declines relative to comparable CMIP5 simulations. There is no consistent reduction in inter-model uncertainty in CMIP6. In fact, with respect to projections of primary production, inter-model uncertainty is substantially increased in CMIP6.

Figure 6 Comparison between CMIP6 and CMIP5 end-of-century changes in upper-ocean impact drivers. (a) Atmospheric CO₂ concentration and radiative forcing derived from the MAGICC6 model in year 2100 for the CMIP6 SSPs and CMIP5 RCPs. (b–e) Global mean anomalies of (b) surface ocean pH, (c) subsurface dissolved O₂ concentration (averaged between 100 and 600 $\mathrm{m};\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}$), (d) ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration (averaged between 0 and 100 $\mathrm{m};\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}$), and (e) integrated net primary production ( %) against anomalies of sea surface temperature (^∘C). Anomalies are 2080–2099 mean values relative to the 1870–1899 baseline period. Error bars represent the inter-model SD.

The projected end-of-century SST increase (2080–2099 minus 1870–1899) in SSP1-2.6, SSP2-4.5, and SSP5-8.5 is higher than in RCP2.6, RCP4.5 and RCP8.5, respectively. This enhanced CMIP6 warming is attributable to generally greater climate sensitivity in the CMIP6 model ensemble relative to the CMIP5 ensemble (Forster et al., 2019). Indeed, the same version of a reduced complexity climate model (MAGICC7.0) run with CMIP5 and CMIP6 forcings, projects marginally greater warming of near-surface air temperatures in the RCPs than comparative SSPs (Meinshausen et al., 2019), underlining that the greater SST increases in CMIP6 are likely driven by changes to models and not forcing datasets.

The enhanced acidification in CMIP6 relative to CMIP5 is consistent across models (Fig. 6, Table 4) and attributable to higher prescribed atmospheric CO₂ levels in the forcing of the SSP scenarios relative to the RCP scenarios with equivalent radiative forcing (Meinshausen et al., 2019). Year 2100 atmospheric CO₂ levels are 1135.2, 602.8, and 445.6 ppm in SSP5-8.5, SSP2-4.5, and SSP1-2.6, respectively. The corresponding levels in RCP8.5, RCP4.5, and RCP2.6 are 936, 538, and 421 ppm, respectively (Meinshausen et al., 2011). Therefore, although the SSP and RCP simulation pairs have analogous end-of-century radiative forcing, the higher CO₂ levels in the SSPs result in greater acidification for the CMIP6 projections.

The greater euphotic-zone ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration declines in SSPs compared to their RCP analogues are likely a consequence of the enhanced surface warming in CMIP6 models. This warming results in a greater increase in upper-ocean stratification than that projected in CMIP5 models (Cabré et al., 2014; Fu et al., 2016). At the global scale, models have been shown to consistently exhibit strong negative correlations between relative stratification anomalies and relative ${\mathrm{NO}}_{\mathrm{3}}^{-}$ anomalies on interannual timescales (Fu et al., 2016). The greater increases in stratification in CMIP6 therefore result in greater reductions in mixing and entrainment of nutrient-rich deep waters into the euphotic zone in comparison with CMIP5.

The enhanced subsurface deoxygenation in the SSPs relative to comparable RCPs is likely the consequence of both physical and biogeochemical processes (e.g. Bopp et al., 2017; Oschlies et al., 2018). The greater warming in CMIP6 projections results in a greater reduction in O₂ solubility, while also affecting the ventilation and transport of O₂ within the ocean interior. In addition, concurrent changes in biological production, export and respiration can either mitigate or exacerbate physically driven subsurface deoxygenation (Oschlies et al., 2018).

The reduced declines in global net primary production projected for the twenty-first century in the SSPs, relative to comparable RCPs, combined with large increases in the associated inter-model uncertainty, is striking (Fig. 6e), particularly given that declines in euphotic zone ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations are typically greater in the SSPs. This suggests that the temporal evolution of phytoplankton resource limitation and grazing pressure under climate change may have significantly altered between CMIP5 and CMIP6. In previous CMIP biogeochemistry intercomparisons, all models projected global primary production declines, albeit with large inter-model uncertainty (Steinacher et al., 2010; Bopp et al., 2013). However, in CMIP6, four of the models (CESM2, CESM2-WACCM, CNRM-ESM2-1, and IPSL-CM6A-LR) consistently project global increases in primary production across the SSPs and are primarily responsible for both the reduced multi-model mean declines and the large increase in inter-model SD. Global increases in net primary production have been previously documented in other model studies (Schmittner et al., 2008; Taucher and Oschlies, 2011; Laufkötter et al., 2015) and attributed to temperature-dependent intensification of the microbial loop increasing regenerated production. Further analysis of the CMIP6 models that project primary production increases is clearly required to determine whether this is also the case or additional processes (e.g. the temporal evolution of nitrogen fixation or iron limitation) explain why they differ from previous generations of the same Earth system model family.

3.7 Global benthic ocean projections

On average, bottom waters are consistently projected to warm, acidify, and deoxygenate across the twenty-first century (Fig. 7). Under SSP1-2.6, the end-of-century model mean changes (2080–2099 relative to 1870–1899) in bottom-water temperature, pH, and dissolved O₂ are $+\mathrm{0.12}\pm \mathrm{0.03}$ ^∘C, $-\mathrm{0.018}\pm \mathrm{0.001}$, and $-\mathrm{5.14}\pm \mathrm{2.04}$ mmol m⁻³, respectively. Under SSP5-8.5 the corresponding changes are $+\mathrm{0.22}\pm \mathrm{0.04}$ ^∘C, $-\mathrm{0.030}\pm \mathrm{0.002}$, and $-\mathrm{6.04}\pm \mathrm{2.19}$ mmol m⁻³ (Table 5). Thus, even for bottom waters, CMIP6 projections highlight that intense mitigation strategies can limit ecosystem exposure to potential warming and acidification stress during the twenty-first century (e.g. Tittensor et al., 2010; Levin and Le Bris, 2015).

Table 5 Global mean projected changes in benthic ocean impact drivers in CMIP6. Global mean anomalies of bottom-water temperature (^∘C), pH, and dissolved O₂ concentration for the CMIP6 SSPs. Anomalies are 2080–2099 mean values relative to the 1870–1899 baseline period. Uncertainty estimates are the inter-model SD.

Figure 7 Global mean projections of benthic ocean impact drivers. CMIP6 global mean projections of benthic (a) temperature (^∘C), (b) pH, and (c) dissolved O₂ concentration (mmol m⁻³). Values are anomalies relative to the 1870–1899 reference period. Mean anomalies for the historical and SSP simulations are shown as solid lines, with shading representing the inter-model SD. The model ensemble size for each scenario is given in parentheses.

The magnitude of projected changes in bottom waters is less than in surface and upper-ocean waters, while bottom-water uncertainties for a given scenario are larger (Fig. 7). This contrast is particularly evident for pH projections with the SSPs, whose ranges of uncertainty fully separate before 2050 in the surface ocean (Fig. 1) but still have some overlap in 2080 for bottom waters. This relative increase in inter-model uncertainty results from surface ocean chemistry being in equilibrium with the same atmospheric CO₂ concentrations for all models. Conversely, benthic pH changes are strongly influenced by ocean circulation, which transports anthropogenic carbon from the upper ocean to the seafloor and is variably impacted by climate change across models (e.g. Gregory et al., 2005; Cheng et al., 2013). The increased uncertainty in pH projections with depth has been previously noted for CMIP5 projections in the North Atlantic (Gehlen et al., 2014) and Arctic Ocean (Steiner et al., 2014). For projected global deoxygenation in bottom waters, inter-model uncertainty is substantially larger than scenario uncertainty in CMIP6. As with projections of subsurface dissolved O₂, this larger model uncertainty results from the isolation of bottom waters from the atmosphere. Thus, bottom waters at a given temperature and salinity may deviate substantially from the value that would be determined by their solubility and air–sea equilibrium due to effects from other physical and biogeochemical processes (e.g. Oschlies et al., 2018).

3.8 Regional patterns of benthic ocean change

In bottom waters, the end-of-century spatial distributions of changes in temperature, pH, and dissolved O₂ are similar between SSPs (Fig. 8) and in broad agreement with CMIP5 projections (Sweetman et al., 2017). The intensity of warming, acidification, and deoxygenation is generally greater in SSP5-8.5 than SSP1-2.6 in benthic waters above 2000 m, while at greater depths the impact is similar for both SSPs.

Figure 8 Projections of multiple benthic ocean impact drivers under SSP1-2.6 and SSP5-8.5. CMIP6 mean historical climatologies and anomalies in benthic (a–d) temperature (^∘C), (e–h) pH, and (j–l) dissolved O₂ concentration (mmol m⁻³). Anomalies are 2080–2099 mean values relative to 1995–2014. Stippling designates areas of projection robustness. For temperature and pH this is defined as the magnitude of the mean anomaly exceeding the inter-model SD. For O₂ this is defined as at least 80 % model sign agreement. Vertical profiles show global mean benthic anomalies in 200 m depth intervals, with error bars denoting the inter-model SD.

The largest projected benthic warming in SSP1-2.6 and SSP5-8.5 occurs in continental shelf waters, the Arctic seas, and the Southern Ocean, where temperature increases can exceed 0.5 ^∘C by the end-of-century (2080–2099 average relative to the 1995–2014 baseline). In contrast, for most of the abyssal benthic ocean projected increases in temperature are less than 0.2 ^∘C. The characteristic North Atlantic “warming hole” present in projections of the surface ocean (Fig. 2) is also evident in benthic layers above 1000 m, such as the mid-Atlantic ridge (Fig. 8b,c). This represents the only major region of multi-model mean benthic cooling across SSP1-2.6 and SSP5-8.5 with high associated uncertainty. As in the surface ocean, this cooling is likely associated with a slow down of the Atlantic meridional overturning circulation (Drijfhout et al., 2012; Menary and Wood, 2018).

Projected end-of-century acidification is highly limited in most bottom waters. However, in the North Atlantic, Arctic seas, and certain continental shelf waters, pH changes can exceed −0.1 in SSP1-2.6 and −0.2 in SSP5-8.5. For shelf waters, the greater bottom-water pH declines can be the result of coupling between surface waters, which experience large changes in carbonate chemistry, and bottom waters (e.g. through mixing and entrainment), as well as benthic remineralisation of organic matter (Bates et al., 2009). In contrast, enhanced bottom-water acidification in the North Atlantic is associated with deep-water formation and high uptake of anthropogenic carbon (Sabine et al., 2004), which rapidly propagates anomalies in surface ocean chemistry to depth. Bottom-water acidification has been previously projected in the North Atlantic by an ensemble of CMIP5 models under RCP8.5 (Gehlen et al., 2014; Sweetman et al., 2017).

In contrast to temperature and pH, projections of benthic dissolved O₂ concentration show changes that are not predominantly confined to shelf waters and specific regions. Most of the global benthic ocean is projected to experience deoxygenation under both SSP1-2.6 and SSP5-8.5, even at depths below 2000 m (Fig. 8). Bottom-water deoxygenation is higher in the Southern Ocean, equatorial Pacific, and North Atlantic, where declines in the multi-model mean can exceed 20 mmol m⁻³.

It should be noted that Earth system models are not explicitly designed to explore the benthic biogeochemical response to climate change and that certain caveats should be considered. Model spin-up simulations, although longer in CMIP6 than CMIP5 (Séférian et al., 2020), are typically insufficient in length to equilibrate biogeochemical conditions in the deep ocean (Séférian et al., 2016), and therefore contemporaneous pre-industrial control simulations are required to correct biogeochemical drift. Vertical thickness of bottom-ocean layers is also highly variable across the CMIP6 ensemble, although for a given model resolution it is typically highest near the surface and decreases dramatically with depth. As such, the extent to which continental shelves are resolved greatly differs and uncertainties associated with resolution are pronounced in the abyssal ocean. Moreover, the representation of biogeochemical processes associated with ocean sediments and benthic ecosystems is typically absent or highly limited in ESMs.

3.9 Depth of maximum acidification

The depth of maximum end-of-century pH and [H⁺] change is often below the surface, and it varies regionally in CMIP6 projections (Fig. 9). Although the maximum pH change is usually found in surface waters in the high latitudes and upwelling regions, it is typically located between 200–400 m in subtropical mode and intermediate waters. Because of its log scale, if the change in pH were identical in surface and subsurface waters it would imply a larger absolute change in [H⁺] in the subsurface, where the mean [H⁺] is higher. Indeed, a change in pH represents a relative change in [H⁺], not an absolute change in that quantity. That relationship, combined with higher [H⁺] at depth, means that the maximum change in [H⁺] is usually deeper than it is for pH. Furthermore, the spatial distribution of the maximum change in pH and [H⁺] also differs.

Figure 9 Magnitude and depth of maximum pH and [H⁺] change under SSP5-8.5. The CMIP6 ensemble mean maximum change in (a) pH and (b) [H⁺] in 2080–2099 of SSP5-8.5 relative to 1995–2014. The mean depth at which the maximum (c) pH and (d) [H⁺] change is projected. Stippling designates robustness, as defined by the mean anomaly exceeding the inter-model SD.

Figure 10 Change in the seasonal amplitude of surface ocean [H⁺] and pH. The CMIP6 historical climatologies and SSP5-8.5 anomalies in peak-to-peak seasonal amplitude of surface ocean (a, b) [H⁺] and (c, d) pH. Anomalies are calculated from the mean seasonal amplitude in 2080–2099 relative to that in 1995–2014. Stippling designates robustness, as defined by at least 80 % model sign agreement.

Enhanced acidification in subsurface mode and intermediate waters has been observed at time series stations (Dore et al., 2009; Byrne et al., 2010; Bates et al., 2012) and in CMIP5 model projections (Resplandy et al., 2013; Bopp et al., 2013; Watanabe and Kawamiya, 2017). Although observational studies have suggested that this enhancement results from changes in circulation and biological activity (Dore et al., 2009; Byrne et al., 2010), model results indicate that it can be explained by the geochemical effect of rising atmospheric CO₂ and the particular carbonate chemistry of these waters (Orr, 2011; Resplandy et al., 2013). Specifically, the enhanced acidification sensitivity in mode and intermediate waters has been attributed to their lower temperatures and their higher ratio of dissolved inorganic carbon to total alkalinity relative to that found in surface waters of the same regions (Orr, 2011; Resplandy et al., 2013).

Figure 11 The seasonal amplitude of sea surface temperature. The CMIP6 multi-model mean peak-to-peak seasonal amplitude of sea surface temperature (^∘C) in (a) 1995–2014 of the historical simulations, (b) 2080–2099 of the SSP5-8.5 simulations, and (c) the change in seasonal amplitude between the two periods. Stippling designates robustness, as defined by at least 80 % model sign agreement.

3.10 Surface ocean seasonality

Changes in the seasonal amplitude of surface ocean temperature, pH, and hydrogen ion concentration ([H⁺]) were determined after detrending by subtracting a cubic spline fit from the monthly time series in each grid cell and then calculating the annual peak-to-peak amplitude for each year of the detrended dataset, following the approach of Kwiatkowski and Orr (2018). Under SSP5-8.5, the seasonal amplitude of global surface ocean [H⁺] is projected to increase by +73± 12 % across the CMIP6 ensemble (2080–2099 average relative to 1995–2014; Fig. 10). Concurrently, the seasonal amplitude of global surface ocean pH is projected to decrease by $-\mathrm{10}\pm \mathrm{5}$ %. Increases in the seasonal amplitude of [H⁺] are robust across all regions but are generally highest in the high latitudes. In contrast, declines in the seasonal amplitude of pH are typically only robust in the low latitudes and mid-latitudes, with inconsistent projections of pH seasonal amplitude projected in the Arctic and Southern Oceans.

The simultaneous amplification of [H⁺] and attenuation of pH seasonal cycles is consistent with previous assessments of CMIP5 projections, with Kwiatkowski and Orr (2018) showing [H⁺] seasonal amplification of +81± 16 % and pH seasonal attenuation of $-\mathrm{16}\pm \mathrm{7}$ % under RCP8.5 (2090–2099 anomalies relative to 1990–1999). Although counterintuitive, this results from the log scale of pH, which means that the seasonal amplitude of pH depends not only on the seasonal amplitude of [H⁺] but also on the inverse of the annual mean [H⁺]. As the projected increase in annual mean [H⁺] is usually greater than the corresponding increase in the seasonal amplitude of [H⁺], the seasonal amplitude of pH declines as a result. Increases in the seasonal cycle of [H⁺] have been shown to be primarily driven by the geochemical effect of increasing atmospheric CO₂. This affects both the seasonal amplitude of the controlling variables dissolved inorganic carbon and alkalinity, as well the sensitivity of [H⁺] to seasonal changes in temperature, dissolved inorganic carbon, and alkalinity (Kwiatkowski and Orr, 2018). Given the near-linear relationship between [H⁺] and pCO₂ on annual timescales (Orr, 2011), projected increases in the seasonal amplitude of [H⁺] are in agreement with historical observations (Landschützer et al., 2018) and twenty-first century projections from CMIP5 models (McNeil and Sasse, 2016; Gallego et al., 2018) of increasing pCO₂ seasonal amplitude.

The multi-model mean seasonal amplitude of global surface ocean temperature is projected to increase by $+\mathrm{0.59}\pm \mathrm{0.21}$ ^∘C across SSP5-8.5 (Fig. 11). Over most of the ocean, the seasonal amplitude of sea surface temperature is projected to show limited but robust increases ($<+\mathrm{0.5}$ ^∘C). However, in the North Atlantic, North Pacific, and Southern Ocean, increases in the seasonal amplitude of SST can exceed +2 ^∘C, and in the Arctic Ocean, the seasonal amplitude typically increases by $>+\mathrm{5}$ ^∘C.

The CMIP6 projections of the changing seasonal amplitude of SST under SSP5-8.5 are consistent with previous projections from the CMIP5 models (Carton et al., 2015; Alexander et al., 2018). The limited increases in SST seasonal amplitude for most of the global ocean have been attributed to greater relative shoaling of the mixed-layer depth in summer than in winter (Alexander et al., 2018). However, in the Arctic Ocean the large increase in SST seasonal amplitude is primarily due to the loss of sea ice. The seasonal melting and refreezing of sea ice accounts for approximately half of the present-day seasonal Arctic Ocean net surface heat flux, buffering seasonal variability in Arctic Ocean heat content and SSTs (Serreze et al., 2007; Fig. 11). The loss of this seasonal melting–freezing cycle under high-emissions scenarios such as RCP8.5 has been shown to account for a doubling of seasonal Arctic Ocean heat content variability. Ice loss further amplifies the seasonal cycle of SSTs by increasing the seasonal cycle of net surface heat fluxes. The net downward radiative flux increases in summer as albedo declines, while the net upward radiative flux increases in winter due to greater evaporative and sensible heat loss (Carton et al., 2015).

4 Conclusions

The latest CMIP6 Earth system models consistently project global surface ocean warming and acidification, subsurface deoxygenation, and euphotic-zone nitrate reductions in the twenty-first century. Multi-model mean projections of global net primary production show declines in the twenty-first century, although this is with large inter-model uncertainty. The projected change in these ocean impact drivers is shown to increase with radiative forcing across the SSPs, highlighting the benefit of emissions reductions to upper-ocean ecosystems. The magnitude of projected warming, acidification, and deoxygenation is lower in the benthic ocean, with greater inter-model uncertainty relative to scenario uncertainty. However, the extent of warming and acidification is still limited under lower-emissions scenarios, demonstrating the potential benefits of mitigation to benthic ecosystems.

In addition to changing mean-state conditions, the CMIP6 models also project changes to the seasonal cycles of temperature and carbonate chemistry under the SSPs. The seasonal amplitude of surface ocean acidity ([H⁺]) nearly doubles over the twenty-first century under SSP5-8.5, with a concurrent reduction in the seasonal amplitude of pH. Over the same period, the seasonal amplitude of temperature is projected to increase, particularly in the Arctic Ocean.

The CMIP6 projections of warming, acidification, deoxygenation, and nutrient reduction are greater than those of previous CMIP5 models under comparable radiative forcing. The enhanced acidification is a consequence of higher atmospheric CO₂ concentrations in the SSPs than their RCP analogues. The enhanced warming, however, reflects the greater climate sensitivity of the CMIP6 models. This increased warming results in greater increases in upper-ocean stratification, which contributes to greater reductions in euphotic nitrate and subsurface oxygen concentration. The CMIP6 multi-model mean projections of primary production declines are less than those of previous CMIP5 models under comparable radiative forcing; however, there is a large increase in inter-model uncertainty that requires further assessment.

Projected changes to the mean state and seasonality of biogeochemical ocean conditions are likely to present major challenges to diverse marine ecosystems from the surface ocean to abyssal depths. Potential organism stress is likely to be exacerbated by simultaneous exposure to multiple biogeochemical changes, emphasising the need for extensive emissions reductions.

Appendix A

Figure A1 The interquartile range associated with CMIP6 projections of upper-ocean impact drivers. The CMIP6 interquartile range (Q₃–Q₁) of projected anomalies in (a, b) sea surface temperature (^∘C), (c, d) surface ocean pH, (e, f) subsurface dissolved O₂ concentration (averaged between 100 and 600 $\mathrm{m};\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}$), (g, h) euphotic-zone ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (averaged between 0 and 100 $\mathrm{m};\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}$), and (i, j) depth-integrated net primary production ($\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{y}}^{-\mathrm{1}}$). The corresponding mean projections are given in Fig. 2.

Figure A2 The interquartile range associated with CMIP6 projections of stratification and mixed-layer depth. The CMIP6 interquartile range (Q₃–Q₁) of projected anomalies in (a, b) the stratification index between 200 m and the surface (kg m⁻³), (c, d) the stratification index between 1000 m and the surface (kg m⁻³), and (e, f) the maximum annual mixed-layer depth (m). The corresponding mean projections are given in Fig. 4.

Figure A3 The interquartile range associated with CMIP6 projections of benthic impact drivers. The CMIP6 interquartile range (Q₃–Q₁) of projected anomalies in benthic (a, b) temperature (^∘C), (c, d) pH, and (e, f) dissolved O₂ concentration (mmol m⁻³). The corresponding mean projections are given in Fig. 8.