# Arctic Ocean annual high in pCO2 could shift from winter to summer

James Orr, Lester Kwiatkowski, Hans-Otto Pörtner

► To cite this version:

James Orr, Lester Kwiatkowski, Hans-Otto Pörtner. Arctic Ocean annual high in pCO2 could shift from winter to summer. Nature, 2022, 610 (7930), pp.94-100. 10.1038/s41586-022-05205-y. hal-03833044

HAL Id: hal-03833044

https://hal.science/hal-03833044

Submitted on 22 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Arctic Ocean annual high in $p_{CO_2}$ could shift from winter to summer

https://doi.org/10.1038/s41586-022-05205-y

Received: 2 December 2021

Accepted: 8 August 2022

Published online: 5 October 2022

Open access

<figure><img src="image_2.png" /></figure>

Check for updates

James C. Orr<fnref n="1" />,²,³, Lester Kwiatkowski² & Hans-Otto Pörtner³

Long-term stress on marine organisms from ocean acidification will differ between seasons. As atmospheric carbon dioxide ($CO_2$) increases, so do seasonal variations of ocean $CO_2$ partial pressure ($p_{CO_2}$), causing summer and winter long-term trends to diverge<fnref n="1" />⁻⁵. Trends may be further influenced by an unexplored factor—changes in the seasonal timing of $p_{CO_2}$. In Arctic Ocean surface waters, the observed timing is typified by a winter high and summer low⁶ because biological effects dominate thermal effects. Here we show that 27 Earth system models simulate similar timing under historical forcing but generally project that the summer low, relative to the annual mean, eventually becomes a high across much of the Arctic Ocean under mid-to-high-level $CO_2$ emissions scenarios. Often the greater increase in summer $p_{CO_2}$, although gradual, abruptly inverses the chronological order of the annual high and low, a phenomenon not previously seen in climate-related variables. The main cause is the large summer sea surface warming⁷ from earlier retreat of seasonal sea ice⁸. Warming and changes in other drivers enhance this century's increase in extreme summer $p_{CO_2}$ by 29 ± 9 per cent compared with no change in driver seasonalities. Thus the timing change worsens summer ocean acidification, which in turn may lower the tolerance of endemic marine organisms to increasing summer temperatures.

More than a decade ago, it was recognized that the ongoing increase in atmospheric carbon dioxide ($CO_2$) would affect not only the annual mean state of ocean $CO_2$ system variables but also their seasonal variations. Theory suggests that seasonal variations in ocean pH and $CO_2$ partial pressure ($p_{CO_2}$) should increase⁹ as more $CO_2$ invades the ocean and its buffer capacity is reduced, consistent with earlier model and mesocosm studies¹⁰,¹¹. Subsequent observations and model projections confirm the proposed increase in the seasonal amplitude of $p_{CO_2}$ across the ocean, although the Arctic Ocean was either excluded from the analysis<fnref n="1" />,⁴ or included while presuming that effects from physical climate change were negligible². For pH, models project that generally the amplitude of seasonal variations will actually decline because its changes are relative, not absolute changes in hydrogen ion concentration [H⁺] owing to the logarithmic transformation; conversely, seasonal variations in [H⁺] increase nearly linearly with those for $p_{CO_2}$ (refs. ³,⁵).

Here we have assessed monthly variations and their potential future changes in timing for surface ocean $p_{CO_2}$ and related ocean $CO_2$ system variables in the Arctic Ocean using a suite of 27 Earth system models (ESMs) that participated in the two most recent phases of the Coupled Model Intercomparison Project (CMIP5 and CMIP6; Supplementary Table 1). Driving mechanisms were assessed with Taylor-series expansions and idealized experiments from CMIP5 to distinguish the direct chemical consequences attributable to the increase in atmospheric $CO_2$ from the indirect consequences of physical climate change (Methods).

Here these tools are used to focus on the timing of the annual cycle and its future change in regards to chemical and biological impacts.

## Timing of annual cycle

For the modern ocean, the timing of the annual cycle of $p_{CO_2}$ and related variables is generally well understood. In the subtropics, where the annual cycle of $p_{CO_2}$ is dominated by temperature-driven variations, observations indicate that $p_{CO_2}$ levels are consistently higher in summer and lower in winter, whereas subpolar regions have the opposite pattern because non-thermal effects dominate¹²,¹³. Similar seasonal patterns and driving mechanisms are simulated by the CMIP5 models under modern forcing³,⁴ except in the Southern Ocean where models struggle¹⁴,¹⁵. However, modelled seasonal variations in $p_{CO_2}$ have not been evaluated in the Arctic Ocean owing to the sparsity of seasonal observations. To fill this observational gap, we exploited a recent neural-network-derived, high-resolution, observation-based product⁶ that includes the non-coastal Arctic Ocean. It indicates that surface-water $p_{CO_2}$ is lower in summer than in winter when averaged across the Arctic, suggesting that non-thermal effects dominate as in the subarctic¹²,¹³ (Supplementary Fig. 1). The models generally show the same tendency for present-day monthly variations of the basin-wide mean, a major component of overall seasonal variability (Supplementary Fig. 2) that correlates strongly with the observation-based variations (Supplementary Fig. 3). This consistency improves confidence in the future projections.

<footnote n="1">¹LSCE-IPSL, Laboratoire des Sciences du Climat et de l'Environnement, CEA-CNRS-UVSQ, CEA Saclay, Gif-sur-Yvette, France. ²LOCEAN-IPSL, Laboratoire d'Océanographie et du Climat: Expérimentations et Approches Numériques, Sorbonne Université-CNRS-IRD-MNHN, Paris, France. ³Integrative Ecophysiology, Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Bremerhaven, Germany. ⁴e-mail: james.orr@lsce.ipsl.fr</footnote>
<figure><img src="image_3.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 1 | Projected seasonal timing of p'_{CO₂} changes little except in the Arctic Ocean.** *a–c*, Curves represent the CMIP5 model mean for 1996–2005 (historical) and 2091–2100 (RCP8.5) detrended decadal climatologies averaged over the Arctic Ocean (*a*), the northern subtropics (*b*) and the southern subtropics (*c*).</figcaption></figure>

Under high-end emissions, the models project that by 2100, the thermally driven, summer maximum in $p'_{CO_2}$ in the low latitudes occurs about a month earlier, as does the biologically driven, summer minimum in the high latitudes (Fig. 1). However, in the Arctic, today's broad summer minimum in the monthly anomaly of $p'_{CO_2}$, that is, relative to the annual mean and denoted as $p'_{CO_2}$, tends to narrow in the future projections of all models, often splitting into a spring–summer minimum and a summer maximum (Fig. 2 and Extended Data Fig. 1). That projected earlier minimum in $p'_{CO_2}$ usually occurs along with an earlier peak in net primary production (NPP), whereas the summer $p'_{CO_2}$ maximum occurs along with a sharp increase in summer sea surface temperature (7). Both appear to be tied to earlier retreat of seasonal sea-ice cover. These projected tendencies are consistent with observations in some Arctic regions that are already experiencing earlier reductions

in sea-ice cover, where phytoplankton blooms occur earlier in the year¹⁶ and there is greater surface warming⁸. Observed summer $p'_{CO_2}$ is higher in low-ice years¹⁷, has grown over recent decades as sea ice retreats¹⁸ and is particularly high under exceptional warming¹⁹⁻²¹.

The projected sign reversal in the summer extreme of $p'_{CO_2}$ occurs in most models under high-end emissions, in many models under mid-range emissions and in some models under low-end emissions (Extended Data Fig. 2).

The evolution in timing of the annual low and high in Arctic Ocean $p'_{CO_2}$ differs from other variables such as NPP and T, which exhibit little or no timing change (Fig. 3a). Throughout the industrial era until present, the annual high in $p'_{CO_2}$ occurs in April for the CMIP5 mean averaged over the Arctic domain. But as atmospheric $CO_2$ continues to increase, the timing of that annual high first retreats to March and then advances rapidly, finishing in September at the end of the century (936 ppm under Representative Concentration Pathway 8.5 (RCP8.5)). Simultaneously, the timing of the annual low retreats from August to June, crossing over the advancing annual high at 827 ppm. The same abrupt transition is found for the CMIP6 model mean but crossover occurs at 571 ppm, and the advance of the timing of the annual high to September reverses itself at about 700 ppm and then recedes to August (Fig. 3b). Although the crossover $CO_2$ level differs by 256 ppm between CMIP5 and CMIP6 means, and by more among models as shown below, it is remarkably consistent across scenarios for a given model. Despite the abrupt transition in the annual maximum $p'_{CO_2}$ from winter to summer, the increase in summer $p'_{CO_2}$ is a gradual evolution (Fig. 3c,d), the range of which seems to contain the different model behaviours seen in 2091–2100 (Fig. 2 and Extended Data Fig. 1).

