Multi-centennial ocean biogeochemical responses to extended Shared Socioeconomic Pathways

Despite their profound consequences for the global carbon cycle, multi-centennial ocean biogeochemical responses to anthropogenic greenhouse gas forcing remain poorly constrained. Here we investigate long-term oceanic responses to atmospheric CO₂ forcing using the Model for Interdisciplinary Research on Climate, Earth System version 2 for Long-term simulations (MIROC-ES2L) Earth system model driven by extended Shared Socio-economic Pathway (SSP) scenarios through to the year 2500. Ocean–atmosphere pCO₂ disequilibrium exhibits strong scenario dependence. Under the low-emission scenario SSP1–2.6, oceanic CO₂ undersaturation gradually weakens and approaches near equilibrium by 2100. By contrast, under the high-emission scenario SSP5–8.5, rapid increase in atmospheric pCO₂ results in increasing undersaturation of CO₂ in the ocean over the 21st century. Thereafter, surface-ocean CO₂ supersaturation progressively expand into the equatorial and subtropical oceans between the mid–23rd and late 25th century. The Southern Ocean remains persistently undersaturated throughout the simulations, consistent with sustained influence from deep waters that are weakly affected by anthropogenic perturbations. Under SSP5–8.5, oceanic CO₂ uptake continues to increase until the late 21st century. Thereafter, it gradually weakens as surface-ocean pCO₂ approaches atmospheric levels, whereas the Southern Ocean remains a persistent carbon sink throughout the simulation. Across all scenarios, early oceanic pCO₂ changes are primarily driven by increases in dissolved inorganic carbon, whereas alkalinity becomes an increasingly important control from the mid–22nd century onward under SSP5–8.5. A progressive decline in the oceanic CO₂ buffering capacity, driven by cumulative CO₂ uptake, increases the sensitivity of surface-ocean pCO₂ to additional carbon uptake. The buffering capacity approaches its effective minimum levels by the late 22nd century. In parallel, enhanced upper-ocean stratification and a slowdown of the deep overturning circulation reduce the upward transport of high alkalinity waters to the surface. This results in a decline in surface-ocean alkalinity that further enhances the influence of alkalinity on surface-ocean pCO₂ during this period. Marine biological productivity and export vary significantly as a function of atmospheric CO₂ trajectory. It increases under low-emission scenarios but declines under sustained high-emission scenarios. The transition between the two behaviors occurs when $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ reaches approximately 3 times the preindustrial level and is closely associated with the weakening of the Atlantic Meridional Overturning Circulation and the resulting nutrient redistribution. Thermal stress and nutrient limitation persist in regions such as the Arctic Ocean until the late 25th century even under SSP1–2.6, indicating a centennial-scale lagged response of marine ecosystems to atmospheric CO₂ forcing. In the surface ocean, marine ecosystem stress (integrating thermal and biogeochemical stressors) continues to intensify through the late 22nd century under SSP5–8.5, despite the stabilization of atmospheric pCO₂, highlighting the long memory and limited reversibility of oceanic ecosystem stress.

1 Introduction

The ocean is a major sink for anthropogenic carbon dioxide (CO₂), currently absorbing about 25 % of annual global emissions; thus it plays a critical role in mitigating ongoing global warming (Friedlingstein et al., 2025). Because carbon–climate feedbacks strongly influence the climate system (Friedlingstein et al., 2006), reliable climate projections require an accurate understanding of how ocean carbon uptake responds to global warming.

Substantial progress has been made in understanding physical and biogeochemical responses of the ocean to climate change up to 2100, particularly through the Coupled Model Intercomparison Project (CMIP). According to CMIP5 multi-model means, marine ecosystems would face concurrent stressors – higher sea surface temperature (SST; +2.7 °C), lower surface pH (−0.33), reduced dissolved oxygen (−3.45 %), and declining net primary production (NPP; −8.6 %) – by the 2090s under Representative Concentration Pathway (RCP) 8.5 (Bopp et al., 2013).

Relative to CMIP5, CMIP6 models generally exhibit higher climate sensitivity (Forster et al., 2019). As a result, the magnitudes of warming, acidification, and nutrient and oxygen declines are higher in CMIP6 than in CMIP5 (Kwiatkowski et al., 2020). In CMIP6, however, despite stronger physical–chemical forcing, projected changes in NPP are smaller in magnitude than those in CMIP5, and are highly uncertain in sign. The variations among models are partly due to differences in the parameterizations of key processes such as nitrogen fixation and associated nutrient limitations; these factors alter the response of NPP even under the same scenario (Bopp et al., 2022). The increase in NPP is also driven by variability in the C:N uptake ratio of phytoplankton (Kwon et al., 2022).

Low-latitude productivity is crucial for the global carbon cycle. However, robust projection of this key variable remains a major challenge because of substantial inter-model spread (Rodgers et al., 2024). A major source of uncertainty lies in the ongoing debate over whether remote nutrient supply or local recycling processes are the dominant mechanism that sustains low-latitude productivity. Moore et al. (2018) argued that a reduction of nutrient supply from the Southern Ocean can strongly suppress low-latitude productivity. By contrast, Rodgers et al. (2024) highlighted that temperature-dependent remineralization enhances local mesopelagic nutrient recycling, together with water mass ventilation, thereby sustaining productivity under high-emission scenarios on centennial timescales. On longer timescales, however, Lee et al. (2025) showed that nutrient depletion within local mesopelagic layers can ultimately drive substantial declines in low-latitude production by 2500. Taken together, current evidence indicates that the drivers and timescales of biogeochemical responses remain uncertain and strongly model-dependent.

Recognizing this uncertainty, studies have emphasized the importance of regional processes. For example, Nakamura and Oka (2019) demonstrated that dominant controls of phytoplankton growth vary by region, underscoring the need for spatially explicit analyses. This spatial heterogeneity also characterizes air–sea CO₂ fluxes and ocean acidification: low- and high-latitude regions respond to distinct drivers. In the Southern Ocean, intensifying westerlies and freshening alter circulation in ways that both enhance surface CO₂ uptake and promote deep-ocean carbon storage (Ito et al., 2015). In addition, future ocean carbon uptake is projected to rely on a hybrid of solubility-driven winter uptake and biologically mediated summer uptake (Mongwe et al., 2024). By contrast, the Arctic Ocean is experiencing rapid acidification. Sea-ice retreat enhances CO₂ uptake (Bates et al., 2006), with models projecting a reversal of seasonal cycles of oceanic partial pressure of CO₂ ($p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$) and intensified summer acidification (Orr et al., 2022). Overall, polar oceans are shaped by freshwater input and sea-ice dynamics, whereas low-latitude responses are governed primarily by stratification and nutrient supply.

Analyses using CMIP models have also advanced our understanding of present-day oceanic uptake of anthropogenic CO₂ and its controls. Tjiputra et al. (2014) identified rising dissolved inorganic carbon (DIC) as the dominant contributor to higher surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ over the 21st century. This is primarily due to anthropogenic CO₂ uptake, whereas SST increase plays a secondary role. However, projections beyond 2100 remain scarce. For example, RCP8.5 simulations of the Community Earth System Model version 1 (CESM1) indicate that global oceanic CO₂ uptake peaks around the 2080s but remains substantial through 2300 (Chikamoto and DiNezio, 2021). More recent work shows that enhanced upper-ocean stratification reduces the upward transport of high alkalinity waters. This leads to declining surface-ocean alkalinity and subsequent reduction in carbonate-system buffering capacity, thereby slowing CO₂ uptake in later centuries (Chikamoto et al., 2023). Furthermore, while shoaling of CaCO₃ dissolution depths weakens the carbonate pump (Planchat et al., 2024).

Given that the deep ocean carbon reservoir is much larger than the atmospheric reservoir of carbon, the long-term adjustment of the ocean carbon reservoir is central to climate trajectories. Most studies, however, do not go beyond 2100 because multi-centennial simulations with comprehensive Earth system models (ESMs) are computationally expensive. Therefore, the mechanisms shaping centennial-to-millennial responses (2100–2500) remain poorly constrained. Although Earth system models of intermediate complexity (EMICs) have enabled simulations over these long timescales (Claussen et al., 2002; Weber, 2010), the multi-centennial mechanisms of ocean carbon uptake have been much less explored using comprehensive ESMs because of their substantially greater computational cost. Recent work has demonstrated that the mechanisms governing the climate–carbon cycle feedback differ between low- and high-emission scenarios on multi-centennial timescales, highlighting the need for long-term mechanistic analyses of marine carbon uptake (Khatiwala et al., 2026). Nevertheless, few studies have focused on such analyses using comprehensive ESMs.

Improving the representation of ocean biogeochemical processes in ESMs is important not only for constraining future carbon-cycle feedbacks but also for increasing confidence in projections of marine ecosystem responses. This is because ESMs provide the physical and biogeochemical forcing for a wide range of marine ecosystem studies (e.g., Masson-Delmotte et al., 2021). For example, the Fisheries and Marine Ecosystem Model Intercomparison Project (FishMIP) is based on CMIP6 ESM outputs. Its results indicate that substantial declines in global marine animal biomass and major regional shifts in marine ecosystems are likely to occur in response to climate change (Tittensor et al., 2021). These applications of ESM output demonstrate that uncertainties in simulated ocean biogeochemistry propagate into downstream impact assessments, influencing estimates of ecosystem change and fisheries impacts.

This study aims to contribute towards improving the representation of ocean biogeochemical processes in ESMs. We analyze extended simulations from the Model for Interdisciplinary Research on Climate, Earth System version 2 for Long-term simulations (MIROC-ES2L) to investigate the centennial-to-millennial evolution of ocean biogeochemistry and the oceanic CO₂ uptake under contrasting atmospheric CO₂ trajectories by extending three baseline Shared Socio-economic Pathways (SSP1–2.6, SSP2–4.5, and SSP5–8.5) through the year 2500.

We address three main questions: (1) How does ocean carbon uptake evolve after atmospheric pCO₂ ($p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$) stabilizes? (2) What processes govern the long-term evolution of ocean carbonate chemistry and buffering capacity? (3) How do biological productivity and ecosystem stress evolve on multi-centennial timescales? To answer these questions, we contrast SSP1–2.6, in which $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ peaks and declines early, with SSP5–8.5, in which $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ continues to rise until 2250.

We further examine the redistribution of biogeochemical tracers across depth horizons (epipelagic, mesopelagic, and deep ocean) and climate zones (tropical, subtropical, subpolar, and polar) to distinguish rapid responses that emerge during the 21st century from slower adjustment processes that continue for centuries thereafter. Using a single comprehensive ESM and multiple emission pathways allows us to investigate the relationships between the trajectory of anthropogenic forcing and long-term adjustments of ocean carbon storage, carbonate chemistry, and marine biogeochemistry.

2 Methods

2.1 Model description of the MIROC-ES2L and its ocean biogeochemical component

