Drivers of future seasonal cycle changes in oceanic pCO2

Recent observation-based results show that the seasonal amplitude of surface ocean partial pressure of CO2 (pCO2) has been increasing on average at a rate of 2–3 μatm per decade (Landschützer et al., 2018). Future increases in pCO2 seasonality are expected, as marine CO2 concentration ([CO2]) will increase in response to increasing anthropogenic carbon emissions (McNeil and Sasse, 2016). Here we use seven different global coupled atmosphere–ocean– carbon cycle–ecosystem model simulations conducted as part of the Coupled Model Intercomparison Project Phase 5 (CMIP5) to study future projections of the pCO2 annual cycle amplitude and to elucidate the causes of its amplification. We find that for the RCP8.5 emission scenario the seasonal amplitude (climatological maximum minus minimum) of upper ocean pCO2 will increase by a factor of 1.5 to 3 over the next 60–80 years. To understand the drivers and mechanisms that control the pCO2 seasonal amplification we develop a complete analytical Taylor expansion of pCO2 seasonality in terms of its four drivers: dissolved inorganic carbon (DIC), total alkalinity (TA), temperature (T ), and salinity (S). Using this linear approximation we show that the DIC and T terms are the dominant contributors to the total change in pCO2 seasonality. To first order, their future intensification can be traced back to a doubling of the annual mean pCO2, which enhances DIC and alters the ocean carbonate chemistry. Regional differences in the projected seasonal cycle amplitude are generated by spatially varying sensitivity terms. The subtropical and equatorial regions (40 S–40 N) will experience a ≈ 30–80 μatm increase in seasonal cycle amplitude almost exclusively due to a larger background CO2 concentration that amplifies the T seasonal effect on solubility. This mechanism is further reinforced by an overall increase in the seasonal cycle of T as a result of stronger ocean stratification and a projected shoaling of mean mixed layer depths. The Southern Ocean will experience a seasonal cycle amplification of ≈ 90–120 μatm in response to the mean pCO2-driven change in the mean DIC contribution and to a lesser extent to the T contribution. However, a decrease in the DIC seasonal cycle amplitude somewhat counteracts this regional amplification mechanism.

Abstract. Recent observation-based results show that the seasonal amplitude of surface ocean partial pressure of CO 2 (pCO 2 ) has been increasing on average at a rate of 2-3 µatm per decade (Landschützer et al., 2018). Future increases in pCO 2 seasonality are expected, as marine CO 2 concentration ([CO 2 ]) will increase in response to increasing anthropogenic carbon emissions (McNeil and Sasse, 2016). Here we use seven different global coupled atmosphere-oceancarbon cycle-ecosystem model simulations conducted as part of the Coupled Model Intercomparison Project Phase 5 (CMIP5) to study future projections of the pCO 2 annual cycle amplitude and to elucidate the causes of its amplification. We find that for the RCP8.5 emission scenario the seasonal amplitude (climatological maximum minus minimum) of upper ocean pCO 2 will increase by a factor of 1.5 to 3 over the next 60-80 years. To understand the drivers and mechanisms that control the pCO 2 seasonal amplification we develop a complete analytical Taylor expansion of pCO 2 seasonality in terms of its four drivers: dissolved inorganic carbon (DIC), total alkalinity (TA), temperature (T ), and salinity (S). Using this linear approximation we show that the DIC and T terms are the dominant contributors to the total change in pCO 2 seasonality. To first order, their future intensification can be traced back to a doubling of the annual mean pCO 2 , which enhances DIC and alters the ocean carbonate chemistry. Regional differences in the projected seasonal cycle amplitude are generated by spatially varying sensitivity terms. The subtropical and equatorial regions (40 • S-40 • N) will experience a ≈ 30-80 µatm increase in seasonal cycle amplitude almost exclusively due to a larger background CO 2 concentration that amplifies the T seasonal effect on solubility. This mechanism is further reinforced by an overall increase in the seasonal cycle of T as a result of stronger ocean stratification and a projected shoaling of mean mixed layer depths. The Southern Ocean will experience a seasonal cycle amplification of ≈ 90-120 µatm in response to the mean pCO 2 -driven change in the mean DIC contribution and to a lesser extent to the T contribution. However, a decrease in the DIC seasonal cycle amplitude somewhat counteracts this regional amplification mechanism.

