Model implementation

The boxes used in our model (Fig. 2) represent both the Amerasian Basin (4,250 m maximum depth) and the Eurasian Basin (4,750 m maximum depth), separated by the Lomonosov Ridge. Each basin is divided into surface (0–200 m), intermediate (200–700 m) and deep (>700 m) water boxes, and there exist flows between these boxes and from external sources, that is, rivers, as well as the Atlantic and Pacific Oceans. Boxes are also included to account for sediment accumulation and benthic CaCO 3 dissolution. The justification for the use of a box model is provided in the Methods section.

Figure 2: The box model for the carbonate system of the Arctic Ocean. The model divides the Arctic Ocean into a collection of six boxes, that is, surface Amerasian Basin (S A ), intermediate Amerasian Basin (I A ), deep Amerasian Basin (D A ), surface Eurasian Basin (S E ), intermediate Eurasian Basin (I E ), deep Eurasian Basin (D E ), as well as an atmospheric box and two sediment reservoirs (S AB and S EB ) connected to the deep water boxes. Water flows between boxes are designated by capital U’s, the capital E’s indicate gas exchanges, the B’s stand for internal fluxes of CaCO 3 and the F’s are the fluxes of CaCO 3 to the sediments. Subscripts refer to the particular ocean (A and E), the different depth boxes (S, I and D) and the nature of the flux, for example, T for thermohaline, M for mixing, D for dissolution and so on. In addition, DBC stands for deep boundary current and PDC abbreviates previously deposited carbonate. Full size image

Our model equations account for and predict the changes in total dissolved CO 2 (∑CO 2 ) and carbonate alkalinity (CA) in each of these boxes, as forced by the increasing CO 2 in the atmosphere and in water feeding into the Arctic. The model subsequently calculates the pH (National Bureau of Standards) and the aragonite saturation state—see the Methods section with regard to the pH scale. Time-varying ∑CO 2 and CA of the Pacific surface and high-latitude Atlantic waters that enter the Arctic are obtained from the output of a previously published global carbon-system model24,25. The atmospheric pCO 2 is also calculated by that same model, in a procedure similar to that used by Spall21.

These calculations start with initial water flows, concentrations and parameter values (Methods) believed to have been in place during pre-industrial times, that is, pre-1850 AD (all dates are AD and that suffix is dropped hereafter), as calculated by correcting and averaging compiled ∑CO 2 and total alkalinity (TA) data15,16 (Supplementary Figs 1 and 2, and Supplementary Table 1). CA is calculated from TA through standard methods (that is, from pH, total boron concentrations and appropriate dissociation constants). Flow values were obtained as explained in the Methods section, and other assigned parameters can be found in Supplementary Table 1.

The increase in pCO 2 , which acidifies the oceans, is driven with an extended version of the IS92a emission scenario. That extended IS92a scenario was used in Boudreau et al.24 and is described and displayed in Supplementary Fig. 3. We emphasize that this is not the original IS92a scenario, but a version modified to be consistent with emissions data and atmospheric CO 2 levels to 2010 and extended to predict emissions over the next millennia. This scenario may overestimate future emissions, but precaution demands that we focus on the worst possible case.

Carbon dioxide absorption from the atmosphere into the surface Arctic Ocean (dark blue E S,A and E S,E in Fig. 2) is thought to be hindered by persistent ice cover3. Thus, as global warming reduces the extent and duration of ice cover, CO 2 adsorption may increase, although that has been debated17,18, as noted above; nevertheless, we include ice melting and corresponding increased CO 2 uptake in our model. The evolution of the ice cover with warming has been modelled with both the fast and slow melting scenarios by Yamamoto et al.3 (Supplementary Fig. 4). Both scenarios were tested and our results were identical in both cases, hence, we only report the results with fast melting.

Model output and observed data

Whereas the model results are in the form of ∑CO 2 and CA, we are primarily interested in the evolving saturation state of the Arctic waters with respect to aragonite, Ω a

where [Ca2+] is the calcium concentration, [CO 3 2−] is the carbonate ion concentration and K* sp is the stoichiometric solubility product under in situ conditions. A similar equation applies to calculation of the saturation state with respect to calcite. [Ca2+] is calculated from the known salinities of the Arctic Ocean basins. Note that, irrespective of the water chemistry, the aragonite (and calcite) saturation state changes with depth (pressure) in response to the increasing solubility of carbonate minerals25,26, that is, K* sp increases with depth (pressure).