Crossover for $p'_{CO_2}$ occurs in most models under mid- and high-level emissions (Extended Data Fig. 3) but never for related drivers (Extended Data Fig. 4). For CMIP5, five of nine models reach $p'_{CO_2}$ crossover, in June or July at atmospheric $CO_2$ levels that vary from 400 ppm to 922 ppm. Two CMIP5 models do not reach crossover but the timing of their annual highs and lows converge. Those seven models all have a new maximum in summer, with four also retaining the spring maximum; the two remaining CMIP5 models fail to converge and have only a spring maximum (Fig. 2). For CMIP6, all models reach positive $p'_{CO_2}$ in summer, with 16 crossing over at atmospheric $CO_2$ levels between 365 ppm and 959 ppm. Models with earlier crossover generally have higher equilibrium climate sensitivity (the global-mean change in surface atmospheric temperature eventually reached after atmospheric $CO_2$ is doubled) and higher area-averaged, summer maximum sea surface temperature in 2091–2100 (Supplementary Table 2).

The tendency for $p'_{CO_2}$ to reverse sign in summer during this century is particularly pronounced in the Arctic shelf seas (Fig. 4). Despite their diversity, most CMIP5 models exhibit this sign reversal over most of the shelf seas, which make up about half of the Arctic Ocean's surface area (Supplementary Fig. 6). For example, in the Barents Sea, all CMIP models assessed already simulate an annual high in regional mean $p'_{CO_2}$ in summer or they project crossover by 2100 under high-end emissions; the same holds in the Kara and Chukchi seas barring one CMIP5 model (Supplementary Figs. 7 and 8). Models also project enhanced summer warming across the Arctic, with generally more in the shelf seas (Supplementary Fig. 9) where sea ice recedes earlier. Positive $p'_{CO_2}$ occurs in summer with enhanced warming in regions where the date of sea-ice retreat occurs before 1 August, referred to as the late-summer transition, after which the declining net air-to-sea heat flux fails to warm surface waters⁸ (Supplementary Fig. 10).

Yet warming from sea-ice retreat does not act alone on $p'_{CO_2}$. Generally opposed in sign is the non-thermal contribution, which varies among models, for example, owing to different biogeochemical components. It is the balance of the non-thermal contribution with the gradually increasing thermal contribution that determines when summer $p'_{CO_2}$ becomes positive and when it surpasses the winter maximum, both abrupt transitions. Models differ in $p'_{CO_2}$ because non-thermal as well as thermal
<figure><img src="image_6.png" /><figcaption>**Fig. 2** | Splitting and inversion of summer low tends to occur for $p'_{CO_2}$ but not its drivers. a–h, Arctic domain averages for decadal climatologies of 1996–2005 (a,c,e,g) and 2091–2100 (b,d,f,h) are shown for nine CMIP5 models (historical and RCP8.5) for surface ocean $p'_{CO_2}$ (a,b), fractional ice concentration (c,d), NPP (e,f) and surface ocean temperature (g,h). The line colours represent individual models, the black dots represent the model mean and the shaded</figcaption></figure>

region is the uncertainty (±1 s.d., n = 9). The red dashes are for modern observational estimates (gridded data products) for $p_{CO_2}$, ice fraction and sea surface temperature (Methods). Extended Data Fig. 1 shows analogous results from CMIP6 (SSP5-8.5). Models fall into three groups for simulated $p'_{CO_2}$ in 2091–2100 but may share a common evolution pathway (Fig. 3).

contributions differ, as suggested by their large differences in T' and NPP⁵ (Fig. 2 and Extended Data Fig. 1). Although modelled NPP varies greatly and is imperfect, for example, not capturing observation-based estimates of regional differences²² and recent changes in phenology²³,²⁴, and observations cannot tell us whether NPP will continue to increase²⁵ or decline²⁶ during this century, the use of many diverse models lends confidence that they would encompass the real ocean response. Weighing these contributions requires a quantitative framework.

Quantifying drivers

A first framework, three idealized experiments run in two CMIP5 models, suggests that the $p'_{CO_2}$ sign reversal is steered by indirect effects

of physical climate change and amplified by direct effects of higher atmospheric CO₂ on ocean chemistry (Supplementary Section 1). A second framework was used to assess contributions in the non-idealized experiments made with all CMIP5 models. Contributions to variations in $p_{CO_2}$ are often deconvolved term by term with a Taylor expansion¹²

$$p'_{CO_2} \approx \frac{\partial p_{CO_2}}{\partial T} T' + \frac{\partial p_{CO_2}}{\partial S} S' + \frac{\partial p_{CO_2}}{\partial A_T} A_T' + \frac{\partial p_{CO_2}}{\partial C_T} C_T' \quad (1)$$

where the four drivers are temperature T, salinity S, total alkalinity A_T and total dissolved inorganic carbon C_T, the prime indicates a monthly
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

**Fig. 3 | Seasonal timing of p'_{CO₂} evolves gradually in three stages and the annual high and low eventually cross.** a–d, Evolution of CMIP means averaged over the Arctic domain for annual highs and lows (a, b), where the angle is the month and the radius is the atmospheric CO₂ level, and the full seasonal cycle (c, d). a, CMIP5 results (historical + RCP8.5) for p'_{CO₂} (black), sea surface temperature (light blue) and NPP (green). b, CMIP6 p'_{CO₂} in four SSP scenarios. Shown are annual highs (solid) and lows (dashed) and p'_{CO₂} crossover points (large open circles). The low for NPP is indistinct (Fig. 2) and not shown. Encircled points are the p'_{CO₂} data product (Methods). Small circles mark the

end of decades (2100, 2090, 2080 and so on) but end in 2090 for SSP5-8.5. c, d, Evolution of the full seasonal cycle (decadal averages) of p'_{CO₂} in CMIP5 (c) and CMIP6 (d) means occurs in three stages: (1) no maximum in summer (black), (2) positive secondary maximum in summer (light blue) and (3) annual maximum in summer (orange). The thick orange line indicates the first decade reaching stage 3; line patterns distinguish different decades. This evolution pathway is the rule among models but some do not make the full journey (Supplementary Figs. 4 and 5).

anomaly relative to the annual mean, and the partial derivatives are the sensitivities of $p'_{CO_2}$ to each driver. Here this equation was used differently to distinguish the indirect effect of physical climate change (radiative effect) from the direct effect of increasing atmospheric CO₂ on ocean chemistry (geochemical effect), assuming that the former affects the driver anomalies, whereas the latter affects the sensitivities. This climate-CO₂ separation was used to distinguish the combined effects of the sensitivities from those of the driver anomalies in terms of how they affect the total change in $p'_{CO_2}$ between 2006–2015 and 2091–2100 (Methods). For brevity, analysis focused on the nine CMIP5 models. The results confirm that the change in sensitivities increases only the seasonal amplitude, whereas the change in driver anomalies alters the timing (Fig. 5). The effect from the change in the driver anomalies is roughly doubled after accounting for the synergy with the change in sensitivities (Extended Data Fig. 5). That synergy is strong enough in most models to produce a sign reversal in the basin-wide mean $p'_{CO_2}$ during at least part of summer. Thus it is the radiative effect that produces the change in timing.

This timing change also affects the seasonal amplitude and trends in seasonal means. With the timing change, the seasonal amplitude of $p'_{CO_2}$ doubles rather than triples during this century under RCP8.5 (Supplementary section 2). It also causes this century's change in $p'_{CO_2}$ to be 20 ± 7% more in summer and 9 ± 2% less in winter, for seasonal averages, compared

with no change in timing, that is, no physical climate change (Supplementary Section 3).

To further distinguish the individual drivers, we used a freshwater Taylor expansion (equation (4), Methods) in the usual way (Supplementary Section 4). For the modern Arctic Ocean, the broad summer minimum is dominated by the $C_T$ term, as in the subarctic$^{4,12,13,27}$ (Extended Data Fig. 6). Conversely, this century's reversal of summer $p'_{CO_2}$ found in most models is driven by the thermal term.