Introduction
Owing to its large chemical capacity to resist changes in [CO 2 ] (referred to as buffering capacity), the ocean has absorbed nearly half of the anthropogenic CO 2 produced by fossil fuel burning and cement production since the industrial revolution (Sabine et al., 2004). While the ocean's absorption of CO 2 lowers the atmospheric concentration, it also increases the ocean's [CO 2 ] and in turn lowers its buffering capacity. This leads to a reduction in the oceanic uptake of CO 2 and an intensification of the pCO 2 seasonal cycle (from now on referred to as δpCO 2 ) (Völker et al., 2002;McNeil and Sasse, 2016). In a recent key observational study by Landschützer et al. (2018), it was demonstrated that the δpCO 2 amplitude increased at a rate of ≈ 2-3 µatm per decade from 1982 to 2015.
Published by Copernicus Publications on behalf of the European Geosciences Union. The pCO 2 already experiences large seasonal fluctuations, which in some regions can reach up to 60 % above and below the annual mean pCO 2 (Takahashi et al., 2002). An intensification of the δpCO 2 amplitude could produce seasonal hypercapnia conditions (McNeil and Sasse, 2016), which together with increased [H + ] seasonality (Hagens and Middelburg, 2016;Kwiatkowski and Orr, 2018) and aragonite undersaturation events (Shaw et al., 2013;Hauri et al., 2015;Sasse et al., 2015) could expose marine life to harmful seawater conditions earlier than expected if considering only annual mean values. Moreover, a projected amplification of δpCO 2 might increase the net CO 2 uptake in some regions, such as the Southern Ocean, thereby further accelerating the decrease in the buffering capacity in that region (Hauck and Völker, 2015).
The pCO 2 seasonal amplitude is controlled mainly by seasonal changes in temperature (T ) and biological activity together with upwelling changes that alter DIC concentrations. Usually, DIC and T changes work in opposite directions (Takahashi et al., 2002;Fay and McKinley, 2017;Sarmiento and Gruber, 2006). In subtropical regions higher pCO 2 values occur in summer when solubility decreases. In subpolar regions, pCO 2 increases in winter when waters upwell that are rich in DIC and when respiration of organic matter takes place. Decreased subpolar pCO 2 occurs in summer when primary productivity is higher and upwelling diminishes. Therefore, we find close relationships of δpCO 2 with the ocean's [CO 2 ], which controls chemical reactions, and with mean pCO 2 , which moderates exchange with the atmosphere. Both factors are related by the solubility constant that depends on temperature and salinity.
Furthermore, the regional differences in the influence of temperature and biology on δpCO 2 are modulated by the ocean's buffering capacity. This is due to the ability of CO 2 to react with seawater to form bicarbonate [HCO  (Zeebe and Wolf-Gladrow, 2001). The larger the buffering capacity, the larger the pCO 2 ability to resist changes in DIC. To quantify this capacity, we can introduce the sensitivity factor γ DIC , which is inversely related to the buffering capacity, defined as γ DIC = ∂ ln(pCO 2 )/∂DIC, (Egleston et al., 2010). Other sensitivity factors are related to the total alkalinity (γ TA ), salinity (γ S ), and temperature (γ T ) changes and are defined in a similar way as ∂ ln(pCO 2 )/∂TA, ∂ ln(pCO 2 )/∂S, and ∂ ln(pCO 2 )/∂T , respectively. It is important to note that the pCO 2 is highly sensitive to temperature due to two factors: first through solubility changes that account for 2 / 3 of the present-day temperature impact, and second through the dissociation constants that control the carbon system reactions (Sarmiento and Gruber, 2006).
While the mechanisms controlling the seasonal cycle of pCO 2 at present day are well documented, the future evolution of these drivers has not been fully elucidated. Current literature suggests that seasonal amplification is a consequence of an increase in the T and DIC contributions to δpCO 2 (Landschützer et al., 2018) and an increased sensitivity of the ocean to these variables (Fassbender et al., 2017).
The aim of our paper is to provide an in-depth analysis of the mechanisms controlling the future strength of δpCO 2 and its regional differences using seven CMIP5 global Earth system models. Our analysis focuses on 21st century evolution using the Representative Concentration Pathway 8.5 (RCP8.5) scenario. We give a comprehensive analysis of the projected evolution of the DIC, TA, T , and S contributions to pCO 2 seasonality. To achieve this goal, we derive explicit analytical expressions for pCO 2 sensitivities in terms of γ DIC , γ TA , γ T , and γ S , thereby extending previous work done by Egleston et al. (2010).