In this study, we use the MIROC-ES2L, which is an Earth system model that participated in CMIP6 (Hajima et al., 2020a). The physical core of MIROC-ES2L is MIROC 5.2, which interactively couples the Center for Climate System Research–National Institute for Environmental Studies atmospheric general circulation model (CCSR-NIES AGCM) (Numaguti et al., 1997), an online aerosol model (Spectral Radiation-Transport Model for Aerosol Species – SPRINTARS) (Takemura et al., 2000; Takemura, 2005), a coupled sea-ice–ocean general circulation model (CCSR Ocean Component – COCO) (Hasumi, 2006), and a land surface physical model (Minimal Advanced Treatments of Surface Interaction and Runoff – MATSIRO) (Takata et al., 2003). The biogeochemical components of land and ocean are represented by the Vegetation Integrative SImulator for Trace gases (VISIT-e) (Ito and Inatomi, 2012; Hajima et al., 2020a) and the Ocean ecosystem component embedded within the ocean circulation model (OECO2) (Yoshikawa et al., 2008; Watanabe et al., 2011; Hajima et al., 2020a), respectively. The atmospheric model has 40 vertical levels up to 3 hPa and a cutoff wavenumber of T42. The ocean model is in a horizontal tripolar coordinate system. A spherical coordinate system is used south of 63 N, with a longitude grid spacing of 1°, but meridian spacing varies from about 0.5° near the equator to 1° at mid-latitudes. The hybrid σ-z coordinate system has 62 vertical layers.

The OECO2 includes two types of phytoplankton (normal non-diazotrophic phytoplankton and diazotrophic phytoplankton), zooplankton, particulate detritus, nitrate, phosphate, dissolved iron, and dissolved oxygen. Atmospheric nitrogen and iron deposition, inorganic nitrogen and phosphorus input from rivers, dissolved iron from continental shelves and hydrothermal vents, and biological nitrogen fixation are considered as external nutrient inputs in OECO2. A land surface biogeochemical model is coupled with a river routine that can compute inorganic nitrogen transport from soil to ocean. In addition to river discharge, fertilizer, atmospheric nitrogen deposition, and land use changes, which are prescribed in the forcing data, affect riverine nitrogen fluxes. Riverine phosphorus flux from land to ocean is diagnosed from nitrogen fluxes using the Redfield ratio ($\mathrm{C}:\mathrm{N}:\mathrm{P}:\mathrm{O}=\mathrm{106}:\mathrm{16}:\mathrm{1}:\mathrm{138}$). Nitrogen deposition to land and ocean is prescribed in forcing data. Atmospheric iron deposition is predicted on the basis of dust simulated by SPRINTARS and data from biomass, fossil fuel, and biofuel combustion emissions datasets.

The equilibrium climate sensitivity (ECS) of this version of the model has been diagnosed using Gregory's method (Gregory et al., 2004) and is reported to be 2.66 K (Tsutsui, 2020; Ohgaito et al., 2021). The transient climate response to cumulative carbon emissions (TCRE) of MIROC-ES2L is 1.3 °C per 1000 Pg C (MacDougall et al., 2020; Hajima et al., 2020a), which lies within the likely range estimated by the Intergovernmental Panel on Climate Change (IPCC). Within the range of ECS values of the CMIP6 models, the ECS of MIROC-ES2L is relatively modest.

2.2 Extended Shared Socioeconomic Pathways scenarios and model experimental design

Anthropogenic greenhouse gas emissions, quantified under the Shared Socioeconomic Pathways (SSP) scenarios, are the primary drivers of current and future climate change in climate model simulations, with CO₂ being the dominant contributor to radiative forcing. The SSPs harmonize projections of $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$, ranging from approximately 393 ppm in 2015 to over 2000 ppm by 2250. They reflect a range of possible futures based on varying societal and policy choices (Figs. 1a and A1) (Meinshausen et al., 2020).

In this study, we use the extended SSPs, which contain projections of greenhouse gas concentrations up to 2500, to investigate the long-term response of the climate and carbon cycle. Eight extended scenarios are employed for concentration-driven simulations of MIROC-ES2L: SSP1–1.9, SSP1–2.6, SSP2–4.5, SSP3–7.0, SSP4–3.4, SSP4–6.0, SSP5–3.4–OS, and SSP5–8.5 (Fig. A1). This paper primarily focuses on three scenarios – SSP1–2.6, SSP2–4.5, and SSP5–8.5 – representing low, intermediate, and high emission pathways, respectively.

By 2500, net-zero emissions will be achieved across all the extended scenarios; $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ begins to decrease after 2250 even in the high-emission scenario. In the MIROC-ES2L experiment driven by these extended scenarios, global mean surface temperatures (GMST) show a tendency to stabilize from around 2250. Whereas the decline in $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ initiates global cooling, the continued warming of the deep ocean induces a gradual rise in surface air temperature, suggesting complex transient behavior beyond simple radiative response. The long-term evolution of the physical ocean state is summarized in Appendix B.

Additionally, changes in GMST levels are associated with non-linear and potentially irreversible responses in the climate system, such as tipping points in sea-ice extent, abrupt changes in El Niño–Southern Oscillation (ENSO) behavior, cloud regimes, and polar climate dynamics, as well as irreversibility in the Atlantic Meridional Overturning Circulation (AMOC) (Lenton et al., 2008; Weijer et al., 2019; Armstrong McKay et al., 2022). Notably, using the same MIROC-ES2L extended SSP simulations, Hayashi et al. (2026) demonstrated non-monotonic changes in ENSO characteristics with increasing global warming levels.

Figure 1 Surface carbon rapidly adjusts to atmospheric forcing, whereas the deep ocean remains committed to centuries of carbon accumulation. (a) Atmospheric pCO₂ forcings prescribed by the extended SSP scenarios that were used in our analyses. (b) Depth-integrated anthropogenic carbon inventory (unit: PgC) in three ocean layers (0–200, 200–1000, and 1000–6000 m) and the global ocean, simulated under the extended SSP scenarios (SSP1–2.6, SSP2–4.5, and SSP5–8.5). The anthropogenic carbon inventory is further partitioned into preformed (dashed lines) and regenerated (dotted lines) components.

2.3 Isolation of factors that change oceanic partial pressure of CO₂

The oceanic partial pressure of CO₂ is determined by temperature, salinity, DIC, and total alkalinity (Sarmiento and Gruber, 2006). Accordingly, the tendency of $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ is governed by the tendencies of these four variables. The tendency of $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ can be decomposed as follows:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}\left[p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}\right]}{\partial t}}& ={\displaystyle \frac{\partial \left[p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}\right]}{\partial \left[T\right]}}\cdot {\displaystyle \frac{\partial \left[T\right]}{\partial t}}+{\displaystyle \frac{\partial \left[p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}\right]}{\partial \left[S\right]}}\cdot {\displaystyle \frac{\partial \left[S\right]}{\partial t}}\\ & +{\displaystyle \frac{\partial \left[p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}\right]}{\partial \left[\mathrm{DIC}\right]}}\cdot {\displaystyle \frac{\partial \left[\mathrm{DIC}\right]}{\partial t}}+{\displaystyle \frac{\partial \left[p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}\right]}{\partial \left[\mathrm{ALK}\right]}}\cdot {\displaystyle \frac{\partial \left[\mathrm{ALK}\right]}{\partial t}}\end{array}\end{array}$$

The terms on the right-hand side of Eq. (1) represent the contributions from temperature, salinity, DIC, and alkalinity, respectively. Contributions for each year are computed using annual mean values of $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$. In this study, we apply this decomposition to analyze the drivers of the temporal changes in $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$.

Further decomposition of the temperature and salinity terms is given by:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial \left[p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}\right]}{\partial \left[T\right]}}\cdot {\displaystyle \frac{\partial \left[T\right]}{\partial t}}=\left[p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}\right]\cdot {\displaystyle \frac{\partial \ln \left[p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}\right]}{\partial \left[T\right]}}\cdot {\displaystyle \frac{\partial \left[T\right]}{\partial t}}\text{(3)}& {\displaystyle \frac{\partial \left[p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}\right]}{\partial \left[S\right]}}\cdot {\displaystyle \frac{\partial \left[S\right]}{\partial t}}={\displaystyle \frac{\left[p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}\right]}{\left[S\right]}}\cdot {\displaystyle \frac{\partial \ln \left[p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}\right]}{\partial \ln \left[S\right]}}\cdot {\displaystyle \frac{\partial \left[S\right]}{\partial t}}\end{array}$$

The sensitivities are $\partial \ln \left[p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}\right]/\partial \left[T\right]=\mathrm{0.0423}$ °C⁻¹ and $\partial \ln \left[p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}\right]/\partial \ln \left[S\right]=\mathrm{1.0}$ (Sarmiento and Gruber, 2006).

3 Results

3.1 Centennial-scale changes in ocean biogeochemistry

3.1.1 Anthropogenic carbon accumulation and dissolved inorganic carbon redistribution

The temporal evolutions of depth-weighted time series of anthropogenic carbon storage and DIC are depth-dependent. They are characterized by a rapid, scenario-sensitive surface adjustment, a delayed yet sustained mesopelagic response, and a monotonic deep-ocean increase across all scenarios (Figs. 1b and A2). In the low-emission (e.g., SSP1–1.9, SSP1–2.6) and overshoot (SSP5–3.4-OS) scenarios, surface DIC stabilizes at approximately historical levels once $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ begins to decline, whereas deep-ocean DIC continues to increase throughout the simulation, reflecting the long adjustment timescales of ocean carbon uptake (Fig. A2).

In the upper ocean (0–200 m), anthropogenic carbon storage and DIC closely follow $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$. Both increase rapidly during the historical period and projections vary according to scenario after approximately 2100. Under SSP1–2.6, surface DIC reaches a maximum during the 21st century before gradually decreasing toward historical levels. Under SSP2–4.5, DIC remains elevated despite a slow decline. By contrast, under SSP5–8.5, anthropogenic carbon storage and DIC increase rapidly throughout the 21st century and remain substantially above historical levels thereafter (Figs. 1b and A2).

Between 200 and 1000 m, DIC concentrations eventually stabilize or gradually decline under the low-emission and overshoot scenarios, whereas they continue to increase or remain elevated under the high-emission scenario. Partitioning anthropogenic carbon storage into preformed and regenerated components using apparent oxygen utilization shows that the long-term increase is dominated by the preformed component. Nevertheless, the regenerated component also increases under all scenarios (Fig. 1b). This delayed response is consistent with slow ventilation and the long residence times of intermediate waters (see Appendix B and Fig. A3).

Below 1000 m, both anthropogenic carbon storage and DIC continue to increase monotonically through 2500, with no indication of reversal in any scenario. This persistent deep-ocean accumulation despite stabilization or decline of $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ reflects substantial oceanic inertia and the slow downward propagation of anthropogenic carbon.

Taken together, these results demonstrate a pronounced vertical contrast in ocean carbon redistribution. Although emission mitigation can halt or reverse changes in the surface ocean, anthropogenic carbon continues to accumulate in the deep ocean over multi-centennial timescales, committing the ocean carbon reservoir to prolonged adjustment long after atmospheric forcing begins to stabilize.