The growing importance of the thermal term is not driven by the changes in the sensitivities of ocean $p'_{CO_2}$ to its driving variables. Although the sensitivity of $p'_{CO_2}$ to temperature nearly triples, $p'_{CO_2}$ sensitivities to $A_T$ and $C_T$ increase by about 30% more (Supplementary Fig. 11). Hence the main cause must be changes in driver anomalies. Indeed, summer $T'$ roughly triples, whereas the magnitude of salinity normalized $A_T$ ($sA_T'$) changes relatively little and that of $C_T$ ($sC_T'$) even declines (Extended Data Fig. 7). Although the effect of warming generally dominates future seasonal variations of $p'_{CO_2}$, it has little effect on those of other CO₂ system variables except for $[H^+]$.

## Other CO₂ system variables

The annual cycle of $[H^+]'$ and its projected changes in CMIP5 under RCP8.5 are much like those of $p'_{CO_2}$, suggesting a nearly linear relation-
<figure><img src="image_11.png" /><figcaption>**Fig. 4** | Future summer $p'_{CO_2}$ is dominated by warming, particularly in shelf seas. a–i, Arctic maps of the summer anomalies of total $p_{CO_2}$ (a–c), the thermal component (d–f) and the non-thermal component (g–i) are shown for the CMIP5 mean (RCP8.5) as decadal averages for 2006–2015 (modern; a,d,g), 2091–2100 (future; b,e,h) and their difference (c,f,i). The summer anomaly</figcaption></figure>

is the average of the monthly anomalies over the three summer months (June, August, September). The components are from the Taylor expansion. The non-thermal component can be further decomposed into its various contributions, as discussed later, showing for instance that its reduction on the shelves is mostly from reduced influence of $C_T$.

ship despite large changes in temperature (Extended Data Figs. 8 and 9). To clarify their ties, we combine Henry's law ($[CO_2^*] = K_0 C_f p_{CO_2}^{28}$, with the first dissociation constant for carbonic acid ($K_1 = [H^+][HCO_3^-]/[CO_2^*]$) and rearrange to

$$p_{CO_2} = \frac{[H^+][HCO_3^-]}{C_f K_0 K_1}, \qquad (2)$$

where $[HCO_3^-]$ is the bicarbonate ion concentration, $K_0$ is the $CO_2$ solubility, $C_f$ is the fugacity coefficient (which remains very near unity), and $[CO_2^*]$ is the sum of aqueous $[CO_2]$ and $[H_2CO_3]$. At constant temperature, the relationship between $p_{CO_2}$ and $[H^+]$ is nearly proportional because the denominator in equation (2) is a constant, and $[HCO_3^-]$, the other term in the numerator, varies little relative to its large background concentration. But even as temperature changes, the relationship between $p_{CO_2}$ and $[H^+]$ remains nearly linear because thermally driven changes in $[HCO_3^-]$ are negligible, while those for $K_0$ and $K_1$ largely cancel (Extended Data Fig. 10). Thus similar to $p_{CO_2}$, the projected change in seasonal timing of $[H^+]$ over this century enhances the summer basin-wide mean, by 15 ± 6% for that season's mean (Supplementary Table 3) and 23 ± 8% for its extreme.

Conversely, changes in seasonal variations of $[CO_2^*]$ do not resemble those for ocean $p_{CO_2}$ (Extended Data Figs. 8 and 9) because of its weaker sensitivity to changes in temperature. For a closed system with no air-sea $CO_2$ exchange, the relative sensitivity of $p_{CO_2}$ to temperature is currently about seven times greater than that of $[CO_2^*]$ (Extended Data Fig. 10), as determined by the solubility $K_0$. As surface water warms, the $CO_2$ solubility changes instantaneously as does ocean $p_{CO_2}$.

Yet even in an open system, the effect of rapid warming between summer months cannot be fully compensated by air-sea $CO_2$ equilibration,

which has an average e-folding time of several months under current Arctic conditions. Thus ocean $p_{CO_2}$ must change more quickly than does ocean $[CO_2^*]$, which can further adjust only through exchange with the atmosphere. These contrasting tendencies are confirmed by Taylor-series expansions (Supplementary Section 4). Seasonal temperature variations largely affect $p'_{CO_2}$ and $[H^+]$, and thus pH, but much less $[CO_2^*]$ and $[CO_3^{2-}]$, whose seasonal variations are controlled largely by variations in $C_T$ and $A_T$ (Extended Data Fig. 6).

## Biological impacts

Ocean acidification is expected to adversely affect life-sustaining processes across various groups of marine organisms$^{29,30}$. In the Arctic Ocean, marine calcifying organisms may be the most sensitive, particularly pteropods, which already show signs of shell degradation$^{31}$. However, some non-calcifiers also have vulnerable life stages, including those of two keystone species, the copepod *Calanus glacialis*$^{32}$ and polar cod (*Boreogadus saida*)$^{33}$, the main link between zooplankton and marine mammals, seabirds and other fish$^{34}$. These organisms are sensitive to different $CO_2$ system variables$^{35}$.

In marine calcifiers, dissolution of $CaCO_3$ depends on $[CO_3^{2-}]$, whereas $CaCO_3$ production may be only correlated with that variable. That is, many marine calcifiers appear to take up $HCO_3^-$, suggesting that increasing $[HCO_3^-]$ from ocean acidification would stimulate calcification ($Ca^{2+} + HCO_3^- \rightarrow CaCO_3 + H^+$) provided the proton is removed from body fluids by acid-base regulation; conversely, increases in body fluid $[H^+]$ would inhibit calcification$^{36-38}$. Thus the $[HCO_3^-]/[H^+]$ ratio in seawater and body fluids may be more physiologically relevant from a thermodynamic perspective; however, kinetic constraints still imply that $[CO_3^{2-}]$ is a fundamental controlling variable$^{39}$.
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 5** | Changes in driver anomalies enhance extreme summer $p'_{CO_2}$ by 29 ± 9%. a,b, Modern and future monthly anomalies of $p_{CO_2}$ from altered sensitivities and driver anomalies are shown in terms of absolute changes (a) and percent changes relative to the modern state (b). The black lines show the CMIP5 mean in 2006–2015 (modern; solid) and 2091–2100 (future; dashed) in terms of monthly anomalies relative to the annual mean; the coloured lines indicate the changes between those decades from changing sensitivities (red dots) and changing driver anomalies including the synergy term (blue dot-dash). Shading indicates model agreement (±1 s.d., *n* = 9) but is omitted for the two black curves in a, for which it was already shown in Fig. 2.</figcaption></figure>

Projected lower $[CO_3^{2-}]'$ and higher $[H^+]'$ in summer are proportionally much larger than the increase in $[HCO_3^-]'$, relative to background levels (Extended Data Fig. 8), suggesting that the long-term decline in calcification in the Arctic will be sharper in summer months. In some calcifiers, high $[CO_2^*]$ may also mediate the impacts of ocean acidification$^{40}$.

For water-breathing animals, metabolic CO₂ is produced in mitochondria and diffuses across membranes and epithelia until it is eventually released to seawater via the gills, but that efflux could slow under higher environmental CO₂ levels until internal levels adjust$^{41}$. Cross-gill transport is usually formulated in terms of the $p_{CO_2}$ gradient based on Fick's first law. Thus more positive summer $p'_{CO_2}$ during this century would result in further build-up of internal CO₂ as needed to re-establish the gradient across the gills. Fish can partially compensate for such build-up by adjusting the chemistry of interior fluids, but that compensation is less efficient for larvae and young juveniles$^{41}$. This analysis would benefit from also considering effects of warming on internal $p_{CO_2}$ and the $[CO_2^*]$ gradient, acid excretion across gills$^{41-43}$ and respiratory plasticity$^{44}$.

Ocean acidification is concerning, but it is warming that is likely to most affect Arctic marine organisms during this century. Between

1996–2005 and 2091–2100, the summer maximum surface temperature averaged over the Arctic shelf seas is projected to increase from 3.0 °C (after model debiasing) to 8.5 ± 1.2 °C in CMIP5 (RCP8.5, *n* = 9) and to 10.8 ± 2.7 °C in CMIP6 (SSP5-8.5, *n* = 18), beyond the upper thermal limits of certain fish and other fauna endemic to the Arctic$^{33,45,46}$. Thermal vulnerability will be heightened further by simultaneous ocean acidification, particularly in summer, based on observations and experiments with boreal to Arctic crustaceans, bivalves and fish$^{33,46-50}$.