The aragonite saturation horizon (Z sat ) is the depth above which the waters are supersaturated with respect to aragonite (Ω a >1) and below which waters are undersaturated (Ω a <1) and in which aragonite will dissolve. To predict the position of that horizon, we coupled our box model to an explicit formula24,27 for Z sat ,

where Z1 sat is a characteristic depth calculated from the solubility equations25,26 for aragonite and K1 sp is the value of K* sp at 1 atm. The short derivation of this latter equation is repeated for completeness in the Methods; note that equation (2) is not dependent on the form of the model being used, that is, a box versus a 3D model.

Central to our interests are the relative roles of atmospheric CO 2 forcing (E) and the input of waters from the Atlantic and Pacific Oceans (U T,IA , U T,SP in Fig. 2) in changing the pH and Ω a of the Arctic Ocean. To facilitate this analysis, we created a hypothetical reference state wherein the pCO 2 in the Arctic atmosphere follows our prescribed CO 2 emissions scenario (Fig. 3a and Supplementary Fig. 3), but ∑CO 2 and CA of the inflowing Atlantic and Pacific waters do not change with time, labelled ‘constant source’ in our figures; these particular inflows are set to year 2010 values. In contrast, a more realistic model results by allowing Atlantic and Pacific water chemistries to change as dictated by the evolving atmospheric CO 2 , and those model results are presented with the label ‘variable source’.

Figure 3: Calculated pH and surface CaCO 3 saturation of the Arctic basins. (a) The pCO 2 created by the extended IS92a carbon emission scenario24 and the ice-free fraction of the surface Arctic Ocean3 as used in these calculated results. (b,c) pH in surface box and deep box (average) in Amerasian Basin and Eurasian Basin, respectively. (d,e) Carbonate saturation state for both aragonite and calcite in surface box in Amerasian Basin and Eurasian Basin, respectively. Full size image

The model predicted evolution of the pH (NBS) in the surface and deep waters is illustrated in Fig. 3, assuming our extended IS92a CO 2 emissions scenario (Fig. 3a) and constant productivity for both organic matter and CaCO 3 (Methods). The surface results are obtained without freshening from melting ice and increased run-off and must be accepted with this caveat.

pH is an inexact measure of the carbonate chemistry of marine waters; consequently, we have also calculated the aragonite (orange lines) and calcite (green lines) saturation states of the surface water, as shown in Fig. 3d,e. Surface waters will become undersaturated with respect to aragonite (Ω a <1), with a minimum between 0.75 and 0.675 attained close to the year 2200 in both basins.

To test our model, we compare our output with available saturation data. Figure 4a is a contour plot of the aragonite saturation state for the Amerasian Basin calculated from the data compilations of Jutterstrom et al.28 and Key et al.15,16, as illustrated in Supplementary Figs 1 and 2. These data indicate a rise of about 500 m over the 1995–2009 period. Figure 4b reproduces the saturation isopleths provided by the 3D model used by Frölicher and Joos14, which are essentially flat over the 1995–2009 data time series. Finally, Fig. 4c displays our model predictions of Ω a isopleths for that same period and reveals that the isopleths are sloped by an amount similar to the data in Fig. 4a, with an upward displacement of ∼300 m, based on the Ω a =1.2 isopleth. Our model predicts slightly more acidic deep water in year 1995 than the data, but the observed and modelled isopleths below 2,000 m are roughly at the same depth in 2009.

Figure 4: Contour plots of aragonite saturation in the Arctic Ocean. Lines are isopleths of constant Ω a . (a) Calculated from the data in Jutterstrom et al.28 and Key et al.15,16. (b) Calculated with a 3D ocean-circulation-biogeochemistry model by Frölicher and Joos14. (c) Calculated with equation (2) and the carbonate chemistry from our model (Fig. 2). (Note the 1-year time shift in the plot for data compared with the plots of the model results; this is an artefact of the contour.m function in Matlab.) Full size image

Finally, Fig. 5 displays our prediction of the long-term evolution of the aragonite saturation horizon (red lines) and that made by Frölicher and Joos14 (black lines). These plots are explained and analysed below.