3.1.2 Biological productivity and nutrient redistribution

The export of carbon from the euphotic zone is a key component of the ocean carbon cycle, influencing the transfer and storage of carbon in the deep ocean. Long-term responses of primary production and the vertical export of detritus at 100 m vary according to scenario, with sustained increase in some scenarios and pronounced declines in others (Fig. 2a and b).

Under low-emission scenarios (SSP1–2.6 and SSP2–4.5), biological productivity increases approximately linearly between the beginning of the 21st century and 2500. By contrast, under the high-emission scenario SSP5–8.5, biological productivity is projected to decline by approximately 5 % by the early 22nd century, with only limited recovery thereafter.

The temporal evolution of detritus export is influenced by a complex interplay of environmental factors, including temperature, light availability, and nutrient concentrations, whereas the basin-scale distributions of these variables are reshaped by changes in ocean circulation. In particular, weakening of the AMOC (Fig. A3) reduces the upward transport of nutrients to the North Atlantic surface ocean, leading to lower surface macronutrient concentrations (Fig. A4) and a suppression of detritus export (Fig. 2d). Meanwhile, nutrients become increasingly trapped within the lower overturning cell associated with Antarctic Bottom Water (AABW). This leads to a pronounced decline in surface and subsurface nutrient concentrations in the Pacific Ocean after 2250 under SSP5–8.5.

The relationships between $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$, detritus export, and AMOC strength are consistent with the mechanisms described above (Fig. 3). Increasing $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ is associated with progressive weakening of the AMOC and reduction in detritus export. In particular, detritus export decreases below its preindustrial level when $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ reaches approximately 3 times the preindustrial concentration (Fig. 3a).

Figure 2 Future biological productivity diverges sharply across emission pathways. (a) Detritus downward flux at the depth of 100 m (unit: Pg C yr⁻¹) simulated under SSP1–2.6, SSP2–4.5, and SSP5–8.5. Thick lines represent 10-year running means, whereas thin lines indicate annual-mean values. (b) Gross primary production (unit: Pg C yr⁻¹) simulated under SSP1–2.6, SSP2–4.5, and SSP5–8.5. Thick lines represent 10-year running means, whereas thin lines indicate annual values. Spatial changes in detritus downward flux at the depth of 100 m (unit: $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), calculated as the difference between the values from the 2240s and the historical values from the 1990s under (c) SSP1–2.6 and (d) SSP5–8.5. The bar plots exhibit integrated air-sea CO₂ fluxes over the regions (unit: Pg C yr⁻¹) used in Mikaloff Fletcher et al. (2007).

Figure 3 Atmospheric pCO₂ exhibits a strong negative relationship with organic carbon export, which is in turn closely linked to the strength of the AMOC. (a) Scatter plot showing the relationship between $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ (ppm) and organic carbon export (PgC yr⁻¹). Each point represents a 50-year mean (2015–2050, 2051–2100, …, 2451–2500) simulated under the extended SSP scenarios. (b) Scatter plot showing the relationship between $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ (ppm) and the maximum AMOC strength (Sv). Each point represents a 50-year mean (2015–2050, 2051–2100, …, 2451–2500) simulated under the extended SSP scenarios.

3.2 Processes controlling the ocean carbon sink

3.2.1 Evolution of surface-ocean partial pressure of CO₂

Across the different scenarios, $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ begins to decline under the prescribed forcing by 2250 (Fig. 1a). Subsequently, global mean temperature stabilizes (Fig. A5). Long-term projections indicate leveling off of both $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ and temperature but continued dynamic responses of the ocean carbon cycle to the prescribed forcing. The oceanic pCO₂ in the surface ocean continues to vary over the long term in response to the physical and biogeochemical changes in the ocean. Under the extended SSP scenarios, surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ generally follows the trajectory of $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$, with relatively small differences across latitudes (Figs. 1 and A6)

3.2.2 Air–sea difference in partial pressure of CO₂

Changes in surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ strongly affect the air–sea difference in pCO₂, defined as $\mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2}}=p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}-p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$. Negative values indicate that the surface ocean is undersaturated with respect to the atmosphere and therefore favors oceanic CO₂ uptake, whereas positive values imply supersaturation and potential outgassing.

At the global scale (Fig. 4a and b), the temporal evolution of ΔpCO₂ differs markedly between scenarios. From the beginning of the 21st century to around 2100, ΔpCO₂ increases from approximately −10 to nearly 0 ppm under SSP1–2.6, reflecting the relatively slow increase in $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$. By contrast, under SSP5–8.5, ΔpCO₂ decreases from approximately −10 to about −15 ppm by around 2080, indicating that the rapid rise in $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ outpaces the increase in surface-ocean pCO₂. After approximately 2080, ΔpCO₂ begins to increase as surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ gradually catches up with atmospheric forcing.

The two scenarios also differ in terms of the evolution of the spatial distribution of ΔpCO₂ (Fig. 4c and d). Under SSP1–2.6, positive ΔpCO₂ remains largely confined to the equatorial Pacific throughout the simulation. By contrast, under SSP5–8.5, positive ΔpCO₂ progressively expands across the equatorial and subtropical oceans during the 21st century. By the end of the 22nd century, supersaturated waters extend over much of the subtropical gyres and subsequently spread into the Northern Hemisphere subpolar oceans during the late 23rd century. This poleward expansion indicates that surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ in increasingly large areas of the surface ocean approaches equilibrium with, or exceeds, $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ under sustained high emissions.

Figure 4 Under SSP5–8.5, oceanic pCO₂ exceeds atmospheric pCO₂ at low latitudes, whereas the opposite applies in the Southern Ocean. Annual-mean changes in the ocean-atmosphere pCO₂ difference (ΔpCO₂, unit: ppm) by latitude band, simulated under (a) SSP1–2.6 and (b) SSP5–8.5. Positive values indicate that oceanic pCO₂ exceeds atmospheric pCO₂, whereas negative values indicate the opposite. Colored lines represent different latitude bands as follows: gray (global mean), dark blue (equatorial: 15° S–15° N), orange (northern subtropics: 15–40° N), yellow (southern subtropics: 15–45° S), pink (northern subpolar: 40–60° N), green (Arctic: >60° N), and light green (Antarctic: <45° S). Annual-mean ΔpCO₂ for (c) SSP1–2.6 and (d) SSP5–8.5. Panels show the decadal means for the 2040s, 2090s, 2240s, and 2490s.

3.2.3 Air–sea CO₂ flux

The temporal evolution of the simulated air–sea CO₂ flux broadly follows that of ΔpCO₂ (Fig. 5). Under SSP5–8.5, oceanic CO₂ uptake increases throughout the 21st century, reaches a maximum near the end of the century, and gradually decreases thereafter. This evolution is consistent with the temporal changes in ΔpCO₂, indicating that changes in the air–sea pCO₂ gradient dominate the long-term evolution of the global ocean carbon sink. By contrast, changes in the air–sea CO₂ flux remain comparatively small under SSP1–2.6.

The spatial distribution of the air–sea CO₂ flux differs in several important aspects from that of ΔpCO₂ (Figs. 5c–d and 4c–d). Under SSP1–2.6, oceanic CO₂ uptake is maintained over much of the extratropical ocean, whereas outgassing remains largely confined to low latitudes. Under SSP5–8.5, positive ΔpCO₂ expands across much of the tropical and subtropical oceans, but this expansion is not accompanied by a corresponding increase in air–sea CO₂ outgassing. Instead, the strongest outgassing is confined in the region near the equatorial oceans.

The Southern Ocean remains a particularly strong carbon sink throughout the SSP5–8.5 simulation. Although the spatial distribution of the air–sea CO₂ flux broadly follows that of ΔpCO₂, the Southern Ocean exhibits substantially stronger uptake than would be inferred from ΔpCO₂ alone. This regional enhancement is consistent with the persistently strong westerly winds over the Southern Ocean (Fig. A7), which may enhance air–sea gas exchange.

The persistence of oceanic carbon uptake under SSP5–8.5 is consistent with the multi-model results of Koven et al. (2022), which showed that the global ocean remains a carbon sink under sustained high $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$; it transitions to a carbon source primarily under concentration-driven overshoot scenarios with declining $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$. The sustained Southern Ocean carbon sink that is present in our simulations is therefore consistent with this broader long-term behavior of the global ocean carbon cycle.

Figure 5 Under SSP5–8.5, the oceanic CO₂ flux indicates gradual shifts toward outgassing at the low latitudes, whereas the Southern Ocean remains a persistent sink. Annual-mean changes in the CO₂ flux (unit: g C yr⁻¹) by latitude band, simulated under (a) SSP1–2.6 and (b) SSP5–8.5. Positive values indicate oceanic CO₂ outgassing, whereas negative values indicate oceanic CO₂ uptake. Colored lines represent different latitude bands as follows: gray (global mean), dark blue (equatorial: 15° S–15° N), orange (northern subtropics: 15–40° N), yellow (southern subtropics: 15–45° S), pink (northern subpolar: 40–60° N), green (Arctic: >60° N), and light green (Antarctic: <45° S). Annual-mean CO₂ flux (unit: $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) for (c) SSP1–2.6 and (d) SSP5–8.5. Panels show the decadal means for the 2040s, 2090s, 2240s, and 2490s.

3.2.4 Carbonate buffering capacity of the surface ocean

The ocean plays a buffering role in the global carbon cycle by absorbing CO₂ and, through seawater carbonate chemistry, buffers the increase in acidity associated with elevated DIC, thereby mitigating the accumulation of CO₂ in the atmosphere. However, ongoing ocean acidification progressively weakens this carbonate buffering capacity, resulting in increased temporal and spatial variability in the carbonate system's response to changes in CO₂ (Orr et al., 2005; Egleston et al., 2010).

To characterize these long-term changes, we first examine the evolution of the Revelle factor, which is the ratio of the relative change in pCO₂ to the relative change in DIC and is the conventional measure of the sensitivity of seawater pCO₂ to changes in DIC. Because the Revelle factor becomes less sensitive under high-CO₂ conditions, we calculate the buffering capacity, $\partial \mathrm{DIC}/\partial p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$, which has been adopted as a practical indicator of seawater buffering capacity (Riebesell et al., 2009). The buffering capacity is not equivalent to the reciprocal of the Revelle factor, but we use it as a proxy that provides a direct measure of the increase in DIC associated with a given increase in seawater pCO₂.

Figure 6a and c shows that the buffering capacity in the surface ocean is projected to decline markedly during the 21st and 22nd centuries even under the low-emission scenario SSP1–2.6. However, gradual recovery is expected to emerge by the end of the 22nd century and the buffering capacity is projected to increase again over multi-centennial timescales. This delayed recovery reflects the stabilization and subsequent decline of atmospheric and surface ocean pCO₂ under SSP1–2.6; $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ reaches levels that are close to those of the present day only after around 2200. The carbonate system re-equilibrates slowly and buffering capacity is partially restored.