Today, surface $p_{CO_2}$ in the Arctic Ocean generally exhibits a broad summer minimum driven by the dominance of biological drawdown of carbon over warming, in both observations and models. Conversely, local observations from an Arctic polynya region exhibit an increase in $p_{CO_2}$ of 110 µatm in early summer, prompted by the disappearance of seasonal ice cover and, among other factors, a 9 °C warming$^{19}$. Thermally driven increases in summer $p_{CO_2}$ are projected here to amplify and their dominance to become more widespread as atmospheric CO₂ increases, seasonal sea ice recedes, and sea surface temperature rises much more in summer than winter, the opposite of surface atmospheric temperature$^{7}$. Resulting changes in the seasonal amplitude and timing of $p_{CO_2}$ will modulate the change expected from the trend in the mean state. The timing changes projected for ocean $p_{CO_2}$ and $[H^+]$ cause their summer extremes to transition from annual lows to annual highs, enhancing their changes during this century by about one-fourth. Hence, in summer, when biological activity is greater and there is usually a seasonal reprieve, Arctic Ocean acidification may instead be elevated beyond the already large projected change in the mean state, thus also making organisms more vulnerable to heightened temperature.

## Online content

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-022-05205-y.

1. Landschützer, P., Gruber, N., Bakker, D. C. E., Stemmler, I. & Six, K. D. Strengthening seasonal marine CO₂ variations due to increasing atmospheric CO₂. *Nat. Clim. Change* **8**, 146–150 (2018).

2. McNeil, B. I. & Sasse, T. P. Future ocean hypercapnia driven by anthropogenic amplification of the natural CO₂ cycle. *Nature* **529**, 383–386 (2016).

3. Kwiatkowski, L. & Orr, J. C. Diverging seasonal extremes for ocean acidification during the twenty-first century. *Nat. Clim. Change* **8**, 141–145 (2018).

4. Gallego, M. A., Timmermann, A., Friedrich, T. & Zeebe, R. E. Drivers of future seasonal cycle changes in oceanic $p_{CO_2}$. *Biogeosciences* **15**, 5315–5327 (2018).

5. Kwiatkowski, L. et al. Twenty-first century ocean warming, acidification, deoxygenation, and upper-ocean nutrient and primary production decline from CMIP6 model projections. *Biogeosciences* **17**, 3439–3470 (2020).

6. Landschützer, P., Laruelle, G. G., Roobaert, A. & Regnier, P. A uniform $p_{CO_2}$ climatology combining open and coastal oceans. *Earth Syst. Sci. Data.* **12**, 2537–2553 (2020).

7. Carton, J. A., Ding, Y. & Arrigo, K. R. The seasonal cycle of the Arctic Ocean under climate change. *Geophys. Res. Lett.* **42**, 7681–7686 (2015).

8. Steele, M. & Dickinson, S. The phenology of Arctic Ocean surface warming. *J. Geophys. Res. Oceans* **121**, 6847–6861 (2016).

9. Eggleston, E., Sabine, C. L. & Morel, F. M. M. Revelle revisited: buffer factors that quantify the response of ocean chemistry to changes in DIC and alkalinity. *Glob. Biogeochem. Cycles* **24**, GB1002 (2010).

10. Riebesell, U., Kortzinger, A. & Oschlies, A. Sensitivities of marine carbon fluxes to ocean change. *Proc. Natl Acad. Sci. USA* **106**, 20602–20609 (2009).

11. Rodgers, K. B. et al. A wintertime uptake window for anthropogenic CO₂ in the North Pacific. *Glob. Biogeochem. Cycles* **22**, GB2020 (2008).

12. Takahashi, T., Olafsson, J., Goddard, J. G., Chipman, D. W. & Sutherland, S. C. Seasonal variation of CO₂ and nutrients in the high-latitude surface oceans: A comparative study. *Glob. Biogeochem. Cycles* **7**, 843–878 (1993).

13. Sarmiento, J. L. & Gruber, N. *Ocean Biogeochemical Dynamics* (Princeton Univ. Press, 2006).

14. Kessler, A. & Tjiputra, J. The Southern Ocean as a constraint to reduce uncertainty in future ocean carbon sinks. *Earth Sys. Dyn.* **7**, 295–312 (2016).

15. Mongwe, N. P., Vichi, M. & Monteiro, P. M. S. The seasonal cycle of $p_{CO_2}$ and CO₂ fluxes in the Southern Ocean: diagnosing anomalies in CMIP5 Earth system models. *Biogeosciences* **15**, 2851–2872 (2018).

16. Kahru, M., Brotas, V., Manzano-Sarabia, M. & Mitchell, B. G. Are phytoplankton blooms occurring earlier in the Arctic? *Glob. Change Biol.* **17**, 1733–1739 (2010).
17. DeGrandpre, M. et al. Changes in the Arctic Ocean carbon cycle with diminishing ice cover. *Geophys. Res. Lett.* **47**, e2020GL088051 (2020).

18. Ouyang, Z. et al. Sea-ice loss amplifies summertime decadal CO₂ increase in the western Arctic Ocean. *Nat. Clim. Change* **10**, 678–684 (2020).

19. Else, B. G. T. et al. Annual cycle of air–sea CO₂ exchange in an Arctic Polynya Region. *Glob. Biogeochem. Cycles* **27**, 388–398 (2013).

20. Geilfus, N.-X. et al. Spatial and temporal variability of seawater $p_{CO_2}$ within the Canadian Arctic Archipelago and Baffin Bay during the summer and autumn 2011. *Cont. Shelf Res.* **156**, 1–10 (2018).

21. Ahmed, M., Else, B. G. T., Burgers, T. M. & Papakyriakou, T. Variability of surface water $p_{CO_2}$ in the Canadian Arctic Archipelago from 2010 to 2016. *J. Geophys. Res. Oceans* **124**, 1876–1896 (2019).

22. Song, H. et al. Strong and regionally distinct links between ice-retreat timing and phytoplankton production in the Arctic Ocean. *Limnol. Oceanogr.* **66**, 2498–2508 (2021).

23. Ardyna, M. et al. Recent Arctic Ocean sea ice loss triggers novel fall phytoplankton blooms. *Geophys. Res. Lett.* **41**, 6207–6212 (2014).

24. Ardyna, M. & Arrigo, K. R. Phytoplankton dynamics in a changing Arctic Ocean. *Nat. Clim. Change* **10**, 892–903 (2020).

25. Lewis, K. M., van Dijken, G. L. & Arrigo, K. R. Changes in phytoplankton concentration now drive increased Arctic Ocean primary production. *Science* **369**, 198–202 (2020).

26. Zhuang, Y. et al. Freshening leads to a three-decade trend of declining nutrients in the western Arctic Ocean. *Environ. Res. Lett.* **16**, 054047 (2021).

27. Takahashi, T. et al. Climatological distributions of pH, $p_{CO_2}$, total CO₂, alkalinity, and CaCO₃ saturation in the global surface ocean, and temporal changes at selected locations. *Mar. Chem.* **164**, 95–125 (2014).

28. Weiss, R. F. Carbon dioxide in water and seawater: the solubility of a non-ideal gas. *Mar. Chem.* **2**, 203–215 (1974).

29. Wittmann, A. C. & Pörtner, H.-O. Sensitivities of extant animal taxa to ocean acidification. *Nat. Clim. Change* **3**, 995–1001 (2013).

30. Kroeker, K. J. et al. Impacts of ocean acidification on marine organisms: quantifying sensitivities and interaction with warming. *Glob. Change Biol.* **19**, 1884–1896 (2013).

31. Niemi, A. et al. Biological impact of ocean acidification in the Canadian Arctic: widespread severe pteropod shell dissolution in Amundsen Gulf. *Front. Mar. Sci.* **8**, 600184 (2021).

32. Thor, P. et al. Contrasting physiological responses to future ocean acidification among Arctic copepod populations. *Glob. Change Biol.* **24**, e365–e377 (2017).

33. Dahlke, F. T. et al. Northern cod species face spawning habitat losses if global warming exceeds 1.5°C. *Sci. Adv.* **4**, eaas8821 (2018).

