Arctic Ocean acidification over the 21st century co-driven by anthropogenic carbon increases and freshening in the CMIP6 model ensemble

The uptake of anthropogenic carbon (C_ant) by the ocean leads to ocean acidification, causing the reduction of pH and the saturation states of aragonite (Ω_arag) and calcite (Ω_calc). The Arctic Ocean is particularly vulnerable to ocean acidification due to its naturally low pH and saturation states and due to ongoing freshening and the concurrent reduction in total alkalinity in this region. Here, we analyse ocean acidification in the Arctic Ocean over the 21st century across 14 Earth system models (ESMs) from the latest Coupled Model Intercomparison Project Phase 6 (CMIP6). Compared to the previous model generation (CMIP5), models generally better simulate maximum sea surface densities in the Arctic Ocean and consequently the transport of C_ant into the Arctic Ocean interior, with simulated historical increases in C_ant in improved agreement with observational products. Moreover, in CMIP6 the inter-model uncertainty of projected changes over the 21st century in Arctic Ocean Ω_arag and Ω_calc averaged over the upper 1000 m is reduced by 44–64 %. The strong reduction in projection uncertainties of Ω_arag and Ω_calc can be attributed to compensation between C_ant uptake and total alkalinity reduction in the latest models. Specifically, ESMs with a large increase in Arctic Ocean C_ant over the 21st century tend to simulate a relatively weak concurrent freshening and alkalinity reduction, while ESMs with a small increase in C_ant simulate a relatively strong freshening and concurrent total alkalinity reduction. Although both mechanisms contribute to Arctic Ocean acidification over the 21st century, the increase in C_ant remains the dominant driver. Even under the low-emissions Shared Socioeconomic Pathway 1-2.6 (SSP1-2.6), basin-wide averaged Ω_arag undersaturation in the upper 1000 m occurs before the end of the century. While under the high-emissions pathway SSP5-8.5, the Arctic Ocean mesopelagic is projected to even become undersaturated with respect to calcite. An emergent constraint identified in CMIP5 which relates present-day maximum sea surface densities in the Arctic Ocean to the projected end-of-century Arctic Ocean C_ant inventory is found to generally hold in CMIP6. However, a coincident constraint on Arctic declines in Ω_arag and Ω_calc is not apparent in the new generation of models. This is due to both the reduction in Ω_arag and Ω_calc projection uncertainty and the weaker direct relationship between projected changes in Arctic Ocean C_ant and changes in Ω_arag and Ω_calc.

1 Introduction