By contrast, under the high-emission scenario SSP5–8.5, the buffering capacity continues to decline throughout the simulation period and reaches its effective minimum values of approximately 0.05 across almost the entire global ocean by around 2200 (Fig. 6b and d). This indicates a state in which carbonate-system buffering is strongly diminished as continued accumulation of anthropogenic carbon drives buffering capacity toward its effective minimum across most of the global surface ocean. As buffering capacity approaches this effective minimum, regional contrasts progressively diminish, resulting in a nearly homogeneous spatial distribution of buffering capacity by around 2200.

When buffering capacity is high, increases in DIC lead to relatively moderate rises in pCO₂. However, as buffering capacity declines, the sensitivity of pCO₂ to a given increase in DIC increases, amplifying the effects of other drivers such as temperature rise and reductions in alkalinity (see Sect. 3.2.5). Changes in buffering capacity thus render the surface ocean more responsive to atmospheric forcing, particularly to the continued uptake of anthropogenic CO₂. A reduced buffering capacity also implies a diminished efficiency in oceanic carbon uptake, increasing the fraction of anthropogenic CO₂ that remains in the atmosphere. This point is discussed in more detail in the Discussion (Sect. 4.3).

Figure 6 Under SSP5–8.5, the ocean's CO₂ buffering capacity reaches a near minimum by the mid–22nd century, whereas even under SSP1–2.6 it takes about a century to return to present-day levels. Annual-mean changes in $\mathrm{dDIC}/\mathrm{d}p{\mathrm{CO}}_{\mathrm{2}}$ – representing the sensitivity of dissolved inorganic carbon (DIC) to changes in surface-ocean pCO₂ – in the surface ocean (unit: ND) by latitude band, simulated under (a) SSP1–2.6 and (b) SSP5–8.5. Colored lines represent different latitude bands as follows: gray (global mean), dark blue (equatorial: 15° S–15° N), orange (northern subtropics: 15–40° N), yellow (southern subtropics: 15–45° S), pink (northern subpolar: 40–60° N), green (Arctic: >60° N), and light green (Antarctic: <45° S). Annual-mean $d\mathrm{DIC}/dp{\mathrm{CO}}_{\mathrm{2}}$ for (c) SSP1–2.6 and (d) SSP5–8.5. A higher value indicates stronger buffering capacity, meaning that the ocean can absorb and retain more CO₂ for a given increase in pCO₂. Panels show the decadal means for the 2040s, 2090s, 2240s, and 2490s.

3.2.5 Drivers of temporal changes in surface-ocean partial pressure of CO₂

Changes in $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ are driven by variations in ocean temperature, salinity, DIC, and total alkalinity. We performed a decomposition analysis of the temporal changes in surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ to identify the contributions of these factors (see Sect. 2.3 for details) (Figs. 7 and A8).

Under the low-emission scenario SSP1–2.6, temporal variations in surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ are largely explained by changes in DIC (Fig. 7a). Under this scenario, surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ increases until around 2080 and subsequently decreases. The positive trend becomes negative as $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ stabilizes. By around 2200, the change in surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ approaches zero, indicating near equilibrium between the ocean and atmosphere.

Similarly, from the beginning of the 21st century to around 2200, the rise in surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ under the high-emission scenario SSP5–8.5 is primarily driven by increasing DIC concentrations, reflecting the rapid uptake of anthropogenic carbon under strong air–sea disequilibrium (Fig. 7b). Before 2200, the contributions from temperature and alkalinity changes remain secondary. After 2200, however, the influence of reductions in surface alkalinity, together with thermal effects, becomes increasingly important in driving further increases in surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$. This influence is associated with enhanced upper-ocean stratification and the subsequent weakening of the meridional overturning circulation (Fig. A9). During this period, although calcium carbonate export decreases, the associated conservation of alkalinity is insufficient to offset the effects of physical changes, resulting in an overall decline in surface-ocean alkalinity. This coincides with the approach of the buffering capacity of the marine carbonate system toward its minimum, consistent with the interpretation of Chikamoto et al. (2023).

Figure 7 Changes in surface-ocean pCO₂ are initially driven by surface DIC, with alkalinity becoming important from the mid–22nd century onward under SSP5–8.5. Contributions of temperature (red), salinity (blue), DIC (yellow), and alkalinity (green) to the temporal change of global-mean surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ (black) under (a) SSP1–2.6 and (b) SSP5–8.5. Thick lines represent 20-year running means, whereas thin lines indicate 5-year mean values. The horizontal axis represents time (years), and the vertical axis shows the time derivative of pCO₂ ($\mathrm{d}p{\mathrm{CO}}_{\mathrm{2}}/\mathrm{d}t$, unit: ppm yr⁻¹).

3.2.6 Deep-ocean carbon storage and remaining carbon uptake capacity

Under both low- and high-emission scenarios, DIC concentrations in the deep ocean continue to increase even after $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ stabilizes or begins to decline (Figs. 1a, A1 and A2). This sustained increase is consistent with the long residence time and slow renewal of deep waters, which allow anthropogenic carbon to accumulate and persist in the ocean interior for centuries once it is transported below the surface. Enhanced stratification and slow, multi-centennial adjustments of ocean biogeochemical cycles further support this continued accumulation.

In particular, pronounced increases in DIC are evident along major deep-water formation and export pathways under SSP5–8.5 (Fig. A10). Along AABW pathways, DIC increases extend into the abyssal and bottom ocean, whereas the North Atlantic Deep Water (NADW) exhibits enhanced DIC accumulation primarily within the deep Atlantic basin. These spatial patterns indicate the efficient transfer of anthropogenic carbon from the surface into the ocean interior along established deep-water circulation pathways.

Changes in the formation and properties of deep waters, together with modifications of the circulation in the ocean interior, regulate the depth at which carbon is stored and the residence time of deep waters. Consequently, they are associated with variations in the vertical distribution and accumulation of DIC. In addition, changes in biological export and remineralization contribute to enhanced DIC accumulation in the ocean interior. Under SSP5–8.5, regenerated carbon increases despite the general decrease in detritus export (Fig. 2d). This is because longer water-mass residence times enhance the accumulation of remineralized carbon in the ocean interior, which is consistent with previous CMIP6 analyses of the future biological carbon pump (Wilson et al., 2022). Interior DIC increase is dominated by preformed carbon, with regenerated carbon providing a secondary but positive contribution (Fig. 1b).