CMIP5 models
For our analysis, pCO 2 , DIC, TA, T , and S monthly mean output variables covering the period from 2006-2100 were obtained from future climate change simulations conducted with seven fully coupled Earth system models that participated in the Coupled Model Intercomparison Project, Phase 5 (CMIP5) (Taylor et al., 2012). The following models were selected based on data availability: CanESM2, CESM1-BGC, GFDL-ESM2M, MPI-ESM-LR, MPI-ESM-MR, HadGEM2-ES, and HadGEM2-CC (see the supplementary material of Hauri et al., 2015). For the purpose of this paper, we used the Representative Concentration Pathway 8.5 (RCP8.5) future climate change simulations (IPCC, 2013). The ocean's surface data sets were regridded onto a 1 • ×1 • grid using climate data operators (CDOs). The Arctic Ocean and the region poleward of 70 • S are removed from the analyses because observational data for model validation are scarce.

Analysis of δpCO 2
To elucidate the underlying dynamical, thermodynamical, biological, and chemical processes controlling δpCO 2 we calculated a first-order Taylor series expansion of δpCO 2 in terms of its four drivers, DIC, TA, T , and S. While T and S are controlled only by physics, DIC and TA are controlled by physical, chemical, and biological processes. Throughout this paper we use salinity-normalized DIC and TA us-ing a mean salinity of 35 psu. This effectively removes the concentration-dilution freshwater effect, following the procedure of Lovenduski et al. (2007). The salinity normalized variables are referred to as DIC s and TA s , corresponding to DIC ·S 0 /S and TA ·S 0 /S, respectively. The freshwater effect on DIC and TA is now included in the S term, renamed as S fw . For the Taylor series expansion, each variable (X = DIC, TA, T , and S) is decomposed into X = X + δX. The term X represents the 21-year mean and δX denotes the seasonal cycle (calculated as the monthly mean deviation from the 21-year average). The Taylor expansion is then computed for an initial  and final (2080-2100) period. We use multi-decade means and eventually multimodel ensemble means to remove the effects of interannual variability. The full first-order series expansion is given by Each term of the right-hand side of Eq. (1) represents the contribution from one of the four drivers of δpCO 2 . The analytical expressions for the derivatives (without the salinity normalization) are given by the following. . The explicit T and S partial derivatives are given in the Supplement (Text S1). The first two derivatives coincide with the results of Egleston et al. (2010) and Hagens and Middelburg (2016), with the exception of the sign of [OH − ] in the Egleston et al. (2010) term S. To verify this approach we compared the sum of the Taylor expansion terms with the full simulated range of δpCO 2 from the model's output. The Taylor expansion reproduces the full seasonal cycle amplitude of the original climate model simulations well (Fig. S1 in the Supplement). The analytical expressions for temperature and salinity presented here are -to our knowledge -the first ones of their kind. Previously the calculation of these terms was based on the approximation given by Takahashi et al. (1993) or on numerical calculations.
To gain more insight into the processes causing the amplification of δpCO 2 we introduce a method based on a second Taylor series expansion described below. Equation (1) can be rewritten using the expressions for the sensitivities γ determined by the relation 1 pCO 2 ∂pCO 2 ∂X = γ X . These sensitivities have been historically used to represent the percentage of change in pCO 2 per unit of DIC, TA, T , or S. With this notation, Eq. (1) can be expressed in the following way.
Each term in Eq. (5) consists of three parts: pCO 2 , the sensitivity γ X , and the corresponding seasonal cycle δX. To understand which component is the main driver for δpCO 2 changes, we perform a second Taylor expansion of the end of the century's δpCO 2 around the initial state of the system in 2006-2026.
To maximize mathematical clarity we will introduce some definitions: first, we introduce the symbol to indicate the difference between the period 2080-2100 and 2006-2026. Therefore, the total future change in δpCO 2 is now referred to as δTpCO 2 . In the same manner, the total changes in sensitivities and seasonal cycles are written as γ DIC s , γ TAs , γ T , γ S fw and δT DIC s , δT TA s , δT T, δTS fw , respectively. Finally, we introduce the vector X formed by the four variables DIC s , TA s , T , and S fw as {X 0 , X 1 , X 2 , X 3 } = {DIC s , TA s , T , S}. With this notation, we can write an expansion of Eq. (5) of the final state of the system by 2080-2100 named X f around the initial state where the first, second, and third terms represent the contributions to δT pCO 2 due to changes in the mean pCO 2 ( pCO 2 ), the pCO 2 sensitivities ( γ X k ), and the seasonal cycles ( δT X k ), respectively; the fourth to sixth rows are the second-order terms. This method is similar to the one used by Landschützer et al. (2018).
3 Results and discussion 3.1 δpCO 2 amplification  (Landschützer et al., 2017). The agreement between models and observations is remarkably good in the equatorial regions, but the initial amplitude is slightly overestimated in the middle and high latitudes (see Fig. S3 in the Supplement). The higher amplitude in models than observations is expected, as the initial period 2006-2026 already experienced an amplification compared to previous years. Moreover, Tjiputra et al. (2014) found that the ocean's pCO 2 historical trend is larger in models than observations when it is estimated in large-scale areas of the ocean. However, they found that model pCO 2 trends agree with observations when the trends are subsampled to the locations where the observations were taken, and therefore they do a good job of reproducing well-known time series. Moreover, differences are expected as Pilcher et al. (2015) suggested that CMIP5 models perform well in reproducing the seasonal cycle timing, but still show considerable errors in reproducing the seasonal amplitude of pCO 2 due to differences in the mechanisms represented in each model, especially in subpolar biomes. By 2080-2100 the annual cycle amplitude attains values of ≈ 197 and ≈101 µatm in the high and middle to low latitudes, respectively (Fig. 1b). These seasonal variations correspond to 20 % and 18 % of annual pCO 2 for the initial and final periods, respectively. Figure 1c shows that the global ocean δpCO 2 will intensify by a factor of 1.5 to 3 for the 2080-2100 period relative to the 2006-2026 reference period. Figure 1d shows the difference in amplitude ( δT pCO 2 ); this pattern differs from the ratio because the ratio overestimates the amplification in areas where the initial amplitude is lower than ≈ 10 µatm. McNeil and Sasse (2016) used observations and a neural-network clustering algorithm to project that by the year 2100, the δpCO 2 amplitude in some regions could be up to 10 times larger than it was in the year 2000. Our mean amplification factor estimation agrees with the mean threefold amplification found for most of the ocean by Mc-Neil and Sasse (2016). However, the high values in this previous study cannot be reproduced here, mainly because we consider 21-year average ratios instead of single-year ratios, which are strongly affected by interannual variability. Using observations, Landschützer et al. (2018) found an increase of 2.2 µatm per decade, which is smaller than our findings of a total 42 µatm increase by the end of the century between 40 • S and 40 • N and a global mean change of 81 µatm on the high latitudes. This difference is again possibly due to the higher mean pCO 2 values in models than observations.
The global ocean mean amplification factor of δpCO 2 roughly coincides with a doubling of pCO 2 (Fig. 2). The direct relationship between these two is explained in Sect. 3.5. Figure 1e-h show the zonal mean of Fig. 1a-d; in general, towards the end of the century the pCO 2 amplifies more in high latitudes, but so does the standard deviation uncertainty among models. This regional pattern agrees with the observation-based findings of Landschützer et al. (2018), which show that high latitudes have already experienced a larger amplification than middle to low latitudes from 1982 to 2015. Furthermore, the same pattern is projected by CMIP5 models for the seasonal amplification of [H + ] by the end of the century (Kwiatkowski and Orr, 2018). This is expected from the near-linear relation between pCO 2 and [H + ]. These regional differences in amplification for pCO 2 can be explained in terms of the relative magnitudes and the phases between the DIC, TA, T , and S contributions, which are explained in subsequent sections.