34. Christiansen, J. S. No future for Euro-Arctic ocean fishes? *Mar. Ecol. Prog. Ser.* **575**, 217–227 (2017).

35. Fabry, V., Seibel, B., Feely, R. & Orr, J. Impacts of ocean acidification on marine fauna and ecosystem processes. *ICES J. Mar. Sci.* **65**, 414–432 (2008).

36. Bach, L. T. Reconsidering the role of carbonate ion concentration in calcification by marine organisms. *Biogeosciences* **12**, 4939–4951 (2015).

37. Jokiel, P. L. Predicting the impact of ocean acidification on coral reefs: evaluating the assumptions involved. *ICES J. Mar. Sci.* **73**, 550–557 (2016).

38. Cyronak, T., Schulz, K. G. & Jokiel, P. L. The omega myth: what really drives lower calcification rates in an acidifying ocean. *ICES J. Mar. Sci.* **73**, 558–562 (2016).

39. Waldbusser, G. G., Hales, B. & Haley, B. A. Calcium carbonate saturation state: on myths and this or that stories. *ICES J. Mar. Sci.* **73**, 563–568 (2016).

40. Vázquez, V. et al. High-CO₂ levels rather than acidification restrict *Emiliania huxleyi* growth and performance. *Microb. Ecol.* https://doi.org/10.1007/s00248-022-02035-3 (2022).

41. Brauner, C. J. et al. in *Carbon Dioxide* (eds Grosell, M. et al.) 69–132 (Academic Press, 2019); https://doi.org/10.1016/bs.fp.2019.08.003.

42. Esbaugh, A. J., Heuer, R. & Grosell, M. Impacts of ocean acidification on respiratory gas exchange and acid-base balance in a marine teleost, *Opsanus beta*. *J. Comp. Physiol. B* **182**, 921–934 (2012).

43. Tseng, Y.-C. et al. CO₂-driven seawater acidification differentially affects development and molecular plasticity along life history of fish (*Oryzias latipes*). *Comp. Biochem. Physiol. A* **165**, 119–130 (2013).

44. Esbaugh, A. J., Ern, R., Nordi, W. M. & Johnson, A. S. Respiratory plasticity is insufficient to alleviate blood acid-base disturbances after acclimation to ocean acidification in the estuarine red drum, *Sciaenops ocellatus*. *J. Comp. Physiol. B* **186**, 97–109 (2016).

45. Franklin, C. E., Farrell, A. P., Altimiras, J. & Axelsson, M. Thermal dependence of cardiac function in arctic fish: implications of a warming world. *J. Exp. Biol.* **216**, 4251–4255 (2013).

46. Kunz, K. L., Claireaux, G., Pörtner, H.-O., Knust, R. & Mark, F. C. Aerobic capacities and swimming performance of polar cod (*Boreogadus saida*) under ocean acidification and warming conditions. *J. Exp. Biol.* **221**, job184473 (2018).

47. Pörtner, H.-O. Oxygen- and capacity-limitation of thermal tolerance: a matrix for integrating climate-related stressor effects in marine ecosystems. *J. Exp. Biol.* **213**, 881–893 (2010).

48. Walther, K., Sartoris, F. J., Bock, C. & Pörtner, H. O. Impact of anthropogenic ocean acidification on thermal tolerance of the spider crab Hyas araneus. *Biogeosciences* **6**, 2207–2215 (2009).

49. Zittier, Z. M. C., Hirse, T. & Pörtner, H.-O. The synergistic effects of increasing temperature and CO₂ levels on activity capacity and acid-base balance in the spider crab, *Hyas araneus*. *Mar. Biol.* **160**, 2049–2062 (2012).

50. Dahlke, F. T. et al. Effects of ocean acidification increase embryonic sensitivity to thermal extremes in Atlantic cod, *Gadus morhua*. *Glob. Change Biol.* **23**, 1499–1510 (2016).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_14.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022
# Methods

## Arctic domain

Following previous studies<sup>7,51</sup>, we adopt the definition of the Arctic Ocean domain as being bounded by four Arctic gateways: the Barents Sea Opening and the Fram, Davis and Bering straits. This domain excludes the Nordic seas, which remain largely ice free even in winter. Ocean grid points external to the domain are masked out, both when showing maps (masking in dark grey) and when computing integrated quantities such as averages over the Arctic domain (basin-wide mean) or the Arctic's individual regional seas.

## Earth system models

To assess potential changes while accounting for regional differences, physical climate change, and the carbon cycle, we used a suite of 9 ESMs that participated in CMIP5 and 18 ESMs that participated in CMIP6 (Supplementary Table 1). Using models from both phases improves statistical robustness and takes advantage of improvements in the community of models over the past decade<sup>52</sup> while providing a test to check whether conclusions hold across model generations. One improvement in CMIP6 that could be important for the Arctic is that some models have much finer lateral resolution.

The CMIP models were used to assess monthly variations and trends in surface ocean $p_{CO_2}$ and related surface-water ocean $CO_2$ system variables. Monthly means of CMIP5 $CO_2$ system variables were previously computed³ from monthly mean model output for total dissolved inorganic carbon $C_T$, total alkalinity $A_T$, temperature $T$, salinity $S$, total dissolved inorganic phosphorus $P_T$ and total dissolved silicon $Si_T$ using mocsy⁵³. For CMIP6, the only $CO_2$ system variable analysed was surface ocean $p_{CO_2}$, which was provided by each model group. The CMIP5 results for 1860–2005 (historical experiment) were combined with three projections for 2006–2100 under the RCPs that reach radiative forcings of 2.6 W m⁻², 4.5 W m⁻² and 8.5 W m⁻² (rcp26, rcp45 and rcp85 experiments), referred to as RCP2.6, RCP4.5 and RCP8.5, respectively. Likewise, the CMIP6 results for 1850–2014 (historical experiment) were combined with those for four projections for 2015–2100 under the Shared Socioeconomic Pathways (SSPs) that reach radiative forcings of 2.6 W m⁻², 4.5 W m⁻², 7.0 W m⁻² and 8.5 W m⁻² (ssp126, ssp245, ssp370 and ssp585 experiments), referred to as SSP1-2.6, SSP2-4.5, SSP3-7.0 and SSP5-8.5, respectively.

Fields were regridded to a 1° × 1° regular grid for model evaluation and comparison. Monthly mean anomalies relative to the annual mean were computed by subtracting a cubic-spline fit⁵⁴ of the monthly mean times series at each grid cell. Decadal mean climatologies were compared between either 1996–2005 and 2091–2100 or 2006–2015 and 2091–2100, and tendencies were assessed as a function of the atmospheric $CO_2$ level. The driving mechanisms were assessed with (1) Taylor-series expansions to quantify contributions from each of the four main input variables ($C_T$, $A_T$, $T$ and $S$) and (2) idealized scenarios from CMIP5 with multiple simulations under different forcing to separate the direct chemical consequences attributable to the increase in atmospheric $CO_2$ (geochemical effect) from the indirect consequences of physical climate change (radiative effect). Error bars given in the text are reported as ±1 s.d about the multimodel mean.

## Model evaluation

CMIP5 and CMIP6 seasonal climatologies constructed from the 1996–2005 model output years of each historical experiment were compared over the Arctic Ocean domain to observation-based products of surface ocean $p_{CO_2}$ (refs. ⁶,⁵⁵), sea surface temperature⁵⁶ and sea-ice concentration⁵⁷. Both the $p_{CO_2}$ and sea-ice concentration⁵⁷ data products are provided on a 0.25° × 0.25° regular latitude–longitude grid. For comparison, these data products and all model fields were regridded to the same World Ocean Atlas 1° × 1° regular latitude–longitude grid, that is, that of the sea-surface-temperature

data product. Regridding was done using the nearest-neighbour algorithm in the cdo package (cdo remapnn)⁵⁸. A common land mask was applied from the World Ocean Atlas over the Arctic Ocean domain.

## Debiasing

Models were debiased only to compute the maximum summer temperature in the shelf seas for 2091–2100 under RCP8.5 as well as under SSP5-8.5. The individual models were debiased by computing the 1996–2005 climatology for each model and then subtracting from that the observational climatology⁵⁶ to obtain the model bias. The bias, at each grid point and month, was then removed from the 2091–2100 climatology. Then the annual maxima were computed and the area-weighted average taken over all grid points in the Arctic shelf seas (bottom depths <500 m). Results are reported as the CMIP5 and CMIP6 multimodel means ±1 s.d. In the same way, the observational database was masked and the maxima computed to obtain the modern data-based reference for the area-weighted summer maximum for the Arctic shelf seas. For other analyses, models were not debiased.