To further examine the remaining capacity of the ocean interior to absorb anthropogenic CO₂, we calculate potential $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ ($P{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$), defined as the surface-equivalent pCO₂ of an interior seawater parcel (Broecker and Peng, 1982; Chen et al., 2022). The deviation of interior $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ from $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ (ΔPCO₂) is used as a diagnostic for the remaining chemical capacity of the ocean interior to absorb anthropogenic CO₂.

Under SSP1–2.6, elevated ΔPCO₂ develops at intermediate depths and progressively extends into the deep ocean (Figs. 8a–d). By the mid–23rd century, positive ΔPCO₂ is evident even within the abyssal ocean, indicating that anthropogenic carbon has penetrated efficiently into the ocean interior and progressively reduced the remaining capacity for additional carbon uptake.

By contrast, under SSP5–8.5, ΔPCO₂ exhibits a pronounced vertical gradient, with elevated values in the thermocline but persistently low values throughout much of the deep ocean (Fig. 8e–h). This pattern indicates that the rapid increase in $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ outpaces the ventilation of the deep ocean, leaving much of the ocean interior substantially undersaturated with respect to anthropogenic CO₂ despite continued increases in interior DIC. In the Southern Ocean, relatively low ΔPCO₂ extends to shallow waters further highlighting the important role of Southern Ocean ventilation in sustaining the long-term carbon storage capacity of the deep ocean. By the 2240s, large positive ΔPCO₂ values develop within the mesopelagic ocean. This is consistent with the nonlinear response of the marine carbonate system to anthropogenic carbon invasion and the continued accumulation of remineralized DIC (Fassbender et al., 2023).

Figure 8 Under SSP5–8.5, surface and subsurface waters become potentially supersaturated in pCO₂ relative to the atmosphere over time, whereas deeper waters remain potentially undersaturated. Zonal mean vertical sections of $\mathrm{\Delta }P{\mathrm{CO}}_{\mathrm{2}}=P{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}-p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ (unit: ppm) for (a–d) SSP1–2.6 and (e–h) SSP5–8.5, where $P{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ denotes ocean potential pCO₂ and $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ denotes atmospheric pCO₂. From left to right, the panels show longitudinally averaged sections for the Atlantic, Southern, and Pacific Oceans. From top to bottom, rows represent the decadal means for the (a, e) 2040s, (b, f) 2090s, (c, g) 2240s, and (d, h) 2490s.

3.3 Ecological consequences of long-term biogeochemical change

As a result of ongoing global warming, marine ecosystems are under increasingly complex and regionally differentiated pressures from concurrent stressors – e.g., warming, acidification, and deoxygenation – that vary in magnitude across time and space. We evaluate ecosystem stress using multiple indicators, including SST, nutrient concentrations, dissolved oxygen, pH, and aragonite and calcite saturation states (Ω_arag and Ω_calc), following Kwiatkowski et al. (2020), over centennial timescales (Figs. 9, 10 and A11).

Ecosystem stress under the low-emission scenario SSP1–2.6 remains relatively limited compared with that under the higher-emission scenarios. Anomalies in temperature, ocean acidification, and deoxygenation persist over several centuries, primarily in the polar and subpolar regions (Fig. 9a–c). Nevertheless, the return of multiple stress indicators toward modern-day levels is slow across a wide range of marine ecosystems, including coastal seas, coral reefs, and kelp forests. Present-day conditions are only reached in the late 22nd century (Fig. 10).

Importantly, several key environmental variables do not fully return to present-day levels even after $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ declines under SSP1–2.6. Sea surface temperatures remain approximately 0.5 °C above present-day values after the 21st century. Surface nutrient concentrations, dissolved oxygen, and alkalinity continue to decrease over multiple centuries (Figs. A12, A13, and A14). These residual anomalies are particularly pronounced in the Arctic, where changes in temperature and nutrient availability persist through to the end of the 25th century (Fig. 10). Together, these results indicate that even under a low-emission pathway, key physical and biogeochemical stress indicators can take several centuries to return to present-day levels.

By contrast, under SSP2–4.5, combined ecosystem stress arising from surface warming and declining nutrient availability expands beyond the polar regions into the subpolar and equatorial regions over multi-centennial timescales (Fig. 9d–f). Notably, this expansion continues well beyond 2100, despite only modest increases in $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ after the end of the 21st century. On multi-centennial timescales, some stressors show partial relaxation; however, conditions do not return to present-day levels and, at best, approach late–21st-century states (Fig. 10). This behavior indicates that ecosystem stress under SSP2–4.5 is driven not only by contemporaneous atmospheric forcing, but also by delayed oceanic responses, including the persistence of elevated surface temperatures and sustained reductions in upper-ocean nutrient supply. As a result, the spatial extent of stressed marine ecosystems continues to grow over subsequent centuries, highlighting a cumulative and lagged response of the ocean–ecosystem system to the forcing under intermediate-emission scenarios.

This trend is further amplified under the high-emission scenario SSP5–8.5 (Fig. 9g–i). Up to around 2090, the evolution of ecosystem stress is broadly consistent with previous assessments focusing on the 21st century (Kwiatkowski et al., 2020), with stress primarily emerging in the tropical and high-latitude regions. After 2100, however, ecosystem stress intensifies markedly, as the combined effects of surface warming and declining nutrient availability expand across nearly all ocean basins. This intensification is accompanied by pronounced surface-ocean acidification on a global scale.

Under SSP5–8.5, declining surface nutrient availability in coastal and kelp systems and strong reductions in surface aragonite saturation in coral reef regions become increasingly pronounced on multi-centennial timescales (Fig. 10). Together, these results highlight a substantial increase in ecosystem stress beyond the end of the 21st century because of interactions between multiple stressors under sustained forcing in high-emission scenarios.

Spatial heterogeneity is also evident in the progression of ocean deoxygenation under SSP5–8.5 (Figs. A13 and A15). The oxygen minimum zones (OMZs) in the equatorial region exhibit relatively little change. However, there are pronounced declines in dissolved oxygen concentrations in the deep North Pacific and Southern Ocean, associated with the weakening of deep ocean ventilation and warming.

Regarding the aragonite saturation state, projections under the low-emission scenario SSP1–2.6 show that aragonite undersaturation (Ω_arag<1) emerges transiently in the polar regions around the end of the 21st century (Fig. A16a and b). Subsequently, much of the area experiencing undersaturated conditions recovers to Ω_arag>1 as $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ declines, consistent with the partial recovery shown in Fig. 10. However, aragonite saturation in parts of the Arctic Ocean remains close to unity on multi-centennial timescales, indicating a slow and incomplete return to present-day conditions.

By contrast, under the high-emission scenario SSP5–8.5, undersaturated conditions (Ω_arag<1) emerge by 2100 across most of the Southern Ocean south of 45° S, as well as over large areas of the Arctic Ocean and the subpolar North Pacific (Fig. A16e and f). These spatial patterns are broadly consistent with previous projections focusing on the 21st century; however, the extended simulations reveal a continued and substantial expansion of undersaturated regions beyond 2100. After 2100, Ω_arag<1 conditions expand further under SSP5–8.5, progressively encompassing more than half of the global ocean by the mid–23rd century and continuing to grow through the late 25th century (Fig. A16g and h). By this time, the global mean aragonite saturation state approaches, and in some periods falls below, unity, indicating a widespread transition toward corrosive conditions for aragonite across much of the ocean.

Regarding the calcite saturation state, Ω_calc remains substantially higher than Ω_arag under most conditions, such that calcite undersaturation provides a more stringent indicator of corrosive conditions for calcifying organisms. Under the low-emission scenario SSP1–2.6, Ω_calc stays well above unity across most of the global ocean, generally exceeding values of about 1.5, with near-unity conditions appearing only locally in limited coastal parts of the Arctic Ocean (Figs. A16i–l). This indicates that, although aragonite may transiently become undersaturated in the polar regions, the carbonate system does not broadly transition into conditions that are corrosive with respect to calcite under strong emission mitigation.

By contrast, under the high-emission scenario SSP5–8.5, calcite undersaturation (Ω_calc<1) emerges by the end of the 21st century over large portions of the Arctic Ocean and in coastal sectors of the Southern Ocean (Fig. A16n). The extended simulations further reveal that these Ω_calc<1 regions continue to expand beyond 2100, spreading into subpolar waters by the mid–23rd century (Fig. A16o and p). This progression indicates that sustained high emissions drive the ocean not only toward widespread aragonite corrosivity, but also toward a more extreme regime in which even calcite becomes thermodynamically unstable over broad high-latitude and subpolar regions.

Figure 9 Even under the intermediate scenario SSP2–4.5, ecosystem stress continues to intensify beyond the 21st century. Ecosystem stressors under (a–c) SSP1–2.6, (d–f) SSP2–4.5, and (g–i) SSP5–8.5. Panels show changes in the (a, d, g) 2040s, (b, e, h) 2090s, and (c, f, i) 2240s relative to the 2010s average under SSP1–2.6. Following Kwiatkowski et al. (2020), stressors are defined using sea surface temperature, surface ocean pH, euphotic-zone nitrate concentration (averaged over 0–100 m), and subsurface dissolved oxygen concentration (averaged over 100–600 m). Stressors include warming (+2 °C), nutrient decline ($-\mathrm{1}\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}$, blue), and their overlap (red). Additional stressors are ocean acidification (surface pH decrease of −0.2, hatched) and deoxygenation (subsurface dissolved oxygen decrease of $-\mathrm{30}\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}$, dotted).

Figure 10 Under SSP5–8.5, nutrient declines in coastal and kelp systems and reductions in aragonite saturation in coral reef regions become pronounced. Projected time series of coastal-ocean temperature, pH, DIC, O₂, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, and Ω_aragonite across coastal seas, coral reefs, kelp systems, upwelling zones, and polar regions (north and south) under SSP1–2.6, SSP2–4.5, SSP3–7.0, and SSP5–8.5. Spatial masks follow IPCC AR6 WGII Fig. 3.7 (Cooley et al., 2022).

4 Discussion

4.1 Mechanisms governing long-term ocean biogeochemical responses

This study examines ocean biogeochemical responses under low (SSP1–2.6), intermediate (SSP2–4.5), and high (SSP5–8.5) emission pathways extended to the year 2500, with the aim of constraining the magnitude, persistence, and reversibility of long-term ocean responses to different atmospheric CO₂ trajectories. In this section, we focus on discussing the contrasting responses between the end-member scenarios SSP1–2.6 and SSP5–8.5 and use SSP2–4.5 to illustrate intermediate and delayed adjustments.

Previous studies have primarily focused on ocean changes during the 21st century, highlighting reductions in carbonate buffering capacity, ocean acidification, and increasing ecosystem stress under continued anthropogenic emissions (Kwiatkowski et al., 2020; Masson-Delmotte et al., 2021). Our simulations of ocean changes between the 21st century and 2500 demonstrate that many of these responses are likely to continue to evolve long after atmospheric forcing has stabilized or begun to decline. Rather than approaching a steady state, the ocean is projected to progressively shift toward distinct long-term biogeochemical regimes. Regime characteristics depend strongly on the emission pathway. In particular, the high-emission scenario is characterized by a progressive approach toward the minimum buffering capacity of the marine carbonate system, continued sequestration of anthropogenic carbon in the deep ocean, and persistent changes in nutrient distributions and marine ecosystem stress. Even under strong emission mitigation, substantial perturbations to ocean chemistry and the ocean interior are projected to remain for centuries.

In the following subsections, we discuss the mechanisms underlying these long-term responses. First, we consider the transient evolution of surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$. This is followed by a discussion of the long-term evolution of carbonate buffering capacity and carbonate mineral saturation. We then consider the mechanisms controlling deep-ocean carbon sequestration and storage capacity, circulation-driven changes in nutrient supply and biological productivity, and finally focus on the implications of multiple stressors and their interactions for marine ecosystems.

4.1.1 Mechanisms controlling the temporal evolution of ocean carbon uptake

The temporal evolution of the ocean carbon sink is governed by the combined effects of atmospheric forcing, carbonate chemistry, and ocean circulation. These processes regulate the air–sea pCO₂ difference (Fig. 4) and consequently the magnitude of air–sea CO₂ exchange (Fig. 5).

Under SSP1–2.6, net oceanic CO₂ uptake is maintained throughout the 21st century because $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ exceeds surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ over most of the global ocean. As $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ subsequently stabilizes and declines, the air–sea pCO₂ difference progressively weakens, leading to a gradual reduction in ocean carbon uptake. By contrast, under SSP5–8.5, sustained atmospheric forcing maintains a large negative air–sea pCO₂ difference well beyond 2100. Consequently, oceanic CO₂ uptake continues to increase until the late 21st century before gradually decreasing as surface-ocean pCO₂ catches up with atmospheric forcing.

Although the temporal evolution of the air–sea CO₂ flux closely follows that of the air–sea pCO₂ difference, regional discrepancies remain. In particular, the Southern Ocean maintains a stronger carbon sink than would be expected from the air–sea pCO₂ difference alone. This behavior is consistent with the enhancement of air–sea gas exchange by persistently strong Southern Hemisphere westerlies (Fig. A7). It is also in alignment with the multi-model analysis of Koven et al. (2022), which indicated that the Southern Ocean remains an important carbon sink under sustained high-emission scenarios.

The long-term decline in carbonate buffering capacity provides an additional control on the temporal evolution of ocean carbon uptake. As anthropogenic carbon accumulates, the Revelle factor increases and the buffering capacity decreases, indicating progressively larger increases in surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ for a given increase in DIC. By the 22nd century under SSP5–8.5, the buffering capacity of much of the surface ocean approaches its minimum. As a result, the efficiency of further oceanic CO₂ uptake is substantially reduced (Fig. 6b).

This simulated slowdown in ocean carbon uptake under SSP5–8.5 is consistent with the late–21st-century saturation of oceanic carbon uptake reported by Chikamoto and DiNezio (2021) who attributed the slowdown primarily to reduced CO₂ solubility and carbonate-system buffering associated with ocean warming and anthropogenic carbon accumulation. Furthermore, Chikamoto et al. (2023) identified declining surface-ocean alkalinity as an additional feedback that operates on multi-centennial timescales. This is consistent with the results of our decomposition analysis, which show reductions in surface alkalinity becoming an increasingly important driver of surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ after the mid–22nd century under SSP5–8.5 (Fig. 7b).

There are, nevertheless, regional differences in terms of the timing and mechanisms of this feedback. Whereas Chikamoto et al. (2023) emphasized alkalinity reductions associated with Southern Ocean processes, our simulations indicate that alkalinity-driven increases in surface-ocean pCO₂ are likely to appear earlier in the Arctic and the Northern Hemisphere subpolar oceans (Fig. A8b). These regional differences are consistent with changes in AMOC strength and associated tracer transport, suggesting that circulation pathways exert an important control on the spatial expression of long-term alkalinity feedbacks. Our simulations indicate that the dominant mechanism underlying these differences is the reduced upward transport of high alkalinity waters to the surface caused by enhanced upper-ocean stratification and changes in overturning circulation, with possible contributions from reduced calcium carbonate export (Fig. A9).

4.1.2 Long-term evolution of ocean buffering capacity and carbonate mineral saturation