Present and future drivers of δpCO 2
To understand the driving factors of δpCO 2 and its spatiotemporal differences, we split δpCO 2 into the four different contributions from DIC s , TA s , T , and S fw for the initial and final periods, following Eq. (1). The results are shown in Fig. 3. For most of the ocean, the ensemble mean estimated contributions from DIC s and T to the present-day δpCO 2 are in good agreement with the data-based estimates of Takahashi et al. (2014b) and Landschützer et al. (2017), particularly in the equatorial regions (see Fig. S3 in the Supplement). However, our T and DIC contributions are slightly larger in middle and high latitudes for the same reasons the pCO 2 seasonal amplitude is overestimated (see Sect. 3.1). Also, differences arise between our DIC s contribution and the observation-based so-called "nonthermal" contribution because the nonthermal contribution also includes the total alkalinity and salinity effects. Nonetheless, between 40 • S and 40 • N our ensemble mean shows that δpCO 2 is dominated by changes in temperature that control CO 2 solubility, which decreases in summer, enhancing pCO 2 ; this is in agreement with observations. The Southern Ocean is controlled by DIC, which responds to changes in upwelling and phytoplankton blooms. Both mechanisms act together to decrease (increase) DIC in summer (winter) (Sarmiento and Gruber, 2006).
The models show that the δpCO 2 in the 40 • N to 60 • N band is controlled by T , which disagrees with the abovementioned observations that show a non-temperature dominance in this band. The difference between models and observations arises from two regions: the North Atlantic basin and the northwestern Pacific, specifically near the Oyashio current and the outflows from the Sea of Okhotsk (see Fig. S3 in the Supplement). Most models show a T dominance in the North Atlantic basin; only CESM1-BGC and GFDL-ESM2M show a DIC dominance (see Fig. S4 in the Supplement). The North Atlantic is one of the major sinks of anthropogenic CO 2 ; however, some models fail to estimate its uptake capacity (Goris et al., 2018). Goris et al. (2018) found that models with an efficient carbon sequestration present a DIC-dominated pCO 2 seasonal cycle in the North Atlantic, but models with low anthropogenic uptake show a T dominance in this region. In the northwestern Pacific, Mckinley et al. (2006) found that coarse models are not able to capture the intricate oceanographic features of this area, and therefore the pCO 2 seasonality is not well captured. Towards the end of the century (Fig. 3, right column), the amplification of δpCO 2 is caused by an increase in the DIC s and T contributions and to a lesser extent due to TA s and S fw . Only in the high latitudes does the TA s contribution reinforce the DIC s effect. The δDIC s and δT relative phase and magnitude play an important role in causing regional differences of future δpCO 2 . For example, between 40 and 60 • , we find a lower amplification factor than at 30-40 • in both hemispheres (Fig. 1c), contrary to what we expected from the general observed larger amplification at higher latitudes. In this band of lower amplification, the warm water from subtropical regions meets the nutrient-rich water from the subpolar regions, but the DIC s and T effects are almost 6 months out of phase, and therefore their cancellation is larger than in the 30-40 • latitude band where, for example, in the North Atlantic, there is 9-month phase difference between the two contributions. A clear illustration of this phase effect is found in the Supplement (Fig. S5).
In the Southern Ocean there is a shift in the maximum δpCO 2 occurring from August-September to March-April (Fig. 3, last row). This shift is generated because the T contribution gains importance over DIC s due to a reduction of δDIC s magnitude at the same time that δT increases (Fig. 5). In the equatorial Pacific region (Fig. 5), T dominates over DIC s but both contributions are small due to their low seasonality. Therefore, this region will experience a low amplification in δpCO 2 . In this region some models underestimate the pCO 2 trend (Tjiputra et al., 2014), and therefore the seasonal amplification might be underestimated too. In the following sections we conduct further analysis by decomposing each contribution as the result of three factors: the mean pCO 2 (pCO 2 ), the regional pCO 2 sensitivities (γ DIC , γ TA , γ T , and γ S fw ), and the seasonal cycles (δDIC s , δTA s , δT , and δS fw ) as determined in Eq. (5).