## Idealized experiments

Three of the CMIP5 models each provided a set of three idealized experiments: 1pctCO₂, esmFixClim1 and esmFdbk1. All three experiments are forced by atmospheric $CO_2$ that increases at the same rate, 1% per year (doubling after 70 years and quadrupling after 140 years, both relative to the pre-industrial level), but how that is felt by the Earth system differs. The 1pctCO₂ simulation considers both the direct ‘geochemical’ effect of increasing $CO_2$ on the carbon cycle and the ‘radiative’ effect of $CO_2$ on physical climate, which drives physical changes, thereby affecting the carbon cycle indirectly. The esmFixClim1 simulation has identical forcing but considers only the direct effect (geochemical), whereas the esmFdbk1 simulation considers only the indirect effect (radiative).

These CMIP5 idealized experiments allow one to deconvolve the geochemical and radiative contributions, but they also come with limitations. The most obvious is that the rate of increase in atmospheric $CO_2$ is larger than in the historical and high-end RCP8.5 scenarios. Second, the separation between the three experiments is imperfect. In the esmFixClim1 simulation, intended to eliminate physical climate change effects, there are slight increases in ocean temperature linked to the response of the terrestrial biosphere, because the higher $CO_2$ reduces stomatal conductance, thus causing greater surface fluxes of sensible heat. Third, only three models provided results for the full set of experiments, and only two of those (IPSL-CM5A-LR and MPI-ESM-LR) continued the 1% per year atmospheric $CO_2$ increase until the end of the 140-year simulation, whereupon atmospheric $CO_2$ had quadrupled relative to the pre-industrial starting point (284 ppm). Conversely, the third model (GFDL-ESM2M) stopped that increase after the first 70 years, the point at which atmospheric $CO_2$ had doubled, holding the same level over the remaining 70 years. Given the limited number of models, no attempt was made to provide quantitative estimates of model uncertainty.

## Climate–$CO_2$ Taylor separation