Anthropogenic uptake of atmospheric CO₂ fundamentally alters the ocean carbonate system by reducing carbonate ion concentrations and the buffering capacity of seawater, which we use as a proxy of defined here as the reciprocal of the Revelle factor. Consistent with earlier coupled climate–carbon cycle simulations (e.g. Orr et al., 2005), our results demonstrate a pronounced decline in buffering capacity throughout the 21st and 22nd centuries, with strong spatial heterogeneity and a clear dependence on emission pathway.

Across all scenarios, the loss of buffering capacity emerges earliest in the Arctic Ocean and the subpolar regions in the Northern Hemisphere (Fig. 6). Low temperatures, relatively low alkalinity, and efficient air–sea gas exchange result in high intrinsic sensitivity to carbon perturbations in these regions, which can, therefore, be considered as early-warning hotspots of carbonate-system vulnerability. Under SSP1–2.6, buffering capacity declines substantially during the 21st and early 22nd centuries, but returns gradually toward present-day values at the global scale as atmospheric and surface-ocean pCO₂ values stabilize and decline (Fig. 6a). Nevertheless, reduced buffering capacity persists at high latitudes. Recovery of buffering capacity is incomplete even under strong emission mitigation.

Under SSP5–8.5, buffering capacity continues to decline throughout the simulation period and converges toward its effective minimum values across almost the entire global ocean by the end of the 22nd century (Fig. 6b). This convergence marks a transition toward a weakly buffered carbonate-system regime, in which further increases in DIC generate disproportionately large increases in surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$. As buffering weakens, the sensitivity of surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ to additional carbon uptake increases, amplifying the effects of ocean warming and alkalinity changes.

Changes in buffering capacity are closely mirrored by the evolution of the aragonite saturation state, Ω_arag, highlighting the causal link between carbonate-system sensitivity and the emergence of corrosive conditions for calcium carbonate minerals. Under SSP1–2.6, transient aragonite undersaturation develops in the polar regions toward the end of the 21st century, but much of this area subsequently returns to supersaturated conditions (Fig. A16a and b). By contrast, under SSP5–8.5, aragonite undersaturation expands beyond the high latitudes after 2100, eventually extending over more than half of the global ocean on multi-centennial timescales (Fig. A16e–h).

Our simulations take into account prognostic alkalinity and dynamic air–sea CO₂ exchange, but not pH-dependent changes in CaCO₃ production and dissolution. In addition, they do not include sediment carbonate dissolution or associated benthic alkalinity feedbacks, which could partially offset surface-ocean alkalinity declines and modify the long-term evolution of ocean buffering capacity and carbon uptake on multi-centennial timescales. Recent work suggests that such carbonate pump feedbacks may partially counteract acidification by increasing alkalinity (Planchat et al., 2024). Consequently, our simulation results may represent overestimates of the long-term reduction in surface-ocean buffering capacity, particularly under sustained high-CO₂ conditions. Under these conditions, sedimentary alkalinity feedbacks are expected to become increasingly important.

4.1.3 Controls on deep-ocean carbon sequestration and long-term storage

Across all scenarios, DIC concentrations in the deep ocean continue to increase even after $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ stabilizes or declines, reflecting the long residence times and slow renewal of deep waters. Once anthropogenic carbon is transported below the surface, it is sequestered in the ocean interior for centuries to millennia and becomes a long-lived memory of past emissions in the climate system.