Future pCO 2 sensitivities
The γ DIC and γ TA are projected to increase by the end of the century due to a lower ocean buffering capacity produced by increasing temperature and larger background concentrations of DIC (Fassbender et al., 2017). This agrees with our results shown in Fig. 4, which shows that all regions will experience an increase in γ DIC and γ TA . Lower buffer factors (higher sensitivities factors) are found in regions where DIC and TA have similar values, and they will decrease (increase) as the DIC / TA ratio in the oceans increases (Egleston et al., 2010). The alkalinity sensitivity is negative, as pCO 2 decreases with increasing alkalinity, but we show here the negative γ TA for better comparison. γ TA will increase (with negative values) more than the DIC sensitivity. However, seasonal changes in open-ocean TA s are small, and therefore the total contribution of alkalinity in our analysis is negligible compared to the DIC s and T contributions. γ S fw decreases everywhere except in the Western Pacific Warm Pool. In this region γ S fw increases, probably due to future changes in precipitation that enhance the freshwater effect. In Fig. 4, the sensitivities (γ ) are expressed as a percentage change in pCO 2 per unit in DIC, TA, T , and S, respectively. This follows the approach of Takahashi et al. (1993); however, in their paper the authors compute the Revelle factor, which is related to γ DIC as R = DIC · γ DIC . To illustrate the meaning of the sensitivities, we will focus on the subtropical North Pacific in the 15-40 • N latitudinal band. In this region γ DIC indicates an average 0.6 % change in pCO 2 per unit of DIC in 2006-2026. Therefore, for a δDIC s seasonal cycle amplitude of 40 µmol kg −1 and pCO 2 ≈ 400 µatm, the total δpCO 2 amplitude equals 96 µatm. Following the same reasoning, by 2080-2100, γ DIC increases to 0.7 % and δDIC s decreases to 30 µmol kg −1 ; therefore, for a pCO 2 equal to 800 µatm, the δpCO 2 amplitude due to δDIC amounts to 168 µatm.
Temperature sensitivity has been experimentally determined by Takahashi et al. (1993), who found a value of 0.0423, meaning that pCO 2 changes by about 4 % for ev- ery • C. This value agrees with our global mean ensemble estimate of 0.0428. However, our analytical expression of γ T shows that this value varies regionally and, for reasons unknown to us, it might decrease in the future to a global mean value of 0.0415 (Fig. 4c, third column). The T sensitivity is larger in colder regions and lower in the warmer tropics; however, colder regions will experience a larger reduction in γ T , which locally prevents a larger amplification of the T contribution to δpCO 2 . In the next section we show that the T seasonality is projected to increase in high latitudes, strengthening the T contribution.
3.4 Future δDIC s , δTA s , δT , and δS fw Towards the end of the century, the global mean amplitude of δDIC s is projected to decrease by ≈ 26 %-28 % in the high latitudes (Fig.5a), according to all the CMIP5 Earth system model simulations used here. In the middle to low latitudi-nal band there is no agreement between models; while some show an increase, others project a decrease in amplitude. As suggested by Landschützer et al. (2018), the larger decrease in the Southern Ocean may be the result of changes in the shallow overturning circulation that prevent CO 2 accumulation in this region. This reduction may be counteracted by the predicted increase in productivity owing to a suppression of light and temperature limitations (Steinacher et al., 2010;Bopp et al., 2013).
According to the CMIP5 models, most of the ocean is projected to experience a slight increase in δT , as shown in Fig. 5b. All models show a slight increase in δT ; only one model showed a slightly decrease in the southern region, and two models showed a decrease in the equatorial region during October to December. It is important to note that Fig. 5 shows the seasonal values with the mean T removed. Therefore, when considering the positive T trends, the absolute summer values show an increase and the absolute winter values a decrease. This agrees with the results of Alexander et al. (2018), who showed that models project a seasonal intensification of T with larger warm extremes and reduced cold extremes. The authors attributed the T seasonality intensification to an increased oceanic stratification and an overall shoaling of the mixed layer depth, which confines seasonal changes in a reduced volume of water, producing larger changes at the surface. They also showed that the intensification trends are stronger in summer than winter, as the mixed layer depth is shallower in summer. Moreover, ice-covered regions will experience the largest increase in T seasonality due to the loss of sea ice because the ice melting and freezing moderates the surface water temperature seasonality (Carton et al., 2015).
The TA seasonality is also projected to increase in the high latitudes according to all models, except CESM1-BGC, which shows a decrease. For δS (see Fig. S6 in the Supplement) there is no agreement among the different CMIP5 models, except in the Southern Ocean where all the models show a slight decrease. Kwiatkowski and Orr (2018) demonstrated that the seasonality of the drivers is important to determine future changes in [H + ] seasonality. In the same fashion, our results show that the four δpCO 2 drivers present changes in seasonality, and in particular δDIC s and δT changes are important to explain future projections of the δpCO 2 amplitude. The increase in δT enhances the δpCO 2 amplification, and the reduction of δDIC s in the Southern Ocean locally prevents a larger amplification.