Human activities such as the burning of fossil fuels, cement production, and land use change have released large amounts of carbon into the atmosphere that cause global warming. The ocean mitigates global warming by taking up around one-quarter of this anthropogenic carbon (C_ant) (Friedlingstein et al., 2019). However, the increase of carbon in the ocean causes ocean acidification, a process that decreases pH, carbonate ion (${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$) concentrations, and in consequence the calcium carbonate (CaCO₃) saturation states of calcite and aragonite minerals (Haugan and Drange, 1996; Orr et al., 2005). The Arctic Ocean is particularly vulnerable to ocean acidification due to its naturally high dissolved inorganic carbon concentrations, its low carbonate ion concentrations, and its thus naturally low saturation states (Orr et al., 2005; Fabry et al., 2009; Gattuso and Hansson, 2011; Riebesell et al., 2013; AMAP, 2018).

In contrast to most of the global ocean, Arctic Ocean acidification is caused not solely by increasing C_ant concentrations (Anderson et al., 2010; Ulfsbo et al., 2018; Terhaar et al., 2020a, b) but also by freshening (Koenigk et al., 2013; Nummelin et al., 2016; Shu et al., 2018; Brown et al., 2020; Woosley and Millero, 2020). Fresh water from rivers, precipitation, and sea ice typically has much lower total alkalinity (A_T) and total dissolved inorganic carbon (C_T) concentrations than the ocean and therefore, in the absence of indirect impacts on other fluxes, dilutes both marine A_T and C_T (Xue and Cai, 2020). As freshwater A_T and C_T concentrations are generally similar, freshwater fluxes into the ocean typically act to reduce the difference between A_T and C_T, decreasing marine ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ concentrations and ocean pH (Bates et al., 2009; Bates and Mathis, 2009; Yamamoto-Kawai et al., 2011; Waldbusser and Salisbury, 2014; Wanninkhof et al., 2015; Xue and Cai, 2020). In the Arctic Ocean, projected freshening over the 21st century is larger than in most other ocean regions due to ongoing sea ice melt, positive precipitation minus evaporation, and large river runoff (Rawlins et al., 2010; Rudels, 2015; Shu et al., 2018).

Due to freshening and increasing C_ant concentrations, the Arctic Ocean is projected to be the first large-scale ocean region to become undersaturated with respect to the metastable CaCO₃ polymorph aragonite (Ω_arag<1) (Steinacher et al., 2009). Under the Representative Concentration Pathway 8.5 (RCP8.5) high-emissions scenario, Arctic Ocean mesopelagic waters may even become undersaturated with respect to the more stable CaCO₃ polymorph calcite (Ω_calc<1) before 2100 (Terhaar et al., 2020a). Aragonite and calcite undersaturation is likely to affect the growth, reproduction, and survival of calcifying organisms, such as sea butterflies (Comeau et al., 2010) and foraminifera (Davis et al., 2017), and could have ramifications for the wider Arctic ecosystem (Armstrong et al., 2005; Karnovsky et al., 2008), including some of its most iconic predators, such as grey whales and walruses (Jay et al., 2011; AMAP, 2018).

Projections of Arctic Ocean acidification over the 21st century had considerable subsurface uncertainties in the simulations conducted as part of the Coupled Model Intercomparison Project Phase 5 (CMIP5) (Steiner et al., 2013; Terhaar et al., 2020a), with projected end-of-century basin-wide Ω_arag in mesopelagic waters ranging from 0.61 to 1.05. This large uncertainty has been attributed to multiple factors, including variable inflow of Atlantic waters and their subsequent subduction; difficulties resolving the narrow passages between the Arctic Ocean and its surrounding basins; and differences in brine rejection during sea ice formation, which is critical to the formation of dense Arctic waters (Terhaar et al., 2019b).

1.1 Emergent constraints on Arctic Ocean carbon uptake and acidification

Emergent constraints are a suite of statistical techniques that relate observable trends or sensitivities across multi-model ensembles to differences in model projections in order to reduce future uncertainties (Allen and Ingram, 2002; Hall et al., 2019). An elegant early application that demonstrates many of the principles of emergent constraints is that of Hall and Qu (2006). In their study of models that contributed to the Intergovernmental Panel on Climate Change's Fourth Assessment Report (AR4), Hall and Qu found a strong positive correlation between the magnitude of a model's snow albedo feedback on present-day seasonal timescales and under future climate change. They concluded that relevant model biases were consistent across these contrasting timescales and therefore observations of the seasonal snow albedo feedback could be used to constrain the ensemble range of the projected snow albedo feedback under climate change. Since the publication of Hall and Qu (2006), emergent constraint approaches have been applied extensively within the Earth sciences to constrain, amongst other things, projections of climate sensitivity (Caldwell et al., 2018), Arctic sea ice extent (Boé et al., 2009), precipitation extremes (O'Gorman, 2012; DeAngelis et al., 2016), carbon cycle feedbacks (Cox et al., 2013; Wenzel et al., 2014), and marine primary production (Kwiatkowski et al., 2017).

Recently Terhaar et al. (2020a) showed that an emergent constraint could be applied to CMIP5 projections of the Arctic Ocean C_ant inventory and coincident acidification over the 21st century. As the C_ant increase in the Arctic Ocean is mainly driven by the inflow of C_ant-rich waters from the Atlantic and their subsequent subduction in the Barents Sea (Midttun, 1985; Rudels et al., 1994, 2000; Jeansson et al., 2011; Smedsrud et al., 2013), the capability of each model to form dense surface waters in the Barents Sea was shown to strongly influence the future Arctic Ocean C_ant inventory. By constraining simulated surface water densities with observations, uncertainties related to the end-of-century Arctic Ocean C_ant inventory in 2100 were reduced by around one-third, and the best estimate under RCP8.5 was increased by 20 % to 9.0±1.6 Pg C (Terhaar et al., 2020a). Along with the projected C_ant inventory, uncertainties in the projected associated basin-wide Arctic Ocean acidification could also be reduced in CMIP5. It should be noted, however, that in CMIP5 projected freshening and reductions in alkalinity were of minor importance for Arctic Ocean acidification over the 21st century. Moreover, the models have been shown to underestimate historical freshwater fluxes (1992–2012) in the Arctic Ocean by around 50 % (Shu et al., 2018), which suggests they might also have underestimated freshwater fluxes over the 21st century.

Given that emergent constraints in many cases conflict with one another (Caldwell et al., 2018; Brient, 2020) and can even be derived from data-mined pseudocorrelations (Caldwell et al., 2014), it is critical to test published constraints, and the mechanisms that underpin them, across Earth system model (ESM) generations (Eyring et al., 2019; Hall et al., 2019). The CMIP6 simulations provide such an opportunity (Schlund et al., 2020).

1.2 From CMIP5 to CMIP6 models and simulations

During the transition from CMIP5 to CMIP6, ESMs have generally improved the simulation of ocean dynamics and marine biogeochemistry (Séférian et al., 2020). Across most ESMs, the horizontal and/or vertical resolution of ocean models has increased, which potentially has large effects on the representation of Arctic Ocean circulation, sea ice dynamics (Docquier et al., 2019), and the carbon cycle (Terhaar et al., 2019b). Ocean biogeochemical model components in CMIP6 also tend to have a more complex representation of the carbon and nutrient cycles than in CMIP5. In particular, the treatment of organic matter carbon cycling has generally evolved, with remineralization of particles in sediments now simulated in 10 out of 14 ESMs. These developments will likely have a large effect on simulating the Arctic Ocean biogeochemistry given that 50 % of the Arctic Ocean is made up of shelf seas (Jakobsson, 2002), where sedimentation and sediment remineralization are crucial components of the carbon and nutrient cycle (Brüchert et al., 2018; Grotheer et al., 2020). Furthermore, the external carbon and nutrient sources from glaciers, atmospheric deposition, and rivers are represented in more models in CMIP6 (Table 1) (Séférian et al., 2020). Riverine inputs in particular have been shown to be of importance for present-day Arctic Ocean acidification (Anderson et al., 2010; Tank et al., 2012) and its future changes (Terhaar et al., 2019a).

Ziehn et al. (2019a, b)(Ziehn et al., 2020)Swart et al. (2019a, b)(Swart et al., 2019e)Swart et al. (2019c, d)Swart et al., 2019eChristian et al., 2021Danabasoglu (2019a, b)(Danabasoglu et al., 2020)Danabasoglu (2019c, d)(Danabasoglu et al., 2020)Seferian (2018a, b)(Séférian et al., 2019)Guo et al. (2018a, b)Held et al., 2019Dunne et al., 2020aDunne et al., 2020bStock et al., 2020Krasting et al. (2018)John et al. (2018)Boucher et al. (2018a, b)(Boucher et al., 2020)(Hajima et al., 2020)Hajima et al. (2019)Tachiiri et al. (2019)Müller et al., 2018Mauritsen et al., 2019Schupfner et al. (2019)Jungclaus et al. (2019)Wieners et al. (2019a, b)(Mauritsen et al., 2019)Yukimoto et al. (2019b, c)(Yukimoto et al., 2019a)(Sellar et al., 2019)Tang et al. (2019)Good et al. (2019)

Table 1 The CMIP6 ESMs used in this study with their ocean–sea ice and marine-biogeochemical (MBG) model components plus their biogeochemical (BGC) riverine input.

^a Only SSP2-4.5 and SSP5-8.5. ^b Only SSP5-8.5.

In this study, we extend recent CMIP6 ocean biogeochemical assessments (e.g. Séférian et al., 2020; Kwiatkowski et al., 2020) and previous attempts to constrain projected Arctic Ocean C_ant uptake by

1. assessing projections of the Arctic Ocean C_ant inventory over the 21st century in CMIP6 simulations

2. exploring the role of C_ant inventory increases and freshening in driving concurrent basin-wide ocean acidification in the Arctic Basin

3. revaluating previous emergent constraints on the Arctic Ocean C_ant inventory and associated acidification using the CMIP6 model ensemble and multiple future-emissions scenarios.

2 Methods

2.1 Arctic Ocean

The Arctic Ocean was defined as the water north of the Fram Strait, the Barents Sea Opening, the Bering Strait, and the Baffin Bay following Bates and Mathis (2009). This is consistent with the previously published emergent constraint on projected Arctic Ocean C_ant and acidification (Terhaar et al., 2020a).

2.2 Earth system models

An ensemble of 14 ESMs from CMIP6 (Table 1) was used with one ensemble member per model. All models follow the biogeochemical protocols outlined in Orr et al. (2017). Riverine input of C_T and A_T is included in six ESMs (Table 1). The absence of C_T and A_T in riverine freshwater input causes an overly strong reduction of ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ concentrations and thus low-biased saturation states in coastal regions but is of minor importance on the pan-Arctic scale (Terhaar et al., 2019a). The spin-up length for each model varied between 500 years (IPSL-CM6A-LR) and 12 000 years (MPI-ESM1-2-LR) (Séférian et al., 2020).

For each model, monthly 3D fields of dissolved inorganic carbon, total alkalinity, dissolved inorganic phosphorus and silicon, temperature, and salinity were used. All 3D fields were regridded to the regular $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid with 33 depth levels used in the GLobal Ocean Data Analysis Project Version 2 (GLODAPv2) observational product (Lauvset et al., 2016) to add simulated changes of these variables over the 21st century to observations of the present-day mean state (see below).

C_ant was defined as the difference between annual means of dissolved inorganic carbon in the historical (1850–2014) simulations merged with the respective shared socioeconomic pathway (SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5; 2015–2100; Riahi et al., 2017) and the concurrent pre-industrial control simulations of each model. Output from 13 models was available for SSP1-2.6, output from 12 models for SSP2-4.5 and SSP3-7.0, and output from 14 models for SSP5-8.5 (Table 1).

Changes in A_T over the 21st century were calculated by subtracting changes in the pre-industrial control simulations from changes in the respective SSP. To quantify the effect of freshening on changes in A_T, the A_T anomalies for each model were further decomposed into changes resulting from freshening and from the combined effect of other biogeochemical processes by calculating the temporal evolution of salinity-corrected alkalinity with a reference salinity of 35 following Lovenduski et al. (2007). A zero-alkalinity endmember was assumed for fresh water. This assumption is correct for models with no alkalinity in fresh water but an overestimation for models with finite alkalinity concentrations in freshwater inputs (Table 1). Unfortunately, information on alkalinity concentrations in fresh water is not available for all models. Moreover, with the available model output, we cannot quantify the individual contributions of land ice melt, sea ice melt, precipitation minus evaporation, and riverine input to freshwater changes. Thus, for simplicity a zero-alkalinity endmember was assumed for all models.

Ocean carbon chemistry variables in ESMs commonly exhibit mean state biases (Orr et al., 2005; Steiner et al., 2013). Therefore, observations of C_T, A_T, dissolved inorganic phosphorus and silicon, temperature, and salinity from GLODAPv2 (Lauvset et al., 2016), which is normalized to the year 2002, were used to calculate present-day Ω_arag and Ω_calc, pH, and pCO₂ using the mocsy2.0 routine (Orr and Epitalon, 2015) and the equilibrium constants recommended for best practices (Dickson et al., 2007). Future Ω_arag and Ω_calc, pH, and pCO₂ were calculated for each model as the sum of the simulated changes in C_T, A_T, dissolved inorganic phosphorus and silicon, temperature, and salinity from 2002 onwards and the observed quantities in 2002.

The present-day maximum sea surface density in the Arctic Ocean was calculated from monthly climatologies over 1986–2005, constructed from the respective salinity and temperature outputs of each model. Maximum sea surface density was calculated, as in Terhaar et al. (2020a), as the mean density of the densest 5 % of Arctic surface waters (95th-percentile waters) over all 12 months of the year.

2.3 Simulations

The simulations performed within CMIP5 were not forced with the same atmospheric CO₂ concentrations as the simulations performed under CMIP6. In CMIP5, historically observed atmospheric CO₂ concentrations were used from 1850 to 2005 (Meinshausen et al., 2011). From 2006 onwards, the CO₂ concentrations follow the different RCPs. In CMIP6, the historical period was extended until 2014, and thereafter CO₂ concentrations follow the different SSPs (Meinshausen et al., 2019).

The different land and energy use assumptions in the SSPs (Riahi et al., 2017) compared to the RCPs (O'Neill et al., 2016) lead to higher atmospheric CO₂ trajectories over the 21st century for the Tier 1 SSPs (Meinshausen et al., 2019) compared to their RCP counterparts (Meinshausen et al., 2011), which results in globally greater surface ocean acidification in CMIP6 compared to CMIP5 (Kwiatkowski et al., 2020). Historical atmospheric CO₂ concentrations were also refined with additional data available since CMIP5 (Meinshausen et al., 2017). This refinement did not change the average atmospheric CO₂ concentration from 1850 to 2005 ($\mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}=\mathrm{0}\pm \mathrm{1}$ ppm) but did change annual CO₂ concentration for single years by up to ±2 ppm. Furthermore, global CO₂ concentrations were additionally provided as monthly latitudinally resolved concentrations, with model groups free to choose the forcing files they use (Meinshausen et al., 2017, 2019).

2.4 C_ant scaling

The different atmospheric CO₂ trajectories over the 21st century between CMIP5 and CMIP6 complicates a comparison of simulated Arctic Ocean C_ant inventories between model generations. To nevertheless compare the simulated Arctic Ocean C_ant, we used the commonly applied scaling approach that assumes that the change in marine C_ant is proportional to the atmospheric C_ant concentration (Gloor et al., 2003; Mikaloff Fletcher et al., 2006; Gerber et al., 2009; Gruber et al., 2009). Under this assumption, the C_ant inventory in 2100 for each scenario (RCP8.5, SSP1-2.6, SSP2-4.5, SSP3-7.0) is linearly rescaled to that of SSP5-8.5 by multiplying the simulated Arctic Ocean C_ant inventory under the respective scenario by the ratio of the mean atmospheric C_ant concentration from 1850 to 2100 in SSP5-8.5 and the respective scenario.

This approximation is likely imprecise when very different scenarios are compared, such as SSP1-2.6 and SSP5-8.5, as the effects of circulation changes, sea ice melt, and warming are not considered. However, when comparing scenarios with the same radiative forcing, such as SSP5-8.5 and RCP8.5, it permits a first-order comparison.

2.5 Observational constraints

As for the ESMs, the maximum sea surface density was calculated based on a monthly sea surface density climatology on a regular $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid, which was constructed from observed monthly salinity and temperature climatologies in the World Ocean Atlas 2018 (Locarnini et al., 2018; Zweng et al., 2018).

The density uncertainty was calculated from the temperature and salinity uncertainties that were reported by the World Ocean Atlas following standard propagation of uncertainty. The total uncertainty is a combination of (1) the standard deviations for sea surface density derived from published standard deviations of sea surface temperature and salinity for each grid cell and each month in the World Ocean Atlas 2018, and (2) the standard deviation from the weighted mean of the 95th-percentile-density waters (Terhaar et al., 2020a).

2.6 Emergent constraints and probability density functions of C_ant

To calculate the emergent constraint, first an ordinary least-squares regression was calculated between the simulated present-day maximum sea surface density and the Arctic Ocean C_ant inventory in 2100 for each ESM of the CMIP6 model ensemble. The uncertainty range was estimated using the 1σ prediction interval. In a second step the probability density function (PDF) from the observations was convoluted with the PDF from the linear regression, assuming a Gaussian distribution in both cases. The convolution of both PDFs is the constrained projection of the Arctic Ocean C_ant inventory following previous studies (Cox et al., 2013; Wenzel et al., 2014; Kwiatkowski et al., 2017). The PDFs for unconstrained projections of the C_ant inventory were calculated using equal weights for each model and assuming a Gaussian distribution.

Extending the analysis of Terhaar et al. (2020a), PDFs for the constrained projections of the C_ant inventory were calculated not only for the year 2100 but for each year from 2002 to 2100 and not only for the highest-emission scenario (RCP8.5) but for the four SSPs (SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5).

3 Results

3.1 Arctic Ocean C_ant inventory

3.1.1 Multi-model mean

Over the historical period from 1850 to 2005, the CMIP6 ESMs simulate an Arctic Ocean C_ant increase of 2.1±0.6 Pg C (inter-model standard deviation; Fig. 1a). Thus, the CMIP6 ESMs simulate an Arctic Ocean C_ant inventory in 2005 that is 36 % below the data-based estimate for the period from 1765 to 2005 of 3.3±0.3 Pg C (Terhaar et al., 2020b).

Figure 1 Arctic Ocean C_ant inventory over the 21st century. (a) Time series of multi-model mean Arctic Ocean C_ant inventory over the 21st century for SSP1-2.6 (blue), SSP2-4.5 (green), SSP3-7.0 (orange), and SSP5-8.5 (red) with the respective standard deviation across the model ensemble (shading) (n=12–14). The bars on the right side of panels indicate the inter-model standard deviation in the year 2100. (b) Time series of Arctic Ocean C_ant inventory for each scenario after the emergent constraint is applied.

Over the 21st century the Arctic Ocean C_ant inventory increases depending on the SSP. Following the low-emission pathway SSP1-2.6 leads to a projected Arctic Ocean C_ant inventory in 2100 of 6.2±1.3 Pg C. With increasing atmospheric C_ant concentrations from SSP1-2.6 to SSP5-8.5, the projected Arctic Ocean C_ant inventory in 2100 also increases, resulting in 8.0±1.7 for SSP2-4.5, 9.6±1.6 for SSP3-7.0, and 10.3±2.2 Pg C for SSP5-8.5. The relatively large uncertainties across the model ensemble (17–21 %) result in an overlap of the simulated C_ant inventories in 2100 for SSP2-4.5, SSP3-7.0, and SSP5-8.5 within ±1 SD.

3.1.2 Constrained results

As was shown for RCP8.5 in the CMIP5 model ensemble (Terhaar et al., 2020a), a linear relationship between maximum sea surface density and the Arctic Ocean C_ant inventory in 2100 is found across the CMIP6 model ensemble for all four of the SSPs (r²=0.63–0.72; Fig. 2a, c, e, and g). By deriving a similar relationship for the projected C_ant inventory in all years from 2000 to 2100 and combining this with observations of present-day sea surface density, the uncertainty of the projected Arctic Ocean C_ant inventory can be reduced throughout the 21st century (Fig. 1b). The emergent relationship is significant over all years from 2015 to 2100 (p<0.05) and the r² increases from 0.39 in 2014 to 0.63–0.72 in 2100 depending on the scenario. In the year 2100, this results in Arctic Ocean C_ant inventory estimates of 6.3±0.9 (SSP1-2.6), 8.2±1.2 (SSP2-4.5), 9.8±1.1 (SSP3-7.0), and 10.7±1.4 Pg C (SSP5-8.5) (Fig. 2b, d, f, and h). Thus, the emergent constraint is shown to slightly increase the CMIP6 multi-model mean projected C_ant inventory for each SSP (+2 % to +4 %) and substantially reduce associated uncertainties (−29 % to −31 %) resulting in greater separation of the SSPs (Fig. 1).

Figure 2 Emergent constraint on Arctic Ocean C_ant inventory in 2100. Maximum Arctic Ocean sea surface density (95th percentile) and Arctic Ocean C_ant inventory in 2100 for 12–14 ESMs for (a) SSP1-2.6 (blue), (c) SSP2-4.5 (green), (e) SSP3-7.0 (orange), and (g) SSP5-8.5 (red). The ordinary least-squares regressions (dashed lines) and the ±1σ prediction intervals (shaded area) are shown for each SSP. Probability density functions before and after applying the emergent constraint for (b) SSP1-2.6, (d) SSP2-4.5, (f) SSP3-7.0, and (h) SSP5-8.5.

3.2 Ocean acidification

3.2.1 Multi-model mean

Over the 21st century, ocean acidification leads to a reduction in Ω_arag and Ω_calc (Fig. 3). Until around 2040, the reduction of both saturation states averaged over the upper 1000 m is independent of the SSP. From 2040 onwards, the rate of the reduction depends on the pathway, and by 2100 basin-wide Ω_arag averaged over the upper 1000 m reaches 1.00±0.06 (SSP1-2.6), 0.82±0.05 (SSP2-4.5), 0.66±0.04 (SSP3-7.0), and 0.59±0.05 (SSP5-8.5), while Ω_calc reaches 1.56±0.09 (SSP1-2.6), 1.29±0.08 (SSP2-4.5), 1.06±0.06 (SSP3-7.0), and 0.92±0.05 (SSP5-8.5).

Figure 3 Basin-wide Arctic Ocean calcium carbonate saturation states over the 21st century. Time series of multi-model mean basin-wide saturation states of aragonite from 2000 to 2014 (black) and from 2015 to 2100 for SSP1-2.6 (blue), SSP2-4.5 (green), SSP3-7.0 (orange), and SSP5-8.5 (red) averaged (a) from 0 to 1000 m, (c) from 0 to 50 m, (e) from 50 to 500 m, and (g) from 500 to 1000 m with ±1 SD (n=12–14) shown as a shaded area. The differences when changing from historical simulations to the SSPs result from the different number of ensemble members that are available for each simulation. (b, d, f, h) The same time series but for calcite.

The upper 1000 m of the Arctic Ocean is projected to be on average undersaturated with respect to aragonite by 2100 under all of the SSPs, with Ω_arag in SSP1-2.6 just below 1 (0.995). Vertically, the entire water column will be undersaturated with respect to aragonite in 2100 for all pathways but SSP1-2.6 (Fig. 4). Under this low-emissions pathway, Arctic Ocean waters above 500 m are projected to remain supersaturated with respect to aragonite, while those below 500 m are projected to be undersaturated. For the more stable calcium carbonate polymorph calcite, basin-wide Arctic Ocean undersaturation is only projected under SSP5-8.5. Under this high-emissions pathway, basin-wide undersaturation is projected for the water masses between 50 and 800 m.

Figure 4 Vertical profile of basin-wide Arctic Ocean calcium carbonate saturation states, pH, and pCO₂ in 2100. Multi-model mean vertical profiles of (a) aragonite and (b) calcite saturation states, (c) pH, and (d) pCO₂ in 2100 for SSP1-2.6 (blue), SSP2-4.5 (green), SSP3-7.0 (orange), and SSP5-8.5 (red) averaged (a) with ±1 SD (n=12–14) shown as a shaded area. Observed profiles in 2002 from GLODAPv2 are shown in black (Lauvset et al., 2016).

Alongside declines in calcium carbonate saturation states, pH and pCO₂ are also projected to change over the 21st century. Basin-wide averaged pH in the upper 1000 m decreases from 8.06 in 2002 to 7.88±0.02 (SSP1-2.6), 7.79±0.03 (SSP2-4.5), 7.68±0.03 (SSP3-7.0), and 7.61±0.03 (SSP5-8.5), while basin-wide averaged pCO₂ in the upper 1000 m increases to 739±38 (SSP1-2.6), 918±75 (SSP2-4.5), 1209±92 (SSP3-7.0), and 1428±124 µatm (SSP5-8.5). The lowest pH values and highest pCO₂ values are projected to occur between 100 and 1000 m from the surface (Fig. 4).

3.2.2 Emergent constraints on acidification

The emergent relationship between observed maximum sea surface density and end-of-century Arctic Ocean acidification that was previously identified in CMIP5 models (Terhaar et al., 2020a) does not exist anymore in the new CMIP6 model generation (Fig. 5). Despite the emergent constraint approach still functioning for projections of the Arctic Ocean C_ant inventory in CMIP6, this no longer translates into an emergent constraint for Arctic Ocean Ω_arag and Ω_calc. However, if all changes in ocean biogeochemistry with the exception of alkalinity, which is kept constant at observed levels from 2002 (Lauvset et al., 2016), are considered when estimating changes in Ω_arag and Ω_calc, emergent relationships are still found between maximum sea surface density and Ω_arag across all SSPs (r²=0.54–0.73). When decreases in alkalinity (A_T) are also taken into account, Ω_arag and Ω_calc decrease further, and the emergent relationship disappears (r²=0.01–0.12).

Figure 5 The CMIP6 ensemble relationship between end-of-century Ω_arag and the C_ant inventory in the Arctic Ocean. Basin-wide Ω_arag averaged over the upper 1000 m with (orange) and without (blue) taking into account concurrent changes in alkalinity (A_T) against the Arctic Ocean C_ant inventory in the year 2100 for (a) SSP1-2.6, (a) SSP2-4.5, (c) SSP3-7.0, and (d) SSP5-8.5.

3.3 Comparison between CMIP5 and CMIP6

Compared to ESMs from CMIP5, the new generation of ESMs (CMIP6) has improved in simulating the maximum Arctic Ocean sea surface density. Specifically, negative density biases have been reduced, and the inter-model range in maximum sea surface density has substantially decreased from 3.6 in CMIP5 (Terhaar et al., 2020a) to 0.9 kg m⁻³ in CMIP6 (Fig. 2). As a result, the inter-model range of the C_ant inventory in the CMIP6 model ensemble is also reduced (Fig. 2). Moreover, without the negative maximum-sea-surface-density bias, the simulated multi-model mean Arctic Ocean C_ant inventory in 2005 is 61 % higher than the inventory that was simulated by the previous model generation (1.3±0.7 Pg C) (Terhaar et al., 2020a).

At the end of the 21st century, the unconstrained simulated Arctic Ocean C_ant inventory under SSP5-8.5 is 37 % larger and the uncertainty is 19 % smaller than the unconstrained simulated C_ant inventory in 2100 under RCP8.5 (7.5±2.9 Pg C) (Terhaar et al., 2020a). After applying the constraint, the Arctic Ocean C_ant inventory in 2100 under SSP5-8.5 (10.7±1.4 Pg C) is 19 % larger than the constrained Arctic Ocean C_ant inventory under RCP8.5 (9.0±1.6 Pg C). This difference is of the same order of magnitude as the difference in prescribed atmospheric CO₂ concentration over the 21st century, which is higher in SSP5-8.5 (CMIP6) than RCP8.5 (CMIP5) (Meinshausen et al., 2011, 2019) and therefore results in greater surface ocean acidification for approximately the same radiative forcing (Kwiatkowski et al., 2020).

To compare the emergent constraint across scenarios with different atmospheric CO₂ concentrations, the simulated Arctic Ocean C_ant inventory in 2100 for each scenario was rescaled to SSP5-8.5 using the mean atmospheric C_ant concentration from 1850 to 2100 as a linear scaling factor (Fig. 6). The relationship remains robust (r²=0.63–0.74) for all five analysed scenarios. The slope of the emergent relationship is however substantially steeper in CMIP6 (9.4–12.6 Pg C kg⁻¹ m³) than in CMIP5 (3.3 Pg C kg⁻¹ m³). However, the slope in CMIP5 increases to 8.9 Pg C kg⁻¹ m³ if the three CMIP5 models with particularly low maximum sea surface densities (<27.5 kg m⁻³) are excluded (dotted line in Fig. 6). The resulting constrained estimate for the rescaled Arctic Ocean C_ant inventory decreases from the low-emission scenario to the high-emission scenario from 12.3 to 10.7 Pg C. When comparing the two high-emission scenarios, the rescaled Arctic Ocean C_ant inventories are 10.7 (SSP5-8.5) and 10.4 Pg C (RCP8.5). The latter remains unchanged if the three CMIP5 models with particularly low maximum sea surface densities (<27.5 kg m⁻³) are excluded.

Figure 6 Arctic Ocean C_ant inventory in 2100, rescaled to SSP5-8.5 atmospheric C_ant concentrations (see Methods), against maximum Arctic Ocean sea surface density (95th percentile) for SSP1-2.6 (blue diamonds; r²=0.70), SSP2-4.5 (green squares; r²=0.63), SSP3-7.0 (orange triangles; r²=0.70), SSP5-8.5 (red circles; r²=0.72), and RCP8.5 (purple stars; r²=0.74) (a) over the entire density range and (b) over the density range from 27.5 to 28.3 kg m⁻³. The simulated Arctic Ocean C_ant inventory in 2100 for each scenario (Fig. 2) was rescaled to SSP5-8.5 using the mean atmospheric C_ant concentration from 1850 to 2100 as a scaling factor. Linear fits use all available models for each scenario (coloured dashed lines). For RCP8.5 an additional fit is shown excluding the three models with density anomalies below 27.5 kg m⁻³ (b) (dotted purple line; r²=0.79).

4 Discussion

4.1 Arctic Ocean C_ant inventory

Across the CMIP6 model ensemble, the Arctic Ocean C_ant storage over the 21st century is highly related to maximum sea surface densities (Fig. 2), which predominately occur in the Barents Sea (Midttun, 1985; Smedsrud et al., 2013; Terhaar et al., 2020a). The inter-model range in maximum sea surface density in the Barents Sea is mainly explained by differences in sea surface salinities (r²=0.93), which are influenced by brine rejection (Midttun, 1985; Årthun et al., 2011) and the strength of inflowing, saltier Atlantic waters through the Barents Sea Opening (Fig. 7). Compared to CMIP5 models, the reduced negative bias of simulated maximum sea surface densities (Fig. 6) indicates model improvement in simulating the circulation of Atlantic and Arctic surface waters. Despite the reduced inter-model range, the robust relationship between maximum sea surface density and C_ant across model generations (Fig. 6) supports evidence that inflowing Atlantic waters through the Barents Sea Opening and their transformation into deep and intermediate waters via brine rejection are the dominant process governing Arctic Ocean C_ant increases (Midttun, 1985; Rudels et al., 1994, 2000; Jeansson et al., 2011; Smedsrud et al., 2013). The application of observational constraints to this emergent multi-model relationship in order to constrain the projected Arctic Ocean C_ant inventory (Fig. 2) and focus efforts on model development therefore remains promising.

Figure 7 Annual mean sea surface salinity in the Barents Sea and its neighbouring waters from (a) the World Ocean Atlas 2018 and from the CMIP6 model with the (b) highest (UKESM1-0-LL) and (c) lowest (CESM2-WACCM) maximum sea surface density and salinity. The 0 ^∘C isotherm indicates approximately the extent of Atlantic waters (black solid line).

However, the slope of the linear relationship between maximum sea surface density and the Arctic Ocean C_ant inventory over the 21st century in the CMIP5 model ensemble (3.3±0.6 Pg C kg⁻¹ m³, scaled to SSP5-8.5 atmospheric C_ant concentrations) is 3–4 times less than that in the CMIP6 model ensemble (12.6±2.6 for SSP1-2.6, 12.1±2.8 for SSP2-4.5, 9.4±1.9 for SSP3-7.0, and 11.3±2.0 Pg C kg⁻¹ m³ for SSP5-8.5) (Fig. 6). The reduced slope in the CMIP5 ensemble is mainly caused by three models with maximum surface density anomalies well below 27.5 kg m⁻³. When these three models are excluded, the remaining CMIP5 models follow a slope of 8.9±2.2 Pg C kg⁻¹ m³, in broad agreement with the CMIP6 model ensemble (Fig. 6). This suggests that the linear emergent relationship does not hold below a certain value of maximum sea surface density below which the impact on deep-water formation and subsequent C_ant storage in the Arctic Ocean is limited. However, as the two linear relationships happen to cross the observed maximum sea surface density at nearly the same location (Fig. 6), the constrained C_ant inventory for the CMIP5 model ensemble remains almost entirely unchanged when the three low-density models are excluded.

Even without the low-density bias in the Barents Sea, the constrained Arctic Ocean C_ant inventory in 2005 in CMIP6 remains 36 % below the data-based estimate (Terhaar et al., 2020b). This underestimation is partly due to the different definition of C_ant in data-based estimates and ESMs. While the historical simulations in CMIP5 and CMIP6 typically start in 1850, data-based estimates account for all C_ant since 1765. This leads to an underestimation of the global ocean C_ant inventory by ESMs of around 30 % (Bronselaer et al., 2017) and of around 20 % in the Arctic Ocean (Terhaar et al., 2019b). Even if we increased the constrained Arctic Ocean C_ant inventory in 2005 by 20 %, an underestimation of around 16 % would remain compared to the data-based estimate. This underestimation of the data-based estimate suggests that all ESMs are missing additional pathways of C_ant entry into the Arctic Ocean, other than the principal pathway via the Barents Sea. Indeed, small-scale density flows along continental slopes can be observed in different regions of the Arctic Ocean (Rudels et al., 1994; Jones et al., 1995) but cannot be simulated by the coarse resolution of most ESMs. Thus, the constrained estimates of the Arctic Ocean C_ant inventory presented here are likely still a lower boundary.

Recent observation of dilution of A_T and C_T in surface waters of the Amerasian Basin caused by freshening have led to the hypothesis that continuous freshening might turn the Arctic Ocean from a sink of C_ant into a source over the 21st century (Woosley and Millero, 2020). However, observations in the Eurasian basins, which receive more saline Atlantic water and less freshwater input, still show increases of C_ant concentrations (Ulfsbo et al., 2018) over a depth of 1500 m in the last 20 years. The CMIP5 and CMIP6 model ensembles both simulate continuous accumulation of C_ant in the Arctic Ocean under all SSPs (Fig. 1), suggesting that the subduction of C_ant-rich Atlantic waters in the Barents Sea remains larger than any loss of C_ant in surface waters over the 21st century. Nevertheless, the reduction of the storage rate of C_ant under SSP5-8.5 (Fig. 1) in combination with constantly increasing atmospheric CO₂ concentrations (Riahi et al., 2017) indicates that dilution may reduce the capacity of the Arctic Ocean to store further C_ant as suggested by Woosley and Millero (2020).

4.2 Arctic Ocean acidification

4.2.1 Ω_arag and Ω_calc, pH, and pCO₂ in 2100

Even under the most optimistic scenario assessed (SSP1-2.6), the Arctic Ocean will become on average undersaturated with respect to aragonite, with possible consequences for calcifying organisms (Comeau et al., 2010) and the food chain (Armstrong et al., 2005; Karnovsky et al., 2008). In this scenario only water masses above 500 m remain supersaturated this century. These findings are in good agreement with idealized estimates of Ω_arag and Ω_calc that project the future C_ant inventory based on the transient time distribution method and observed CFC-12 concentrations (Anderson et al., 2010; Terhaar et al., 2020b).

In addition to widespread Arctic Ocean undersaturation with respect to aragonite, water masses between 50 and 800 m are even projected to become undersaturated with respect to calcite by the end of the century under SSP5-8.5. Calcite undersaturation in Arctic Ocean subsurface waters will likely further enhance the pressure on the Arctic Ocean ecosystem as calcite-forming organisms, such as foraminifera (Davis et al., 2017) and coccolithophores (Kottmeier et al., 2016), experience potential impacts on growth and survival.

Furthermore, the projected increases in pCO₂ (Fig. 4c and d) alongside projected increases in its seasonal amplitude (McNeil and Sasse, 2016; Kwiatkowski and Orr, 2018) are likely to lead to hypercapnic conditions that might directly affect the growth and survival of Arctic fish (Frommel et al., 2012; Schmidt et al., 2017; Kunz et al., 2018) under high-emissions scenarios.

4.2.2 Driving processes of acidification

The CMIP6 model ensemble shows astonishingly good agreement with respect to projections of Arctic Ocean acidification over the 21st century. Compared to the CMIP5 model ensemble, the uncertainties in projected Ω_arag averaged over the upper 1000 m in 2100 have been reduced from 0.13 to 0.04–0.06, and those for Ω_calc from 0.21 to 0.05–0.09.

The main reason for the reduced uncertainty appears to be that the ESMs that have a lower maximum sea surface density and that thus take up less C_ant in the Arctic Ocean over the 21st century (Fig. 2) are the same models that simulate a stronger reduction in A_T (Fig. 8). To explain this negative correlation between C_ant uptake and reductions in A_T, we propose the following mechanism: over the 21st century, melting of sea ice, land ice, and increased river runoff are projected to freshen the Arctic Ocean (Koenigk et al., 2013; Nummelin et al., 2016; Shu et al., 2018) and to reduce A_T (Fig. 8) (Woosley and Millero, 2020). In the CMIP6 ensemble, this freshening and reduction of A_T tends to be stronger in models with lower sea surface salinities and densities; i.e. the models that simulate less inflow of saline Atlantic waters, such as CESM2-WACCM (Fig. 7), have lower Barents Sea surface salinity and density, weak deep-water formation, and therefore less C_ant storage (Fig. 2). In contrast, models with a larger inflow of saline Atlantic water have stronger deep-water formation, greater C_ant storage, less freshening, and a smaller reduction in A_T. This compensation then results in a similar simulated reduction of Ω_arag and Ω_calc over the 21st century across all models in the CMIP6 ensemble, through either a reduction in A_T or an increase in C_T. Other possible drivers, such as changes in temperature and salinity, are likely of minor importance given the relatively good correlation between the projected C_ant inventory and basin-wide reduction of Ω_arag and Ω_calc over the 21st century when A_T is not taken into account (r²=0.54–0.73).

Figure 8 CMIP6 basin-wide changes in alkalinity (A_T) averaged over the upper 1000 m from 2002 to 2100 caused by freshwater dilution (blue) and changes in biogeochemistry (green) against the Arctic Ocean C_ant inventory for (a) SSP1-2.6, (b) SSP2-4.5, (c) SSP3-7.0, and (d) SSP5-8.5. A_T anomalies were decomposed by calculating the temporal evolution of salinity-corrected alkalinity following Lovenduski et al. (2007).

In the CMIP5 model ensemble, this compensation effect did not exist, and projected saturation states and pH were mainly driven by the projected C_ant storage in each model (Terhaar et al., 2020a), with changes in A_T, temperature, and salinity being of minor importance. The relatively small changes in A_T in CMIP5 compared to CMIP6 might have been caused by an underestimation of Arctic Ocean freshening over the 21st century given that the CMIP5 models also underestimated historical Arctic Ocean freshening by around 50 % in the Arctic Ocean (Shu et al., 2018). In the CMIP6 model ensemble this negative bias with respect to freshwater fluxes and subsequent dilution of A_T appears to have been reduced. This is probably a consequence of better representation of riverine fluxes as well as increased model resolution and hence improved circulation (Séférian et al., 2020).

To estimate the likely impacts of C_T increases and A_T decreases to Arctic Ocean acidification over the 21st century, one can assess their relative contributions to acidification in constrained estimates of end-of-century C_ant (dotted line in Fig. 5). In SSP5-8.5 for example, increasing C_T along with changes in temperature and salinity reduce Ω_arag from 1.38 to 0.75, while decreasing A_T further reduces this to 0.59. Thus, the future storage of C_ant in the Arctic Ocean still remains the main driver of Arctic Ocean acidification, with acidification further enhanced by decreases in A_T.

5 Conclusions

Earth system models tend to have improved their performance in the Arctic Ocean from CMIP5 to CMIP6. The negative bias with respect to present-day maximum sea surface density is substantially reduced in CMIP6, and the historical C_ant inventory is closer to the data-based estimate (Terhaar et al., 2020b). Furthermore, the inter-model range of projected Ω_arag and Ω_calc and the C_ant inventory have been reduced. Uncertainties in the future C_ant inventory under four shared socioeconomic pathways can be further reduced by applying the previously identified emergent constraint between present-day maximum Arctic Ocean sea surface density and the future C_ant inventory (Terhaar et al., 2020a).

Hall et al. (2019) define “confirmed” emergent constraints as fulfilling the following criteria: (1) a plausible mechanism, (2) verification of that mechanism, and (3) out-of-sample testing. The mechanism underpinning the relationship between maximum sea surface densities and the projected C_ant inventory in the Arctic Ocean is consistent with observations and has been verified in an ocean-only model at different resolutions (Terhaar et al., 2020a). The testing of this emergent constraint in the CMIP6 ensemble can be considered out of sampling, although the extent of model independence between generations remains questionable (Knutti et al., 2013; Sanderson et al., 2015). Despite the overall emergent constraint on projected Arctic Ocean C_ant being similar between CMIP6 and CMIP5, the change in the multi-model emergent relationship indicates that further validation is required. Specifically, it appears that the linear relationship between maximum sea surface densities and the projected C_ant inventory breaks down in models with extremely low-biased sea surface densities, highlighting a potential limitation to this emergent constraint that requires further assessment.

While the mechanisms underlying emergent constraints on the future Arctic Ocean C_ant inventory appear to be consistent between CMIP5 and CMIP6, this is not the case for the previously identified emergent constraint on ocean acidification variables (Ω_arag and Ω_calc, pH, and pCO₂). In the CMIP6 model ensemble, projected model uncertainties in Ω_arag and Ω_calc, pH, and pCO₂ are dramatically reduced compared to the uncertainties in CMIP5 and are driven not only by the increase in C_ant but also by a reduction in A_T due to pan-Arctic freshening. A weak inflow of saline Atlantic waters in ESMs results not only in lighter surface waters and less C_ant storage over the 21st century but also in a stronger reduction of alkalinity caused by a stronger freshening. A strong Atlantic inflow, however, appears to limit freshening and reductions in alkalinity but results in greater C_ant storage. Although both C_ant storage and alkalinity reductions contribute considerably to Arctic Ocean acidification, our results suggest that C_ant remains the dominant process. In the CMIP5 model ensemble, the influence of freshening on emergent constraints on acidification variables might have been of limited importance because (a) freshening was largely underestimated (Shu et al., 2018) and (b) differences in simulating the inflow of Atlantic waters through the Barents Sea Opening were considerably larger across the model ensemble (Terhaar et al., 2020a). Thus, the absence of an emergent constraint on projected Arctic acidification in CMIP6 could be viewed as a consequence of successful model development. Indeed, there is clearly no need for such constraints when models are in broad agreement.

Independent of the driving mechanism, the projected Arctic Ocean acidification over the 21st century has possibly grave consequences for the wider Arctic Ocean ecosystem (Gattuso and Hansson, 2011; Riebesell et al., 2013; AMAP, 2018). Under all socioeconomic pathways, the upper 1000 m of the Arctic Ocean will be on average undersaturated with respect to aragonite. Thus, keystone species like the sea butterfly (aragonitic pteropod Limacina helicina) may lose most to all of their suitable habitat in the Arctic Ocean (Comeau et al., 2010). However, under the low-emissions pathway SSP1-2.6 aragonite undersaturation is projected to be avoided in the upper 500 m of the water column, providing a potential refuge. Under the high-emissions pathway SSP5-8.5, water masses between 50 and 800 m are additionally projected to become undersaturated with respect to calcite. While isolated regions of the surface Arctic Ocean that are strongly influenced by riverine inputs are already seasonally undersaturated with respect to calcite (Bates and Mathis, 2009), subsurface Arctic waters would be the first to exhibit annually averaged and basin-wide calcite undersaturation over a depth of several hundred metres. Vertically migrating organisms that form calcite shells and skeletons, such as coccolithophores (Kottmeier et al., 2016) and foraminifera (Davis et al., 2017), may lose their natural refugia to which they migrate during the day (Berge et al., 2015). These hostile conditions for calcifying organisms will likely cause their decline, with unknown consequences for the wider Arctic food web (Armstrong et al., 2005; Karnovsky et al., 2008) and its iconic species (Jay et al., 2011; AMAP, 2018).