The implications of this storage are clarified using $P{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$. Potential oceanic pCO₂ isolates the biogeochemical imprint on subsurface carbon and can be used as a diagnostic of the ocean interior's potential to absorb or release CO₂ upon ventilation.

Under SSP1–2.6, $P{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ develops pronounced maxima at intermediate depths and exceeds $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ even in the deep ocean by the mid–23rd century (Fig. 8a–d). This indicates that, although strong emission mitigation limits perturbations in the surface ocean, much of the carbon stored in the ocean interior would have the potential to outgas if deep waters were re-exposed to the surface, highlighting the delayed reversibility of oceanic carbon uptake.

Under SSP5–8.5, $P{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ exhibits strong vertical gradients (Fig. 8e–h). Elevated values are confined to the thermocline and values in the deep ocean remain persistently low, indicating that the rapid rise in $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ outpaces deep ocean ventilation. As a result, the deep ocean remains far from chemical equilibrium with the atmosphere despite continued DIC accumulation. Consequently, the deep ocean retains substantial capacity for further CO₂ uptake, highlighting that carbon storage and chemical equilibration proceed on different timescales under sustained high emissions.

Examination of the preformed and regenerated components of deep-ocean carbon storage further highlights the complementary roles of the physical and biological pumps in long-term carbon sequestration (Fig. 1b). Although detritus export declines under the high-emission scenarios, regenerated carbon storage continues to increase. This contrasting behavior suggests that long-term regenerated carbon storage is influenced not only by biological export but also by changes in ocean circulation and water-mass residence time. Nevertheless, preformed carbon remains the dominant contributor to the long-term increase in deep-ocean carbon storage. These results suggest that changes in ocean circulation remain a primary control on long-term carbon sequestration, whereas biological processes continue to influence the redistribution and storage of carbon within the ocean interior.

The accumulation of DIC along AABW pathways further underscores the importance of deep-ocean circulation in regulating long-term carbon storage. Freshwater input from ice-sheet melt, which is not included in the present simulations, is expected to amplify the long-term impact of anthropogenic carbon storage on seawater chemistry by reducing AABW formation and deep ocean ventilation. The resulting increase in deep-water residence time would prolong the isolation of carbon-rich waters from the atmosphere, reinforcing deep-ocean acidification (Lago and England, 2019; Chen et al., 2023).

4.1.4 Ocean circulation controls on biological productivity

Simulation results for detritus export at the depth of 100 m vary substantially across emission scenarios. Detritus export is projected to increase over the long term under SSP1–2.6 and SSP2–4.5, but to decrease under SSP3–7.0 and SSP5–8.5 (Fig. 2a). Most of this divergence develops during the 21st century and closely follows the rapid adjustment of the large-scale ocean circulation. Detritus export declines as the weakening of the AMOC becomes more pronounced under stronger forcing; this transition occurs when $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ exceeds approximately 3 times the pre-industrial level (Fig. 3). Consequently, the long-term behavior of detritus export is largely established by the end of the 21st century and remains comparatively stable thereafter. This behavior is difficult to identify from end-of-century analyses alone. Detritus export decreases across the high-latitude North Atlantic under all scenarios. The largest decline occurrs under SSP5–8.5, in which the magnitude of AMOC weakening is also the largest (Fig. 2d).

Under SSP5–8.5, enhanced nutrient trapping within the AABW cell after 2250 leads to a sustained reduction in surface nutrient availability, particularly at low and mid-latitudes (Fig. A4). This nutrient sequestration in the deep ocean coincides with declining biological productivity in regions such as the equatorial Pacific, contributing to long-term reductions in biological export (Fig. 2d). On the basis of the available model diagnostics, it is unclear whether this low-latitude decrease is a consequence of the changes in the Southern Ocean (Moore et al., 2018) or is associated with local changes resulting from the interplay between water mass renewal and remineralization under large perturbations (Rodgers et al., 2024). Attribution of the cause of these low-latitude changes will require model sensitivity studies in the future and is beyond the scope of this study.

4.1.5 Ecological implications of multiple stressors

We found pronounced contrasts across ecosystem stress indicators and emission scenarios, reflecting interactions between physical, chemical, and biogeochemical changes and cumulative effects. Considering that Ω_arag and Ω_calc integrate the effects of anthropogenic CO₂ uptake, we have included them as indicators of ecosystem stress in this study. Their values can indicate increased energetic costs and reduced fitness in calcifying organisms and have often been used as ecologically relevant thresholds.

Under SSP1–2.6, ecosystem stress remains comparatively limited and several indicators gradually recover toward present-day conditions (Fig. A11a–c). However, recovery is slow and spatially heterogeneous, with polar and subpolar ecosystems remaining close to critical thresholds for centuries (Fig. 10). Persistent anomalies in temperature, nutrient availability, oxygen, and alkalinity indicate that ecological recovery may substantially lag behind the stabilization of atmospheric forcing.

As illustrated by the results from the intermediate-emission scenario SSP2–4.5, ecosystem responses are not determined solely by the magnitude of atmospheric forcing. Under SSP2–4.5, warming, ocean acidification, and altered nutrient distributions continue to cause ecosystem stress to intensify although the increase in $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ slows after 2100 (Fig. 10). These results demonstrate that marine ecosystem stress depends on the persistence of ocean biogeochemical change as well as on the trajectory of $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$.

Under SSP5–8.5, ecosystem stress becomes widespread throughout nearly all ocean basins (Fig. A11g–i). Consistent with previous 21st-century projections (e.g., Kwiatkowski et al., 2020), stress first emerges in tropical and high-latitude regions, but our extended simulations reveal continued expansion and intensification well beyond 2100. The coexistence of persistent warming, declining nutrient availability, widespread ocean acidification, and prolonged carbonate mineral undersaturation indicates a transition toward ecosystem states with multiple chronic stressors and limited potential for recovery.

The extended simulations also reveal distinct long-term behaviors of Ω_arag and Ω_calc (Fig. A16). Because calcite generally remains supersaturated at lower saturation states than aragonite, the emergence of calcite undersaturation represents a more severe perturbation of seawater carbonate chemistry. Under SSP1–2.6, aragonite undersaturation develops transiently in polar regions (Fig. A16a–d), whereas calcite remains largely supersaturated (Fig. A16i–l). By contrast, under SSP5–8.5, calcite undersaturation expands into subpolar regions (Fig. A16m–p), indicating that high-emission pathways may push marine ecosystems beyond chemical thresholds associated with substantially lower resilience.

Together, these results demonstrate that emission pathways determine not only the magnitude but also the persistence and reversibility of marine ecosystem stress. Strong emission mitigation substantially limits cumulative impacts but does not eliminate the long-term ecological legacy of anthropogenic forcing.

These findings also have implications for marine ecosystem impact assessments. Outputs from ESMs are widely used as environmental forcing for higher-trophic-level ecosystem models. For example, the CMIP6-driven FishMIP simulations project substantial declines in global marine animal biomass together with pronounced regional redistributions of marine ecosystems under future climate change (Tittensor et al., 2021). Our results suggest that ecosystem responses may continue to evolve for centuries after atmospheric forcing has stabilized because key stressors, including ocean acidification, carbonate mineral saturation, nutrient availability, and warming, continue to change beyond the 21st century. Extending ESM simulations beyond the conventional 2100 time horizon may therefore improve assessments of long-term marine ecosystem responses and their implications for fisheries and ecosystem services.

4.2 Implications for future ocean carbon uptake

Here, we synthesize the implications of the multi-centennial ocean responses that we have identified in this study for the long-term behavior of the climate–carbon cycle system. Our results demonstrate that even after the stabilization of $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ and surface temperatures around the mid–22nd century, the ocean continues to absorb, redistribute, and store carbon over subsequent centuries. The deep ocean therefore does not reach a passive equilibrium once surface conditions stabilize, but remains an active component of the carbon cycle because of ongoing circulation and biogeochemical processes.

In particular, DIC accumulating along AABW pathways is physically isolated from the surface and retained in the ocean interior over long timescales, contributing to progressive deep-ocean acidification. This transport and storage of carbon persists on millennial timescales (Yamamoto et al., 2018; Hajima et al., 2020b), providing a clear illustration of emission legacy – the long-term impacts of past and present greenhouse gas emissions that endure even after atmospheric forcing stabilizes. Such behavior reflects the inertia of the ocean and the broader Earth system, underscoring the potential for long-lasting changes in deep-ocean chemistry and climate.

With respect to ocean carbon chemistry, our analysis of buffering capacity indicates that the ocean approaches a limit in its ability to absorb additional CO₂. Although these simulations are concentration-driven and $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ is prescribed, the emergence of a near-minimum buffering regime implies that, under emission-driven scenarios, declining oceanic uptake efficiency could substantially amplify climate–carbon cycle feedbacks.

By extending simulations beyond the conventional cutoff of 2100 and evaluating ocean responses under scenarios extended to 2500, this study expands the temporal scope of future projections. A robust understanding of ocean response timescales is essential for assessing long-term Earth system stability, particularly in the context of climate mitigation strategies, emission legacy, and the potential deployment of negative emissions technologies.

4.3 Future challenges and uncertainties

In this study, we have presented simulation results that provide new insight into multi-centennial ocean biogeochemical responses. Because of limitations associated with model structure, experimental design, and the representation of key Earth-system processes, several uncertainties remain. Addressing these uncertainties will be essential for improving the robustness of long-term climate projections and the application of long-term projections to marine ecosystem and climate-impact assessments.

4.3.1 Model structural uncertainty

A major source of uncertainty arises from structural differences among ESMs, especially in terms of model representations of the carbon cycle and marine biogeochemistry. Previous multi-model analyses have demonstrated substantial spread in projections of NPP, biological export, and ocean acidification, reflecting differences in nutrient limitation, ecosystem parameterizations, biological pump formulations, and carbonate chemistry (Bopp et al., 2013; Kwiatkowski et al., 2020; Bopp et al., 2022).

The simulations presented in this study show contrasting responses of detritus export across emission scenarios, which can largely be explained by the combined influence of atmospheric CO₂ forcing and changes in large-scale ocean circulation (Fig. 3). Under high-emission scenarios, weakening of the AMOC reduces nutrient resupply to the North Atlantic and suppresses detritus export, whereas detritus export is maintained or enhanced under lower-emission scenarios because the overturning circulation remains comparatively stable.

However, CMIP-class ESMs vary in terms of their simulations of the mean state, and deep convection of the overturning circulation (Weijer et al., 2020; Masson-Delmotte et al., 2021). These uncertainties and those associated with the magnitude and timing of AMOC weakening propagate through to projections of nutrient supply, marine productivity, and carbon export. Accordingly, uncertainties in biological export projections are associated with those in the simulated physical climate response as well as model parameterizations of the ecosystem.

The concentration-driven experimental design also represents an important source of uncertainty. Because $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ is prescribed, feedbacks between anthropogenic emissions, the coupled climate system, and the carbon cycle are not represented explicitly. Previous comparisons between concentration- and emission-driven experiments have shown differences in ocean carbon uptake trajectories and stronger climate–carbon feedbacks in emission-driven simulations (Friedlingstein et al., 2014, 2025). In this study, we have identified future declines in buffering capacity and ocean carbon uptake, which may have larger magnitudes under fully coupled emission-driven scenarios.

4.3.2 Missing Earth-system processes

Several processes that become increasingly important on centennial-to-millennial timescales, such as interactive ice-sheet dynamics, carbonate dissolution in marine sediments and coral reefs, and associated alkalinity feedbacks, are not represented in the MIROC-ES2L model framework used in the present study.

Freshwater input from ice sheet melt has long been recognized as a key driver of large-scale ocean circulation changes, particularly through its influence on buoyancy-driven overturning (Rahmstorf, 2002). Recent studies have shown that Antarctic ice sheet melt can weaken AABW formation, alter Southern Ocean circulation, and substantially modify ocean carbon uptake and storage (Bronselaer et al., 2018, 2020), whereas freshwater perturbations from the Greenland Ice Sheet have been shown to reduce ocean carbon uptake through their influence on the AMOC (Swingedouw et al., 2007). The omission of dynamic freshwater forcing in the present study therefore represents an important source of uncertainty in projections of long-term ocean carbon storage and acidification.

Additionally, the magnitudes of the reductions in surface-ocean alkalinity and buffering capacity may have been overestimated in our simulations, which did not take into account sediment carbonate dissolution declining coral reef calcification, and benthic alkalinity feedbacks (Kwiatkowski et al., 2025). Sediment dissolution supplies alkalinity to the ocean, partially compensates for ocean acidification, and weakens the long-term decline in carbonate buffering capacity. On longer timescales, carbonate compensation provides one of the principal negative feedbacks regulating atmospheric CO₂ and restoring ocean carbonate chemistry (Archer et al., 1998; Archer and Brovkin, 2008).

Furthermore, continental shelf and coastal oceans and are strongly influenced by terrestrial carbon and nutrient inputs but are poorly resolved in coarse-resolution global models (Yamamoto et al., 2022). River-derived carbon may drive oceanic CO₂ outgassing that is equivalent to approximately 20 % of present anthropogenic ocean carbon uptake (Lacroix et al., 2020; Regnier et al., 2022). Therefore, the biogeochemical responses in the continental shelf and coastal oceans represent another important source of uncertainty in the long-term projections reported in this study.

4.3.3 Future directions

Future modeling efforts should prioritize improving the coupling between ice sheets and the ocean, as well as targeted sensitivity experiments that explore the impact of freshwater perturbations on long-term carbon cycle dynamics. Improved land–ocean coupling and higher-resolution global simulations are also needed. Beyond these physical and Earth system interactions, further progress will require improved representation of biogeochemical and ecological processes.

At present, ecosystem responses in ESMs are commonly assessed using aggregated physical and biogeochemical stress indicators, rather than explicit representations of species- or community-level dynamics. Such simplified representations limit the ability of current models to capture ecosystem complexity, biodiversity responses, and trophic interactions, particularly at higher trophic levels (Stock et al., 2017; Lotze et al., 2019). Recent studies incorporating multiple plankton functional types demonstrate that accounting for greater ecosystem complexity can substantially alter projected community composition and functional structure under climate change (Henson et al., 2022). Addressing these limitations will require the integration of ecosystem models that explicitly represent multiple plankton functional types and higher trophic levels, together with sustained observational constraints.

In addition to ecosystem complexity, seasonal variability represents another important dimension of future ecological change. Although this study focuses on annual mean states and long-term trends, changes in the timing and amplitude of seasonal cycles can substantially alter the seasonal exposure of marine organisms to acidification stress (Kwiatkowski and Orr, 2018; Orr et al., 2022).

Another important source of uncertainty concerns the strength and nonlinearity of climate–carbon feedbacks associated with ocean carbonate chemistry. Although projections of ocean acidification are relatively robust, substantial uncertainty still remains with respect to the timing and magnitude of changes in buffering capacity and the resulting limits on long-term ocean carbon uptake (Orr et al., 2005; Egleston et al., 2010). Biases in the simulated Revelle factor have also been linked to biases in projected ocean carbon uptake across CMIP6 models (Terhaar et al., 2022).

Future progress will require coordinated multi-model experiments that extend beyond the conventional 2100 horizon while incorporating improved representations of ocean circulation, carbonate chemistry, sediment processes, ice-sheet interactions, and marine ecosystem dynamics.

5 Conclusions

In this study, we presented comprehensive analyses of ESM simulations from MIROC-ES2L and examined the sensitivity of ocean biogeochemistry and ocean biogeochemical fluxes under a range of anthropogenic emission scenarios extended to the year 2500. This timescale is much longer than the timescale that is typically used in climate projection studies. Our results confirm that changes can persist over multiple centuries, with broad implications for human resource management and ecosystem health. Importantly, our results underscore the value of choosing a low-emission strong-mitigation pathway (rather than a high-emission weak-mitigation pathway) in promoting both ocean health and limiting marine ecosystem stress over these extended timescales. Emission mitigation substantially reduces many of the deleterious impacts of the forcing produced by high emissions, including warming, deoxygenation, and reduced availability of nutrients in the euphotic zone.

In terms of the carbon cycle, the ocean is projected to continue to absorb CO₂ from the atmosphere beyond the 21st century, leading to sustained accumulation of DIC in the deep ocean. This is especially evident along pathways of AABW, where acidification in the deep ocean also intensifies. Analysis of the buffering capacity reveals a substantial weakening of the ocean's carbonate buffering system under SSP5–8.5, indicating a decline in the efficiency of oceanic CO₂ uptake in the future. This weakening is delayed in the Southern Ocean because of the increase in biological export and influence of deep water masses.

Biological export exhibits a pronounced dependence on the trajectory of $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$. Detritus export increases under low-emission scenarios but declines under sustained high-emission scenarios, with the transition occurring when $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ reaches approximately 3 times the preindustrial level. Although this threshold behavior reflects the combined influence of multiple physical and biogeochemical processes, it is closely associated with changes in large-scale ocean circulation, particularly the weakening of the AMOC and the resulting redistribution of nutrients.

Ecosystem stressors – warming, acidification, and deoxygenation – are projected to expand in both space and time. Even under the low-emission scenario SSP1–2.6, return to modern-day levels is delayed by roughly 100 years, whereas under SSP2–4.5, more areas are stressed and ecosystem stress is projected to intensify through the 25th century. In the Southern Ocean, the aragonite saturation state is projected to fall below 1.0 by 2100 under SSP5–8.5, posing serious risks to calcifying organisms. Furthermore, reduced deep-ocean ventilation in the North Atlantic, North Pacific, and Southern Ocean is expected to drive oxygen loss outside current OMZs.

These findings underscore the importance of accounting for sustained and delayed ocean responses in climate projections beyond the 21st century. Long-term simulations extending to the year 2500 show that persistent anthropogenic forcing leads to profound and long-lasting alterations in ocean biogeochemistry, particularly through changes in the polar regions and the deep ocean. Weakening and reorganization of the overturning circulation result in the redistribution of biogeochemical tracers over multi-centennial timescales, reducing the ocean's capacity to buffer $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$ and limiting future carbon uptake. As a consequence, deep-ocean acidification intensifies, especially along AABW pathways, raising concerns about persistent impacts on benthic ecosystems and highlighting the need to place high priority on these processes in future model development and climate policy planning.

These findings highlight long-term (multi-centennial) memory of the ocean interior of climate perturbations and have a number of broad implications. One of them is the reversibility of such changes under climate stabilization. Because perturbations persist in the deep ocean under high-emission scenarios, identification of the timescales of the persistence of marine biogeochemical perturbations should be a priority in future studies. Fully understanding reversibility will likely require developing attribution methods that can deconvolve the direct effects of warming on biogeochemical cycling and the indirect effects as manifested through changes in circulation. This will provide insight into the indirect contributions of weakening overturning circulation structures to global biogeochemical change and facilitate assessments of the relative contribution of the direct effect of warming on biogeochemical processes beyond 2100 or 2150 and its temporal variability.

Appendix A

Figure A1 Atmospheric pCO₂ forcings prescribed by historical and extended SSP scenarios.

Figure A2 Surface DIC follows atmospheric pCO₂, whereas deep-ocean DIC continues to increase gradually over centuries. Depth-integrated DIC concentrations (unit: mmol m⁻³) in three ocean layers (0–200, 200–1000, and 1000–6000 m), simulated under the extended SSP scenarios.

Figure A3 Under SSP5–8.5, the AMOC becomes much weaker and shallower by the end of the 21st century. (a) Time series of AMOC and Pacific Meridional Overturning Circulation (PMOC) strength under SSP1–2.6, SSP2–4.5, and SSP5–8.5. AMOC and PMOC strengths are defined as the maximum meridional overturning streamfunction over 31–90° N and 500–5400 m depth in the Atlantic and 30° S–90° N and 1000–5400 m depth in the Pacific, respectively. Meridional overturning streamfunctions in the global, Atlantic, and Pacific Oceans for (b) SSP1–2.6 and (c) SSP5–8.5.

Figure A4 Nutrients accumulate in AABW cells, whereas subsurface nutrients decrease in low latitudes. Zonal mean vertical sections of nitrate concentration anomalies (unit: mmol m⁻³) under SSP5–8.5. From left to right, the panels show longitudinally averaged sections for the Atlantic, Southern, and Pacific Oceans. From top to bottom, rows represent anomalies for the (a) 2040s, (b) 2090s, (c) 2240s, and (d) 2490s, relative to the 2010s average under SSP1–2.6.

Figure A5 Ocean temperatures at the surface stabilize, whereas temperatures below the depth of 200 m continue to rise over centuries across all scenarios. Depth-integrated temperature (unit: °C) in three ocean layers (0–200, 200–1000, and 1000–6000 m), simulated under the extended SSP scenarios.

Figure A6 Surface-ocean pCO₂ broadly tracks atmospheric pCO₂, despite regional air–sea differences. Annual-mean changes in pCO₂ in the surface ocean (unit: ppm) by latitude band, simulated under (a) SSP1–2.6 and (b) SSP5–8.5. Colored lines represent different latitude bands as follows: gray (global mean), dark blue (equatorial: 15° S–15° N), orange (northern subtropics: 15–40° N), yellow (southern subtropics: 15–45° S), pink (northern subpolar: 40–60° N), green (Arctic: >60° N), and light green (Antarctic: <45° S).

Figure A7 The Southern Hemisphere westerlies strengthen under SSP5–8.5. 10 m wind speed anomalies (unit: m s⁻¹) under SSP5–8.5, relative to the 2010s mean under SSP1–2.6. Panels show the (b) 2040s, (c) 2090s, (d) 2240s, and (e) 2490s.

Figure A8 Surface DIC is the primary driver of pCO₂ changes, but under SSP5–8.5 in the mid–22nd century, alkalinity also makes a secondary yet notable contribution. Regional contributions of temperature (red), salinity (blue), DIC (yellow), and alkalinity (green) to the temporal change of surface $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ocn}}$ (black) under (a) SSP1–2.6 and (b) SSP5–8.5. Thick lines represent 20-year running means, whereas thin lines indicate 5-year mean values. The horizontal axis represents time (years), and the vertical axis shows the time derivative of pCO₂ ($\mathrm{d}p{\mathrm{CO}}_{\mathrm{2}}/\mathrm{d}t$, unit: ppm yr⁻¹).

Figure A9 The vertical gradient in alkalinity between the surface and deep ocean intensifies. Zonal mean vertical sections of alkalinity concentration anomalies (unit: mmol m⁻³) under SSP5–8.5. From left to right, the panels show longitudinally averaged sections for the Atlantic, Southern, and Pacific Oceans. From top to bottom, rows represent anomalies for the (a) 2040s, (b) 2090s, (c) 2240s, and (d) 2490s, relative to the 2010s average under SSP1–2.6.

Figure A10 Deep-ocean DIC continues to increase over centuries along the AABW pathway. Zonal mean vertical sections of DIC concentration anomalies (unit: mmol m⁻³) under SSP5–8.5. From left to right, the panels show longitudinally averaged sections for the Atlantic, Southern, and Pacific Oceans. From top to bottom, rows represent anomalies for the (a) 2040s, (b) 2090s, (c) 2240s, and (d) 2490s, relative to the 2010s average under SSP1–2.6.

Figure A11 Even under the low-emission SSP1–2.6 scenario, it takes until around the year 2200 for surface ocean temperature, acidification, deoxygenation, and aragonite saturation state to return to near modern levels. Projected time series of open-ocean temperature, pH, DIC, O₂, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, and Ω_aragonite under SSP1–2.6, SSP2–4.5, SSP3–7.0, and SSP5–8.5. Results are shown for four depth ranges (0–200, 200–1000, 1000–2000, and >2000 m).

Figure A12 Surface-ocean nitrate stabilizes and slowly recovers under low emissions, while deep-ocean nitrate continues to accumulate. Depth-integrated nitrate concentrations (unit: mmol m⁻³) in three ocean layers (0–200, 200–1000, and 1000–6000 m), simulated under the extended SSP scenarios.

Figure A13 Deep-ocean oxygen declines persistently across emission scenarios. Depth-integrated oxygen concentrations (unit: mmol m⁻³) in three ocean layers (0–200, 200–1000, and 1000–6000 m), simulated under the extended SSP scenarios.

Figure A14 Surface-ocean alkalinity declines under climate forcing, while deep-ocean alkalinity continues to increase. Depth-integrated alkalinity concentrations (unit: mmol m⁻³) in three ocean layers (0–200, 200–1000, and 1000–6000 m), simulated under the extended SSP scenarios.

Figure A15 Subsurface oxygen concentrations decrease substantially at high latitudes, whereas changes within equatorial OMZs are small. Dissolved oxygen anomalies (unit: mmol m⁻³) under SSP5–8.5, relative to the 2010s mean of SSP1–2.6. Panels show the (b) 2040s, (c) 2090s, (d) 2240s, and (e) 2490s.

Figure A16 Under SSP5–8.5, the aragonite saturation state falls below 1 after the 21st century, except in parts of the subtropics. Aragonite saturation state (Ω_arag) under (a–d) SSP1–2.6 and (e–h) SSP5–8.5 for the upper ocean (0–200 m) in the (a, e) 2040s, (b, f) 2090s, (c, g) 2240s, and (d, h) 2490s. Calcite saturation state (Ω_calc) under (i–l) SSP1–2.6 and (m–p) SSP5–8.5 for the upper ocean (0–200 m) in the (i, m) 2040s, (j, n) 2090s, (k, o) 2240s, and (l, p) 2490s.

Figure A17 Under SSP5–8.5, summer sea ice reaches a minimum by the mid–21st century, and winter sea ice does so within the 22nd century. Changes in sea-ice concentration by latitude band under four scenarios: (a) SSP1–2.6 (February–April; FMA), (b) SSP1–2.6 (August–October; ASO), (c) SSP5–8.5 (February–April; FMA), and (d) SSP5–8.5 (August–October; ASO). Colored lines represent different latitude bands as follows: gray (global mean), dark blue (equatorial: 15° S–15° N), orange (northern subtropics: 15–40° N), yellow (southern subtropics: 15–45° S), pink (northern subpolar: 40–60° N), green (Arctic: >60° N), and light green (Antarctic: <45° S).

Appendix B: Multi-centennial ocean physical responses

Biogeochemical changes are strongly modulated by the physical state of the ocean, not only through source and sink processes but also via advection and diffusion, which transport biogeochemical tracers throughout the ocean interior. Therefore, in this section, we provide a brief overview of the key physical responses of deep ocean circulation, temperature, salinity, stratification, and sea-ice, focusing on aspects that are the most relevant to the interpretation of biogeochemical changes.

Global-scale ocean heat uptake and enhanced stratification lead to a pronounced weakening of the AMOC (Fig. A3). Under high-emission scenarios in particular, NADW formation declines substantially by the year 2100, with the associated deep ocean volume transport decreasing from approximately 30 Sv (where $\mathrm{1}\mathrm{Sv}={\mathrm{10}}^{\mathrm{6}}{\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$) in the preindustrial era to about 5 Sv. This pronounced decline occurs primarily within the first 100 years of the simulation. In the Pacific Ocean, the deep circulation driven by AABW also weakens considerably, with transport rates approaching zero as global temperature reaches its maximum. As a result, the structure of the deep ocean circulation is substantially different from that of the present day.

With increasing $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$, SST rises accordingly, and tends to remain elevated once the CO₂ concentration stabilizes (Fig. A5). Under SSP5–8.5, the global mean SST increases by approximately +2 °C by 2100 and reaches around +4 °C by 2200. This warming exhibits substantial regional variability; SST increases more rapidly in the Arctic. Warming progresses more slowly in the Southern Ocean (south of 60° S) and the subantarctic region (30–60° S), where the rate of increase is roughly half that observed in the Arctic.

In contrast to the surface ocean, deep ocean temperatures continue to increase beyond 2300 under SSP5–8.5 approximately linearly (Fig. A5). This deep ocean warming is particularly pronounced in the Arctic, where the rate of temperature increase is roughly 3 times that in other regions. By 2500, the temperature increase below the depth of 1000 m in the Arctic Ocean exceeds 7 °C, compared to a global mean increase of approximately 2 °C at the same depth. Such regional differences are likely attributable to variations in the propagation of heat into the deep ocean and differences in deep ocean circulation. Under the low-emission scenario SSP1–2.6, the temperature increase in the upper 200 m of the ocean peaks in the latter half of the 21st century, reaching approximately +0.8 °C. By contrast, below 1000 m, warming continues gradually and reaches a comparable magnitude of about +0.8 °C by 2500.

Moreover, there is significant warming in the mid-depth layer (approximately 500–1000 m) even under low-emission scenarios, and this warming persists through to the year 2500. As a result, although surface warming enhances stratification, the subsurface is influenced by warming from intermediate depths, potentially weakening the stratification beneath the surface layer in some regions.

Timescales of sea-ice change are relatively short compared with those of changes of other physical variables. Under SSP5–8.5, wintertime sea-ice concentration reaches its minimum in the Arctic Ocean by the mid–22nd century and in the Southern Ocean around 2200 (Fig. A17). By contrast, summertime sea ice nearly disappears by the late 21st century.

The projected future loss of sea ice is expected to further amplify surface warming in the polar regions through reduced albedo. It is expected to enhance surface stratification and potentially impact biogeochemical cycles through changes in both the light environment and vertical nutrient supply.