Regional dominant factors
To identify the main cause of the δpCO 2 amplification we use the Taylor series expansion method. With this method we consider the system's final state (δpCO 2 by 2080-2100) as a perturbation of the initial state (δpCO 2 by 2006-2026), as shown in Eq. (6). The expansion is done in three groups of variables: the seasonal cycles of DIC s , TA s , T , and S (δX), the sensitivities of pCO 2 to the same four variables (γ x ), and Black arrows point out that while the T seasonal cycle is projected to increase in most of the ocean, global DIC s is projected to decrease. The shading represents 1 standard deviation across the models. It is important to note that the scale is different for some of the latitudinal bands. the mean pCO 2 (pCO 2 ). Therefore, each term of the expansion represents how much of the total δpCO 2 change (indicated by δT pCO 2 and calculated as the 2080-2100 value minus the 2006-2026 value) is due to the change in each of these factors. We also add the second-order terms that come from their combination. The results are shown in Fig. 6a and they indicate that the leading cause of the δpCO 2 amplification is the change in pCO 2 ( pCO 2 ), which confirms previous findings by Landschützer et al. (2018).
It is important to note that our linear Taylor expansion approach neglects one aspect of the highly nonlinear carbonate chemistry of the ocean: it assumes pCO 2 and the sensitivities as independent variables and therefore does not include the positive feedback between larger pCO 2 and increasing γ DIC (decreasing buffering capacity). Hence in the following, we use changes in pCO 2 and changes in seawater carbonate chemistry synonymously, overall resulting in an enhanced response of δpCO 2 to seasonal changes in DIC, TA, T , and S.
Considering regional differences, we note that amplification increases as we move poleward in spite of decreasing pCO 2 (see Figs. 1 and 2). This characteristic geographical pattern of stronger high-latitude amplification is the result of larger present-day sensitivities (γ DIC s , γ T ) and seasonal am-plitudes (δDIC s , δT ) in the high latitudes that amplify the effect of pCO 2 even when its value is small compared to other regions (see Eq. 6, first row term). Some exceptions can be found south of Greenland and near the subtropical gyres, where pCO 2 reaches higher values and therefore also presents large amplification. We also found spatial differences on smaller scales; for example, the western equatorial Pacific presents lower initial δpCO 2 and amplification than the eastern equatorial Pacific (see Fig. 1). This is because the eastern side of the basin has larger DIC s and T contributions than the western side (see Fig. S2) as a consequence of the upwelling of cold, CO 2 -rich waters in the east, which lower the buffering capacity and induce larger δpCO 2 amplitude due to the seasonal effects of productivity and solubility (Valsala et al., 2014).
To further disentangle which of the two main drivers (DIC s or T ) is most affected by pCO 2 , we decomposed the DIC s and T contributions into their sensitivity, seasonal cycle, and pCO 2 components. Figure 6b shows the total DIC and T components together with the pCO 2 and seasonal cycle effects on them. The effects from the sensitivities are not depicted, as they only play a minor role. Only the γ DIC term gains importance in the Southern Ocean (not shown). The total change in seasonal pCO 2 ( δT pCO 2 ) is depicted as a dashed black. This change is decomposed into changes in seasonalities ( δT X, purple), sensitivities ( γ X , green), mean pCO 2 ( pCO 2 , red), and second-order terms (blue) summed over the four variables that control pCO 2 (DIC, TA, T and S). For comparison with the expansion, δT pCO 2 is calculated from model output (yellow). Column (b) shows the total change in DIC (dashed red) and T (dashed blue) contributions. Also shown are two components of the total change in these contributions: the pCO 2 effect on the DIC (solid orange) and T (solid blue) contributions and the δT DIC (yellow) and δT T (light blue) effects. In column (a), the δpCO 2 change follows the pCO 2 effect. Column (b) shows that, actually, the leading cause of amplification is the pCO 2 effect on the T contribution. It is important note the different scale between column (a) and (b). Also, the scale was reduced in the 15 • S-15 • N region to highlight its features.
In most of the ocean, the pCO 2 effect on T contribution is the leading cause of amplification. This effect is the result of seasonal solubility changes acting over a larger [CO 2 ] (Gorgues et al., 2010). In the northern high latitudes, an increase in δT reinforces the amplification. In general, the δT T contribution gains importance as we move poleward in both hemispheres and therefore the second-order terms originating from pCO 2 · δT also reinforce the amplification. Interestingly, in the high latitudes, the amplification through second-order terms is as important as the change in the seasonality of the drivers.
The Southern Ocean is an exception to the T dominance; in this region the pCO 2 effect on the DIC s contribution dominates, and the regional amplification is reinforced by low values of the mean buffering capacity (high γ DIC s ). This result agrees with the findings of Hauck and Völker (2015). In this area the amplification is counteracted by a reduction in δDIC s .