As $p'_{CO_2} = f(T, S, A_T, C_T)$, a Taylor-series expansion yields equation (1), after neglecting second-order terms, covariances and minor contributions from $P_T$ and $Si_T$. Typically, that equation is used to compare contributions from each of the four drivers¹², but here we use it in a broader way to distinguish the effects of the atmospheric $CO_2$ increase from those of physical climate change. Thus the effects from the combined sensitivities were separated from those of the combined driver anomalies to ascribe general causes for differences between the modern reference state, defined as the 2006–2015 decadal climatology of surface $p'_{CO_2}$, and the future state, defined as the 2091–2100 climatology. Mathematically, that separation takes the following form
$$ \Delta p'_{CO_2} = \left( \sum_{i=1}^{4} \Delta \mathbf{v}_i \cdot \mathbf{X}_{i,0} \right) + \left( \sum_{i=1}^{4} \mathbf{v}_{i,0} \cdot \Delta \mathbf{X}'_i \right) + \left( \sum_{i=1}^{4} \Delta \mathbf{v}_i \cdot \Delta \mathbf{X}'_i \right), \quad (3) $$

where the prime is the monthly anomaly, Δ is the change between the modern (2006–2015) and future (2091–2100) decades, the 0 subscript refers to the modern decade, and two vectors represent the four drivers $\mathbf{X} = (T, S, A_T, C_T)$ and the corresponding sensitivities $\mathbf{Y} = (\partial p_{CO_2}/\partial T, \partial p_{CO_2}/\partial S, \partial p_{CO_2}/\partial A_T, \partial p_{CO_2}/\partial C_T)$. The first right-hand side term (in parentheses) in equation (3) characterizes the effect of increasing atmospheric CO₂ (without physical climate change), which affects the sensitivities, whereas the second and third terms (each in separate parentheses) would be null without the effect of physical climate change, which affects the driver anomalies.

In practice, the normal Taylor expansion for the modern state is made using modern sensitivities and driver anomalies, whereas that for the future state is made using future sensitivities and driver anomalies. The total difference between those two states is due to changes in both the sensitivities and the driver anomalies. To isolate the effect of the changes in the sensitivities (without physical climate change), the sum of the four terms in equation (1) is computed using the future sensitivities with the modern driver anomalies and then the modern reference state is subtracted to get the perturbation (ΔSensitivities, first term in equation (3)). Likewise, the effect of the changes in driver anomalies is computed using the modern sensitivities with the future driver anomalies and then subtracting the modern state (ΔAnomalies, second term in equation (3)). However, the sum of the modern reference state and those two perturbations does not add up to the future state because it does not account for the synergy between the change in driver anomalies and the increase in sensitivities (third term in equation (3)). Thus that synergy is accounted for along with the change in driver anomalies (ΔAnomalies*, second and third terms in equation (3)) by subtracting the state computed with only increased sensitivities (first term) from the future state (all three terms).

This climate–CO₂ separation is a simplification of a more elaborate regrouping of terms that was derived to analyse contributions to the amplitude of the annual cycle of $p_{CO_2}$ in the CMIP5 models in a study⁴ that excluded the Arctic Ocean and did not address seasonal timing. Nor have previous studies emphasized that changes in sensitivities come mainly from the increase in atmospheric CO₂, whereas changes in driver anomalies come from physical climate change. This climate–CO₂ Taylor-series expansion requires results from just one model experiment. Hence, we were able to use it here to assess all nine CMIP5 models forced under the RCP8.5 scenario, unlike the approach described in the previous section that relies on multiple idealized simulations carried out with less realistic forcing and for which only two CMIP5 models have provided a complete set of results.

## Freshwater Taylor-series expansion

To assess the contributions to $p_{CO_2}$ variations from individual terms, a Taylor expansion that accounts for effects from freshwater fluxes is adopted. After normalizing the $A_T$ and $C_T$ terms in equation (1) to a reference salinity $S_0$, ref. ⁵⁹ noted that interannual variations in $p_{CO_2}$ driven by the normalized $A_T$ term became negligible whereas the normalized $C_T$ term declined, essentially becoming equal to the sum of the last two terms in equation (1), which are not normalized. Building on that finding and the work of ref. ⁶⁰, ref. ⁶¹ introduced a modified equation that separated out the effects on $A_T$ and $C_T$ into those that are biogeochemically driven and those that are physically driven from variations in freshwater fluxes (precipitation minus evaporation, river input, and sea-ice melt and formation). Thus equation (1) can be rewritten as

$$ p'_{CO_2} \approx \frac{\partial p_{CO_2}}{\partial T} T' + \frac{\partial p_{CO_2}}{\partial S} S' + \left[ \frac{A_T}{S} \frac{\partial p_{CO_2}}{\partial A_T} S' + \frac{C_T}{S} \frac{\partial p_{CO_2}}{\partial C_T} S' \right] + \left( \frac{S}{S_0} \frac{\partial p_{CO_2}}{\partial A_T} s A_T' + \frac{S}{S_0} \frac{\partial p_{CO_2}}{\partial C_T} s C_T' \right), \quad (4) $$

where $sA_T$ and $sC_T$ are the salinity-normalized quantities ($sX = XS_0/S$). In our case, $S_0$ is the annual mean salinity in each grid cell because the focus is on monthly anomalies relative to the annual mean, a choice also adopted previously¹ that should minimize known problems with salinity normalization⁶². Thus each of the two original terms for $A_T$ and $C_T$ in equation (1) are split into two components: one driven by variations in salinity (freshwater fluxes, in square braces) and another driven by variations in salinity-normalized quantities (biogeochemical, in parentheses).

Many subsequent studies have used this freshwater Taylor expansion. However, the $S/S_0$ ratio before the two terms in parentheses was subsequently dropped⁶³, a simplification that is often adopted¹,⁶⁴⁻⁶⁶. That is, the $S/S_0$ ratio is assumed to be equal to 1, for example, for seasonal anomalies relative to the annual mean¹. Here this simplification is avoided because in the Arctic Ocean, substantial short-term variations in salinity are expected.

In practice, the deconvolution was performed locally and resulting terms were area-weighted for basin-wide averages. The partial derivatives (sensitivities) were computed numerically using derivnum⁶⁷ from mocsy⁵³. For the 2006–2015 mean, we adopted a basic approach: for each term and month, the monthly mean anomaly relative to the annual mean was computed and multiplied by the average of the corresponding monthly mean and annual-mean sensitivities. The sum of all terms generally agreed well with the actual simulated variable (for example, $p_{CO_2}$) for that modern decadal average. Conversely, at the end of the century under RCP8.5 (2091–2100), the basic approach led to poor agreement when there were dramatic changes between months, such as between the $p_{CO_2}$ minimum in early summer and its maximum in late summer. To improve agreement, we revised the approach for the end-of-century deconvolution following three steps: (1) anomalies were instead computed for each month between consecutive years (between Januaries, between Februaries and so on) and multiplied by the corresponding average sensitivity between each pair of years; (2) each of those products (for each month and each term) were then summed up across years to have a decomposition of the total change between the modern and future decades; and (3) finally, the total change for each term and month were added to each term of the monthly deconvolution for 2006–2015 to obtain the deconvolution for 2091–2100. Agreement for that decade then became similar to that found when using the basic approach for 2006–2015.

## Spatial average timing of highs and lows

Plots are shown detailing the evolution in timing of the annual high and low for multiple variables. That timing (month of annual high and low) is represented as an average across a basin or region, which is computed in one of two ways. For monthly anomalies of $p_{CO_2}$ and other variables, the area-weighted mean of the variable was first computed for each month of the annual-cycle decadal climatology, and then the months of the maximum and minimum were selected from the resulting 12 points. This approach gives less weight to regions with low monthly anomalies such as in ice-covered regions for $p'_{CO_2}$. A second approach was used for sea-ice retreat and growth dates, defined as when sea-ice concentration first drops below 0.15 and when it first rises back above 0.15, respectively. In this case, the timing (month index) at the different grid cells was recorded and used to compute the area-weighted mean month index. When showing the evolution of this timing of annual highs and lows as a function of increasing atmospheric CO₂, curves were fit with a cubic spline to suppress interannual variations.
## CO₂ system equilibrium calculations

Equilibrium calculations for [H⁺], pCO₂ and [CO₂*] shown in Extended Data Fig. 10 were made with mocsy⁵³ and the constants recommended for best practices with Aᵣ = 2,130 µmol kg⁻¹ and Cᵣ = 2,000 µmol kg⁻¹. Total dissolved inorganic phosphorus and silicon were set to zero.

## Data availability

Output from the CMIP5 and CMIP6 Earth system models used in this study is available for download from the Earth System Grid Federation, for example, from the node at Lawrence Livermore National Laboratory (https://esgf-node.llnl.gov/projects/esgf-llnl/). The gridded observation product of surface ocean temperature is from the World Ocean Atlas 2018⁵⁶ and is available at https://www.ncei.noaa.gov/products/world-ocean-atlas. The gridded observation product for sea-ice concentration from the National Snow and Ice Data Center (NSIDC), referred to as Gridded Monthly Sea Ice Extent and Concentration, 1850 Onward, Version 2 from ref. ⁵⁷, is available at https://nsidc.org/data/g10010. The gridded observation product for surface ocean pCO₂ from refs. ⁶,⁵⁵ is available at https://accession.nodc.noaa.gov/0209633 and that from ref. ², which was used only in Supplementary Fig. 1b, is taken from their supplementary information.

## Code availability

For analysis, we used open-source Python, versions 2.7 and 3, namely the numpy, sciPy, pandas and matplotlib modules. To regrid model data, we used climate data operators available at https://code.mpimet.mpg.de/projects/cdo/. To compute ocean CO₂ system variables and their sensitivities to driving variables, we used mocsy⁵³,⁶⁷, available from https://github.com/jamesorr/mocsy. To detrend and smooth data, we relied on a fast cubic spline⁵⁴ translated for use in Python and available at https://github.com/jamesorr/fastspline.

51. Tsubouchi, T. et al. The Arctic Ocean in summer: a quasi-synoptic inverse estimate of boundary fluxes and water mass transformation. *J. Geophys. Res. Oceans* **117**, C01024 (2012).

52. Séférán, R. et al. Tracking improvement in simulated marine biogeochemistry between CMIP5 and CMIP6. *Curr. Clim. Change Rep.* **6**, 95–119 (2020).

53. Orr, J. C. & Epitalon, J.-M. Improved routines to model the ocean carbonate system: mocsy 2.0. *Geosci. Model Dev.* **8**, 485–499 (2015).

54. Weinert, H. L. A fast compact algorithm for cubic spline smoothing. *Comput. Stat. Data Anal.* **53**, 932–940 (2009).

55. Landschützer, P., Laruelle, G. G., Roobaert, A. & Regnier, P. A combined global ocean pCO₂ climatology combining open ocean and coastal areas (NCEI Accession 0209633). *NCEI* https://doi.org/10.25921/qb25-f418 (2020).

56. Locarnini, M. et al. *World Ocean Atlas 2018, Volume 1: Temperature Atlas NESDIS 81* (NOAA, 2018).

57. Walsh, J., Chapman, W. & Fetterer, F. *Gridded Monthly Sea Ice Extent and Concentration, 1850 Onward, Version 2* (2019); http://nsidc.org/data/g10010

58. Schulzweida, U. CDO User Guide (Version 1.9.8). Zenodo https://doi.org/10.5281/zenodo.3539275 (2019).

59. Wetzel, P., Winguth, A. & Maier-Reimer, E. Sea-to-air CO₂ flux from 1948 to 2003: a model study. *Glob. Biogeochem. Cycles* **19**, GB2005 (2005).

60. Keeling, C. D., Brix, H. & Gruber, N. Seasonal and long-term dynamics of the upper ocean carbon cycle at station ALOHA near Hawaii. *Glob. Biogeochem. Cycles* **18**, GB4006 (2004).

61. Lovenduski, N. S., Gruber, N., Doney, S. C. & Lima, I. D. Enhanced CO₂ outgassing in the Southern Ocean from a positive phase of the Southern Annular Mode. *Glob. Biogeochem. Cycles* **21**, GB2026 (2007).

62. Friis, K., Körtzinger, A. & Wallace, D. W. R. The salinity normalization of marine inorganic carbon chemistry data. *Geophys. Res. Lett.* **30**, 1085 (2003).

63. Doney, S. C. et al. Mechanisms governing interannual variability in upper-ocean inorganic carbon system and air-sea CO₂ fluxes: Physical climate and atmospheric dust. *Deep Sea Res. Part II* **56**, 640–655 (2009).

64. Hauri, C. et al. Spatiotemporal variability and long-term trends of ocean acidification in the California Current System. *Biogeosciences* **10**, 193–216 (2013).

65. Turi, G., Lachkar, Z. & Gruber, N. Spatiotemporal variability and drivers of pCO₂ and air-sea CO₂ fluxes in the California Current System: an eddy-resolving modeling study. *Biogeosciences* **11**, 671–690 (2014).

66. Arruda, R. et al. Air-sea CO₂ fluxes and the controls on ocean surface pCO₂ seasonal variability in the coastal and open-ocean southwestern Atlantic Ocean: a modeling study. *Biogeosciences* **12**, 5793–5809 (2015).

67. Orr, J. C., Epitalon, J.-M., Dickson, A. G. & Gattuso, J.-P. Routine uncertainty propagation for the marine carbon dioxide system. *Mar. Chem.* **207**, 84–107 (2018).

**Acknowledgements** The World Climate Research Programme's Working Group on Coupled Modelling is responsible for the Coupled Model Intercomparison Project (CMIP). The US Department of Energy's Program for Climate Model Diagnosis and Intercomparison provided coordinating support for CMIP and led the development of its software infrastructure in partnership with the Global Organization for Earth System Science Portals. We thank the modelling centres that participated in CMIP5 and CMIP6 and provided the model output used here. That output is freely available from the Earth System Grid Federation (ESGF) (https://esgf-node.llnl.gov/search/cmip6/). We thank O. Torres for local organization and restructuring of the CMIP6 data. Support for this analysis came from the EU H2020 projects COMFORT (grant 820989) and 4C (grant 821003), the French Agence Nationale de la Recherche (ANR) projects SOBUMS (grant ANR-16-CE01-0014) and CONVINCE (grant ANR-18-ERC2-0001-0), and the Deutsche Forschungsgemeinschaft (Po 278/16-1 and -2) as part of the Research Unit Tersane (FO 2332). Analysis was carried out with computational resources from the IPSL Prodiguer-Ciclad cluster, a facility supported by CNRS, UPMC and Labex L-IPSL, the latter which is funded by the ANR (grant ANR-10-LABX-0018) and by the European FP7 IS-ENES2 project (grant 312979), as well as resources from TGCC and IDRIS through a GENCI/DARI grant (gen0040).

**Author contributions** The study arose from discussions between J.C.O. and L.K. J.C.O did the analysis and initial writing, and all authors contributed to scientific discussions and revisions.

**Competing interests** The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41586-022-05205-y.

**Correspondence and requests for materials** should be addressed to James C. Orr.

**Peer review information** Nature thanks Claudine Hauri, Zhangxian Ouyang and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permissions information** is available at http://www.nature.com/reprints.
<figure><img src="image_15.png" /></figure>

**Extended Data Fig. 1** | CMIP6 also reveals splitting and inversion of the summer low for *p*'$_{CO_2}$, but not its drivers. Monthly variations averaged over the Arctic domain for decadal climatologies of 1996–2005 (a, c, e, g) and 2091–2100 (b, d, f, h) are shown for 18 CMIP6 models (historical and SSP5-8.5) for a, b, surface ocean *p*'$_{CO_2}$, c–d, fractional ice concentration, e, f, NPP, and

g, h, surface ocean T. Line colours represent individual models and black dots
the model mean, while the shaded region is the uncertainty (±1 s.d., n = 18). Red
dashes are for modern observational estimates (gridded data products) for
p<sub>CO2</sub> ice fraction, and SST. Fig. 2 shows analogous results from CMIP5 (RCP8.5).
<figure><img src="image_16.png" /><figcaption>Extended Data Fig. 2 | Summer *p*'$_{CO₂}$ turns positive in most models under mid-range to high-end emissions. All panels show climatological averages of surface ocean *p*'$_{CO₂}$, averaged over the Arctic domain. Columns are for CMIP5 (a, c, e, g) and CMIP6 (b, d, f, h). Rows are for a, b, 1996–2005 from the historical</figcaption></figure>

experiments and for 2091–2100 from the RCP and SSP scenarios at **c, d, 2.6 e, f, 4.5**, and **g, h, 8.5 W m⁻²** (see Methods). Line colours represent individual models and black dots the multimodel mean, while the shaded region is the uncertainty (±1 s.d., n = 9 for CMIP5, n = 18 for CMIP6).
<figure><img src="image_17.png" /></figure>

**Extended Data Fig. 3** | Most models exhibit crossover but at varying levels of atmospheric CO₂. Evolution of timing of annual highs and lows in *p*<sub>CO₂</sub> averaged over the Arctic domain, combining historical and high-end scenarios. Columns are for the **a**, **c**, annual low and **b**, **d**, annual high; rows are for **a**, **b**, CMIP5 RCP8.5 and **c**, **d**, CMIP6 SSP5-8.5 models. In **c**, IPSL-CM6A-LR

and MRI-ESM2-0 are not shown due to noise from averaging over a broad minimum spanning January (Dec-Mar). Plot characteristics are as in Fig. 3. Italic model labels indicate thin lines. Red points are the observation-based product for *p*<sub>CO₂</sub>. Small circles mark the end of decades (2100, 2090, 2080 ...), but 2090 is reached in SSP5-8.5 at the 1000-ppm limit of the plot domain.
<figure><img src="image_18.png" /></figure>

**Extended Data Fig. 4** | Drivers of $p_{CO_2}$ do not exhibit crossover in any CMIP model despite large changes in sea ice. Seasonal evolution of $p_{CO_2}$-related drivers in the CMIP models for their timing of a, d, sea-ice retreat and regrowth (see Methods), b, e, annual high of net primary production, and c, f, annual low

and high of temperature are shown for both a-c, CMIP5 RCP8.5 and d-f, CMIP6 SSP5-8.5. Curves indicate Arctic domain averages for each model, while red points are for the observation-based product of sea-ice concentration (Methods). Plot characteristics are as in Fig. 3 and Extended Data Fig. 3.
<figure><img src="image_19.png" /><figcaption>**Extended Data Fig. 5** | Summer *p*'₂₀₅ₐ increases because effects from changing driver anomalies outweigh those from changing sensitivities.</figcaption></figure>

Climatologies are shown for each CMIP5 model under RCP8.5 for *p*'₂₀₅ₐ in 2006–2015 (modern, black solid) and 2091–2100 (future, black dashes) as well as contributions to the change between the two periods, assessed with the climate–CO₂ Taylor separation, from changing sensitivities (red dots) and

changing driver anomalies alone (cyan dot-dot-dash). The latter is also shown combined with the synergy term, i.e. its amplification from the change in sensitivities (dark blue dash-dots). The dashed black curve for *p*'₂₀₅ₐ in 2091–2100 is the sum of the black solid, red dotted, and dark blue dash-dotted curves (equation 3).
<figure><img src="image_20.png" /></figure>

**Extended Data Fig. 6** | Thermally driven variations increase and become dominant for *p*'<sub>CO₂</sub>, unlike for [CO₂*]⁺ and [CO₃²⁻]*. Contributions to seasonal variability deconvolved with the freshwater Taylor expansion are shown as climatologies of the CMIP5 mean under RCP8.5 for 2006–2015 (left), 2091–2100 (middle), and their difference (right). Modern monthly variations in *p*'<sub>CO₂</sub> (top) are driven mostly by changes in *C*<sub>T</sub>, while in the future, thermally driven

variations become dominant. Conversely, thermally driven variations are always small for [CO₂*]⁺ (middle) and [CO₃²⁻]⁻ (bottom), both now and in the future. Agreement between the simulated variable (black dots, True) and the sum of the components (solid black) indicates the consistency of the deconvolution.
<figure><img src="image_21.png" /></figure>

**Extended Data Fig. 7** | Reversal in summer *p*'_{CO₂} is driven by a threefold increase in summer T', while other changes are much smaller. Arctic domain averages of monthly anomalies of the CMIP5 mean under RCP8.5 are shown for the modern (2006–2015, solid) and future (2091–2100, dashed) climatologies

for the four driving variables: sA'_{T} (cyan), sC'_{T} (green), T' (blue), and S' (red). The shading indicates model spread (±1 s.d., n = 9), which for sA'_{T} and sC'_{T} is shown only for the future.
<figure><img src="image_22.png" /><figcaption>**Extended Data Fig. 8** | Changes in seasonal variations of $p_{CO_2}$ and $[^{13}C]CO_2$ are similar but differ from other CO$_2$ system variables. Climatologies of monthly anomalies for the CMIP5 mean under RCP8.5 averaged over the Arctic domain</figcaption></figure>

are shown for $p_{CO_2}$ and other surface CO$_2$ system variables for 2006–2015 (modern, blue solid) and 2091–2100 (future, orange dashed).
<figure><img src="image_23.png" /></figure>

**Extended Data Fig. 9** | Summer [H⁺]⁻ reverses sign during this century mainly in the shelf seas, resembling *p'CO₂*. Arctic maps of summertime anomalies of [H⁺] and [CO₂*] for the CMIP5 mean under RCP8.5 are shown for climatologies of 2006–2015 (left, modern), 2091–2100 (middle, future), and

their difference (right). The similarity of the sign reversal in summer anomalies of *p*<sub>CO₂</sub> and [H⁺] is explained by their common sensitivity to temperature, while [CO₂*] does not reverse sign because it remains relatively insensitive to temperature.
<figure><img src="image_24.png" /></figure>

**Extended Data Fig. 10** | Seasonality of $p_{CO_2}$ and $[CO_2^{*}]$ differ because of the $CO_2$ solubility $K_0$, which depends on T. Shown is the relative effect of T on $K_0$, $K_1$, and their product as well as on $[H^+]$, $p_{CO_2}$, and $[CO_2^{*}]$. Values between 0 to 10 °C were normalized by dividing by the value at 0 °C to compute their relative

<figure><img src="image_25.png" /></figure>

effect. All $K_0$ and $K_1$ are computed at S = 34. The $CO_2$ system variables were computed with standard equilibrium calculations at the same salinity and temperatures as described in the Methods. Although not shown, $[HCO_3^-]$ declines by 0.3% over the 0–10 °C range.