Conclusions
In this study, we used output from seven CMIP5 global models, subjected to the RCP8.5 radiative forcing scenarios, to provide a comprehensive analysis of the characteristics and drivers of the intensification of the seasonal cycle of pCO 2 between present (2006-2026) and future (2080-2100) conditions. By 2080-2100 the δpCO 2 will be 1.5-3 times larger compared to 2006-2026. The projected amplification by the Earth system models and the possible causes of it are consistent with observation-based amplification for the period from 1982 to 2015 (Landschützer et al., 2018). However, the models slightly overestimate the present-day amplification, probably due to the larger pCO 2 trends in models than observations (Tjiputra et al., 2014).
The models confirm the well-established mechanisms controlling present-day δpCO 2 (Takahashi et al., 2002;Sarmiento and Gruber, 2006;Fay and McKinley, 2017). DIC s and T contributions are the main counteracting terms dominating the seasonal evolution of δpCO 2 . Furthermore, the models show that under future conditions the controlling mechanisms remain unchanged. This result confirms the findings of Landschützer et al. (2018) that identified the same regional controlling mechanism for the past 30 years. The relative role of the DIC and T terms is regionally dependent. High latitudes and upwelling regions, such as the California current system and the coast of Chile, are dominated by DIC s and the temperate low latitudes are driven by T . Only in the North Atlantic and northwestern Pacific do the models show a dominance of thermal effects over nonthermal effects, which is in disagreement with observations. This further illustrates the urgent need for models to accurately represent regional oceanographic features to accurately reproduce the δpCO 2 characteristics.
In agreement with Landschützer et al. (2018), the model projections towards the end of this century also demonstrate that the global amplification of δpCO 2 is due to the overall long-term increase in anthropogenic CO 2 . A higher oceanic background CO 2 concentration enhances the effect of Tdriven solubility changes on δpCO 2 and alters the seawater carbonate chemistry, also enhancing the DIC seasonality effect. The spatial differences of δpCO 2 amplification, however, are determined by the regional sensitivities and seasonality of pCO 2 drivers. For example, polar regions show larger sensitivity to DIC and T and larger seasonal cycles of DIC and T . Therefore, these areas present a strong enhancement of δpCO 2 in spite of smaller changes in mean pCO 2 .
Moreover, the pCO 2 seasonal cycle amplitude depends on the relative magnitude and phase of the contributions. The models ensemble mean reproduces the highly effective compensation of DIC s and T contributions when they are 6 months out of phase, confirming previous studies (Takahashi et al., 2002;Landschützer et al., 2018). The compensation of DIC and T prevents a larger amplification of δpCO 2 , even when both contributions are largely amplified.
The amplification of the TA and S contributions has a small impact on δpCO 2 in most regions, except in the high latitudes at which the TA contribution complements the DIC one, enhancing the nonthermal effect in this region.
The use of Earth system models allowed us to state the importance of including future changes in driver seasonalities for future δpCO 2 projections. The T seasonality is projected to increase in most of the ocean basins, thereby reinforcing the δpCO 2 amplification. The δT increase is consistent with an increase in stratification that will confine the seasonal changes in net heat fluxes to a shallower mixed layer (Alexander et al., 2018). The DIC s seasonality decreases in some cold areas and its reduction prevents a larger amplification. For the sensitivities, while γ DIC increases, γ T decreases. The latter phenomenon needs further study.
The increasing amplitude of δpCO 2 might have implications for the net air-sea flux of CO 2 , in particular in regions where there is an imbalance between winter and summer values (Gorgues et al., 2010). Examples of such behavior can be found in the Southern Ocean (between 50 and 60 • S) (Takahashi et al., 2014a) and in the latitude band from 2-40 • in both hemispheres (Landschützer et al., 2014). Moreover, seasonal events of high pCO 2 could have an impact on acidification, aragonite undersaturation events (Sasse et al., 2015), and hypercapnia conditions (McNeil and Sasse, 2016). Therefore, understanding the drivers of future δpCO 2 may help to better assess the response of marine ecosystems to future changes in carbonate chemistry. Finally, our complete analytical expansion of δpCO 2 in terms of its four variables provides a practical tool to accurately and quickly diagnose temperature and salinity sensitivities from observational or modeling data sets.