Seasonal Variability of the Surface Ocean Carbon Cycle: A Synthesis

The seasonal cycle is the dominant mode of variability in the air-sea CO 2 flux in most regions of the global ocean, yet discrepancies between different seasonality estimates are rather large. As part of the Regional Carbon Cycle Assessment and Processes Phase 2 project (RECCAP2), we synthesize surface ocean p CO 2 and air-sea CO 2 flux seasonality from models and observation-based estimates, focusing on both a present-day climatology and decadal changes between the 1980s and 2010s. Four main findings emerge: First, global ocean biogeochemistry models (GOBMs) and observation-based estimates ( p CO 2 products) of surface p CO 2 seasonality disagree in amplitude and phase, primarily due to discrepancies in the seasonal variability in surface DIC. Second, the seasonal cycle in p CO 2 has increased in amplitude over the last three decades in both p CO 2 products and GOBMs. Third, decadal increases in p CO 2 seasonal cycle amplitudes in subtropical biomes for both p CO 2 products and GOBMs are driven by increasing DIC concentrations stemming from the uptake of anthropogenic CO 2 (C ant ). In subpolar and Southern Ocean biomes, however, the seasonality change for GOBMs is dominated by C ant invasion, whereas for p CO 2 products an indeterminate combination of C ant invasion and climate change modulates the changes. Fourth, biome-aggregated decadal changes in the amplitude of p CO 2 seasonal variability are largely detectable against both mapping uncertainty (reducible) and natural variability uncertainty (irreducible), but not at the gridpoint scale over much of the northern subpolar oceans and over the Southern Ocean, underscoring the importance of sustained high-quality seasonally resolved measurements over these regions.

• pCO 2 seasonal cycle amplitude changes over 1985-2018 are detectable against both mapping uncertainty and natural variability uncertainty • The dominant driver of pCO 2 amplitude increases over decadal timescales is attributed to the direct effect of C ant invasion • A discrepancy is found with surface dissolved inorganic carbon (DIC) seasonality being systematically less in global ocean biogeochemistry models than in surface DIC observation-based products

Supporting Information:
Supporting Information may be found in the online version of this article.

Introduction
How is the ocean carbon cycle changing as a consequence of sustained increases in emissions of carbon to the atmosphere?Important steps toward answering this question over the last several decades have been provided via estimates of ocean carbon uptake from both interior hydrographic measurements (Gruber et al., 2019;Müller et al., 2023;Sabine et al., 2004), surface fluxes inferred from measurements of the sea surface partial pressure of CO 2 (pCO 2 ) (Chau et al., 2022(Chau et al., , 2023;;Gregor et al., 2019;Landschützer et al., 2014;Rödenbeck et al., 2015), global ocean biogeochemistry model simulations (Friedlingstein et al., 2022a(Friedlingstein et al., , 2022b;;Hauck et al., 2020;Orr et al., 2001) and ocean inverse models (Gloor et al., 2003;Gruber et al., 2009).A first global synthesis was performed roughly a decade ago through the REgional Carbon Cycle Assessment and Processes (RECCAP) project (https://www.globalcarbonproject.org/reccap), highlighting inevitable forced carbon cycle changes, while also identifying sources of uncertainty.
In a parallel direction of inquiry, it has also become clear that the seasonal cycle in surface ocean pCO 2 (and thereby air-sea CO 2 flux) has been changing, with this first identified in modeling studies (Gallego et al., 2018;Gorgues et al., 2010;Hauck & Völker, 2015;Joos et al., 2023;Lerner et al., 2021;Riebesell et al., 2009;Rodgers et al., 2008;Rustogi et al., 2023) and subsequently inferred from observational constraints (Fassbender et al., 2018;Landschützer et al., 2018) and more recently also identified from changes in the seasonal cycle of atmospheric CO 2 induced by air-sea CO 2 fluxes in the Southern Ocean (Yun et al., 2022).Studies to date show that the dominant drivers of this trend toward increasing pCO 2 seasonal amplitude are the increasing surface concentrations of dissolved inorganic carbon (DIC), due to the invasion of anthropogenic CO 2 into the ocean, and the associated decline in surface ocean buffering capacity.The trend toward increased seasonal amplitude is also expected to be modulated by natural variability and by warming (Schlunegger et al., 2019).The net effect of warming in high emission projections by the end of the 21st century is to enhance increases in the seasonal amplitude of CO 2 fluxes over the expansive subtropical domain by approximately 15% (Rodgers et al., 2020).
Our synthesis effort builds on the earlier studies of Fassbender et al. (2018), Landschützer et al. (2018), and Fassbender et al. (2022), but here we rely on a greatly expanded set of pCO 2 products and GOBMs, spanning the time period 1985 to 2018, to assess the drivers of multi-decadal changes in the seasonal cycle of surface pCO 2 .Such changes are of interest for three primary reasons.First, they may impact stressors of marine ecosystems, as emphasized in the analysis of future changes in pH, hydrogen ion concentration, and calcium carbonate mineral saturation state seasonality by Kwiatkowski and Orr (2018).Second, they can have potential impacts for climate feedback mechanisms, as demonstrated by Fassbender et al. (2022) for 21st century climate change under both strong and moderate CO 2 emissions pathways.Third, they will have implications for optimizing the design of the global marine carbon observing system, as summer-biased measurements can lead to significant errors for a system where the seasonal cycle itself is not stationary.
For the majority of the global ocean, the seasonal cycle in sea surface pCO 2 exhibits strong variability that is up to two orders of magnitude larger than the seasonal cycle in atmospheric pCO 2 (e.g., Takahashi et al., 2002).Seasonal variations in pCO 2 reflect the interplay of four principal underlying drivers: sea surface temperature (SST), sea surface salinity (SSS), dissolved inorganic carbon concentrations (DIC), and total alkalinity (TA) (Takahashi et al., 1993(Takahashi et al., , 2002)).Variations in these drivers are the result of changes in ocean circulation and mixing, biological processes, as well as exchanges of heat, freshwater, and carbon with the atmosphere.While the temperature-driven or thermal component of observed pCO 2 seasonality is well constrained by laboratory studies and by high-quality satellite-derived SST products, this is not the case for the other drivers.Therefore, observational studies often refer to changes related to DIC, TA, or salinity collectively as changes in the nonthermal component (Takahashi et al., 1993(Takahashi et al., , 2002)).Importantly, thermal and nonthermal drivers of pCO 2 seasonality are typically found to be opposed in phase and large in amplitude relative to the resulting net pCO 2 variations.For this reason, skillful modeling of seasonal pCO 2 cycling can be compromised if any of the underlying processes are not well-represented, both for the amplitude and the phase/timing of seasonality changes.Therefore, state-of-the-art models still face challenges in skillfully simulating the seasonal cycle of pCO 2 in regions that are known to be sustained, seasonally unbiased measurements over the high latitudes as part of an optimized marine carbon observing system. 10.1029/2023GB007798 3 of 34 important for CO 2 exchange with the atmosphere (Anav et al., 2013;Goris et al., 2018;Kessler and Tjiputra, 2016;Mongwe et al., 2018;Pilcher et al., 2015).For the important case of the North Atlantic, it has been shown that future carbon uptake in ESM projections can be constrained by diagnosing the timing and magnitude of modeled seasonal variations in pCO 2 (Goris et al., 2018).Similarly, Nevison et al. (2016) demonstrated, for the Southern Ocean, that the way an ESM represents the seasonal carbon cycle is related to the model's projected carbon uptake.A much weaker seasonal cycle than in higher latitudes is observed in the tropics, where seasonal pCO 2 amplitudes are obscured by interannual variability.In polar latitudes where there is seasonal ice cover, the challenge of data sparsity associated with the difficulty of sampling pCO 2 through the seasonal cycle presents challenges and uncertainty that are beyond the scope of this study.
An important source of uncertainty with detecting anthropogenic trends in the carbon cycle is natural variability.As a matter of definition, we consider the forced anthropogenic signal to represent the total sum of the effects of anthropogenic forcing on atmospheric radiation (greenhouse gases, anthropogenic aerosols, etc.) and thereby warming and the ensuing physical and biogeochemical state changes in the ocean, as well as the effect of the invasion flux of anthropogenic carbon into the ocean.For the analysis presented here, even if pCO 2 products were able to perfectly resolve multi-decadal changes in the seasonal amplitude of sea surface pCO 2 over 1985 to 2018, it would not be possible to confidently isolate the anthropogenic trend from the impacts of natural variability.Given that natural variability is an intrinsic property of the Earth system, natural variability uncertainty is understood to be an irreducible form of uncertainty that can nevertheless be quantified.To address the degree to which natural variability obscures detection of anthropogenic trends, we include analyses of large ensemble (LE) simulations with Earth system models (ESMs) in this study.LE simulations were considered previously in the study of Schlunegger et al. (2019) for the case of 21st century climate change projections under strong emissions, and they showed that over the subtropics forced changes in ΔpCO 2 seasonal amplitude are more emergent than forced trends in annual mean ΔpCO 2 .
Previous research has pointed toward significant model biases in representing the seasonal cycle of pCO 2 and its underlying drivers (Goris et al., 2018;Mongwe et al., 2018;Rodgers et al., 2014).The question naturally arises as to whether models have improved since the evaluation of seasonal variability in an early generation of prognostic biogeochemistry models (McKinley et al., 2006;Tjiputra et al., 2012).For the case of the North Pacific, the forward models in McKinley et al. (2006) had their seasonal cycles in pCO 2 evaluated directly against time series measurements from observing stations spanning the subpolar and subtropical gyres to identify non-trivial biases in models.
Our study differs in that, rather than direct comparisons with time series stations, our observation-based data sources consist of mapped global sea surface pCO 2 products constructed from the same underlying ungridded SOCAT (Surface Ocean CO 2 Atlas; Pfeil et al., 2013) data, as well as ocean interior seasonal climatological fields of DIC.
More extensive analyses for the regions with seasonal ice cover, which are not considered here due to their sparsity of carbon observations for estimating the seasonal cycle, are provided by the Arctic Ocean (Yasunaka et al., 2023) and Southern Ocean (Hauck et al., 2023) regional contributions to RECCAP2.Our analysis and synthesis is developed using a diverse set of resources: nine observationally based pCO 2 products, three of which include associated time-varying monthly surface DIC and TA fields; two independent three-dimensional seasonally resolved DIC climatologies; hindcast simulations from eleven GOBMs; five large ensemble simulations; and individual realizations of 11 CMIP6 models.
This paper is structured as follows.We begin in the methods section by introducing the RECCAP2 pCO 2 products and GOBMs, along with ancillary products and Earth system model simulations used for this study.In the main results section, we present a descriptive account of changes in pCO 2 seasonality as well as discrepancies between pCO 2 products and GOBMs, followed by attribution of these changes and discrepancies to their drivers.We finally present an analysis focused on identifying detectability of forced/anthropogenic signals, as well as an assessment of the robustness of our chosen biome aggregation for detectability of the forced changes in seasonality to changes due to natural variability.Our objective is to present a synthesis of the carbon cycle community's current knowledge of how the seasonal cycle of sea surface pCO 2 and the air-sea CO 2 exchange has been changing over the last three decades.

Considered Regions
Ocean biomes offer a number of advantages as aggregation regions for studying the global carbon cycle (Fay & McKinley, 2014).For our interests in seasonality, biomes appropriately reflect the structures that are determined by real-world interactions between ocean circulation and biogeochemistry.For all of the oceanic studies within RECCAP2, a discrete number of ocean biomes based on Fay and McKinley (2014) are used to facilitate consistent intercomparison between regions (described in the supplement to the Müller (2023) publication of the RECCAP2 data).This study analyzes pCO 2 and CO 2 flux seasonality for six biomes (Table 1 and central panel of Figure 1) aggregated from the aforementioned.Our decision to aggregate some biomes and exclude others is based on the following: First, this is a synthesis study and the seasonal cycle in surface ocean pCO 2 exhibits important similarities across the subtropics of each hemisphere.Second, we have elected to not include either the Arctic, the Bering Sea, or the regions of the Southern Ocean that are impacted by ice cover, with the analysis of those regions  left to the individual regional RECCAP2 studies, as mentioned in the Introduction.Third, we have chosen to not include the equatorial regions, given that the seasonal cycle there is relatively weak and often obscured by large interannual variability.
Thus, our six aggregated biomes (Table 1 and Figure 1) (their precise boundaries given in Supporting Information S1 of the RECCAP2 data release of Müller, 2023) consist of North Atlantic subpolar seasonally stratified (NA-SPSS), North Pacific subpolar seasonally stratified (NP-SPSS), Northern Hemisphere subtropical seasonally stratified (NH-STSS), Northern Hemisphere subtropical permanently stratified (NH-STPS), Southern Hemisphere subtropical permanently stratified (SH-STPS), and Southern Hemisphere seasonally stratified (subpolar and subtropical combined, SH-SS).We have elected to keep the North Atlantic and North Pacific subpolar biomes separate in our analysis in light of important differences that arise in their seasonal cycles.
Besides the aggregation as designated in Table 1, adjustments have been taken to avoid confounding changes in observational coverage with seasonality or changes in seasonality.We limit our analysis to regions that have valid observation-based (and GOBM) fields for the whole time period 1985-2018 for all pCO 2 products and for all variables considered in this study (CO 2 fluxes, pCO 2 , SST, SSS, surface DIC, and surface TA).Thereby we make sure that the same region is considered for all variables over the whole seasonal cycle but also over the whole observational period (see Figure 1).

Temporal Filtering of RECCAP2 Products
In general, the monthly mean pCO 2 products and GOBMs considered here can be expected to represent pCO 2 seasonality modulations that are sustained through the combined effects of the invasion flux of C ant , natural variability, and warming.In order to identify the amplitude of the seasonal cycle (with its modulations), we follow Landschützer et al. (2018) in using a quadratic polynomial fit to remove the decadal trend from monthly time series over 1985-2018.Our analysis through the majority of this manuscript characterizes multi-decadal changes in pCO 2 seasonal variations by considering differences between the 5-year intervals 1985-1989 and 2014-2018.We also investigated the degree to which the main results are sensitive to the choice of 5-year versus 10-year versus 15-year intervals for describing climatologies (Figure S3 in Supporting Information S1), and upon finding that the results are relatively insensitive to the approach we opted to use the difference between the two 5-year intervals 1985-1989 and 2014-2018.As such, our approach is consistent with that used by Landschützer et al. (2018), who also used 5-year climatologies in calculating decadal changes.

Surface Observation-Based Products of pCO 2 , DIC, and TA
Observation-based pCO 2 products that are compared for their representation of pCO 2 seasonality are listed in  database (SOCAT, Bakker et al., 2016) yet they differ significantly in the way they spatiotemporally interpolate or map the sparse observations.Methodologically they include multiple linear regression (Iida et al., 2021), gradient boosted decision trees (Gloege et al., 2022), neural networks (Chau et al., 2022;Landschützer et al., 2016;Zeng et al., 2022), ensembles of various machine learning approaches (Gregor et al., 2019;Gregor & Gruber, 2021), and a Bayesian approach that is additionally constrained by mixed-layer dynamics (Rödenbeck et al., 2013(Rödenbeck et al., , 2022)).Additionally, one neural network product applies adjustments to the SOCAT data to account for temperature gradients between the ocean surface and the depth of the seawater intake on the observing ships and for the cold skin layer effect (Watson et al., 2020).
Of central importance to our study is that three of these products (JMAMLR, OceanSODA-ETHZ, and CMEMS-LSCE-FFNN) additionally provide time-varying surface DIC and TA products spanning 1985-2018.
While we use all available pCO 2 products listed in Table 2 for our analysis of seasonal cycles of CO 2 -fluxes and pCO 2 , our analysis greatly benefits from using these three products to attribute changes of pCO 2 seasonality to its drivers.We will refer to them as the pCO 2 /TA products.We provide a summary of how DIC and TA are derived for these three pCO 2 products in Supporting Information S1, but we refer the reader to the papers cited in Table 2 for a more complete description of methods for the individual observation-based products.

Three-Dimensional DIC Climatologies
For our evaluation of seasonal variability in surface DIC concentrations, we consider two three-dimensional DIC climatologies constructed from observational products: the Mapped Observation-Based Oceanic DIC (MOBO-DIC) product by Keppler et al. (2020aKeppler et al. ( , 2020b)), and the NNGv2LDEO monthly climatology of interior DIC by Broullón et al. (2020).Unlike the pCO 2 /TA products, these products are based on direct observations of DIC and as such provide an independent estimate of surface DIC seasonality, with the caveat that they offer a climatological view only.Both of these three-dimensional DIC climatologies are described in greater detail in Supporting Information S1.

Global Ocean Biogeochemistry Models (GOBMs)
Global ocean biogeochemistry models (GOBMs) are compared against the pCO 2 products for modulations of CO 2 fluxes and pCO 2 for the six aggregated biomes over the period 1985-2018.Many of the GOBMs considered here have contributed to the Global Carbon Project's estimate of the ocean carbon sink (Friedlingstein et al., 2022a(Friedlingstein et al., , 2022b)), and have been evaluated by Hauck et al. (2020).Both the phase and amplitude of seasonal variability in the GOBMs will be evaluated against the pCO 2 products, as will the behavior of the drivers (SST, DIC, SSS, and TA).The GOBMs are listed in Table 3, and their forcing fields are heterogenous, including JRA-55 (Kobayashi et al., 2015), CORE (Large & Yeager, 2009), and forcing derived from the NCEP/NCAR reanalysis (Kalnay et al., 1996;Kistler et al., 2001).All but one of the GOBMs listed in Table 3 provide the full suite of variables and timescales needed for this analysis.The exception is the CCSM model for which output is available only through 2017 (rather than 2018).For this model alone we construct climatologies over 2014-2017.
Additionally we have included the abiotic data-assimilation model OCIM (DeVries, 2014(DeVries, , 2022) ) where appropriate in our analysis.Although more commonly applied to the uptake of anthropogenic carbon, the abiotic nature of the OCIM model is of value as an endmember case in evaluating potential biotic biases in the GOBMs.
Air-sea heat fluxes for the GOBMs used throughout RECCAP2 are calculated with bulk formulae following the protocols of Large and Yeager (2009).This method does not impose a specific nudging or restoring of simulated sea surface temperature (SST) to observed SST, but rather imposes a negative feedback to SST through heat fluxes determined using observed surface boundary layer atmospheric temperatures.In this way, the simulation of seasonal variations in SST shows strong fidelity to observations when aggregated over biome scales.

Attribution of Drivers of the Climatological Seasonal Cycle in pCO 2
The first stage of the attribution analysis has the goal of identifying the mechanisms responsible for discrepancies between pCO 2 products and GOBMs in representing the climatological seasonal cycle of pCO 2 .We will consider both a decomposition into thermal and nonthermal drivers, using the analysis methods presented by Fassbender et al. (2022).The thermal/nonthermal decomposition provides a valuable means to identify spatially regions where seasonal SST variations are or are not the dominant driver of seasonal variations in pCO 2 .We will also apply the Taylor Series decomposition method previously considered by Sarmiento and Gruber (2006), where we use the formulation: For this formulation, we consider salinity-normalized DIC and TA (nDIC and nTA, respectively), to remove the dilution effects from the DIC and TA contributions.For our interests here, the terms on the right-hand-side of Equation 1 will provide estimates of the sensitivity of pCO 2 seasonality, and thus its amplitude and phasing, to the seasonal evolution of nDIC, nTA, SST, and SSS.In this way we will evaluate any discrepancies that may call into question the fidelity of modeled pCO 2 to the real ocean.This analysis is possible for all GOBMs and the three observation-based products that include DIC and TA.The details of the coefficients used in Equation 1 are provided in Equations S1-S9 of the Supporting Information S1, with averaged values for the pCO 2 /TA products provided in Table S1 in Supporting Information S1.
Additionally mixed layer depth (MLD) products will be used as part of our attribution analysis.For the sake of consistency with our choice of the relatively data-rich period 2014-2018 for evaluating climatologies of pertinent variables, we present a MLD product for this period 2014-2018 derived from the gridded monthly Argo product of Roemmich and Gilson (2009) for temperature and salinity.We define MLD by a density threshold through application of the MLD definition of Holte and Talley (2009).

Attribution of Decadal Changes in the Seasonal Cycle of pCO 2
The second stage of attribution analysis consists of applying the method of Fassbender et al. (2022) to identify the drivers of decadal changes in pCO 2 seasonality over 1985-2018.The method provides a means to separate the impacts of the invasion flux of anthropogenic carbon from the atmosphere from other climate signals that can modulate the seasonal cycle in pCO 2 .More specifically, the goal for this work will be to identify the direct impact of an increase in the anthropogenic carbon (C ant ) content of surface waters (i.e., a decrease of the surface ocean buffering capacity) relative to climate change and natural variability impacts on pCO 2 seasonality changes over 1985-2018.The original application by Fassbender et al. (2022) was for 21st century LE projections with an ESM, where the ensemble mean of the LE was considered to identify the forced component of change.Here our intention is to apply the method to the period 1985-2018 for all GOBMs (Table 3) and the three gridded pCO 2 /TA products that include associated monthly DIC and TA fields along with SST and SSS (Table 2).

Note.
Model names are given in the first column, references are given in the second column, and comments pertinent to the application of the products are given in the third column.The CCSM model differs from the others in that output only extends to 2017 rather than to the RECCAP2 protocol end-year of 2018.Table 3 Global Ocean Biogeochemical Models (GOBMs) and an Assimilation Model Used for This Study The method of Fassbender et al. ( 2022) is a modified version of the framework originally developed by Takahashi et al. (1993Takahashi et al. ( , 2002)), where for each grid point local temperature and pCO 2 are used as follows to define the thermal (pCO 2 T ) and nonthermal (pCO 2 NT ) components of pCO 2 seasonal variability: (2) (3) Within this Takahashi et al. (1993) approach of Equations 2 and 3, the subscripts am and mm represent annual means and monthly means, respectively.However, one of the limitations of the method of Takahashi et al. (1993Takahashi et al. ( , 2002) ) is that the two components in Equations 2 and 3 do not sum to reproduce the full pCO 2 seasonal cycle (e.g., pCO 2 mm ≠ pCO 2 T + pCO 2 NT − pCO 2 am ).This is because the thermal sensitivity of pCO 2 (expressed with a factor of 0.0423/°K in the equations above) in fact varies slightly with background chemistry (Wanninkhof et al., 1999(Wanninkhof et al., , 2022)).
In this study, we use the alternative approach of Fassbender et al. (2022), where an important focus of the analysis is to identify asymmetries in seasonal pCO 2 variations.For this case, rather than working with an annual mean pCO 2 to characterize asymmetries, one defines a neutral pCO 2 with respect to the annual cycle as being the pCO 2 one would obtain using annual means of its drivers (T, S, DIC, TA, PO 4 , and SiO 4 , where we use T and S as shorthand for SST and SSS), with this reflecting the mean pCO 2 if pCO 2 were to respond linearly to the seasonal variability of its drivers.With this in mind, pCO 2 AM is given by: where the function f() represents a carbon system calculator.The time-varying (monthly mean) thermally driven pCO 2 component can be identified using the same approach but by using monthly varying output temperature for the carbonate system calculator: where the subscript FASS and the term on the left-hand side denote that the calculation follows the method developed by Fassbender et al. (2022), as this is the approach that will be applied in this study.Within this approach, the thermal pCO 2 component seasonal cycle anomaly is then determined as: making it possible to calculate the nonthermal component by difference: The subscript FASS here as well denotes that this is calculated following the method described in Fassbender et al. (2022), rather than Takahashi et al. (1993Takahashi et al. ( , 2002)).These definitions for thermal and nonthermal decomposition in Equations 6 and 7, are applied in Figure 6.To attribute decadal changes in pCO 2 seasonality to C ant and identify whether decadal changes project differently onto summer and winter conditions, we independently reconstruct the nonthermal pCO 2 term following the method outlined in Supporting Information S1.

Uncertainty Analysis for Detecting the Anthropogenic Signal
As in the study of Fassbender et al. (2022), we apply large ensemble (LE) simulations with ESMs to aid in our identification of forced changes in pCO 2 seasonality in the presence of non-negligible background variability.Whereas Fassbender et al. (2022) focused on centennial timescale projections and could rely on the ensemble mean to isolate forced signals, we are here faced with a different challenge of identifying forced changes over shorter historical records (namely the RECCAP2 period spanned by 1985-2018).Stated differently, for the decadal changes in pCO 2 seasonality that we intend to identify for the nine pCO 2 products in this study, we wish to apply LE simulations to estimate quantitatively our degree of confidence that the observed decadal changes represent anthropogenically forced signals.
The LE models considered here are listed in Table 4.The first set of models were run under CMIP5 protocols using historical/RCP8.5 forcing (CanESM2, ESM2M, and CESM1) and the second set were run under CMIP6 protocols with historical/SSP3-7.0 forcing (CanESM5 and CESM2).Such a mix of CMIP5 and CMIP6 is appro-priate to our interest in the historical period spanning 1985-2018, as the total radiative forcing component for the CMIP5 and CMIP6 forcing pathways is very similar through 2018 (Riahi et al., 2017).These LE outputs are not part of the resources made available through RECCAP2 but were rather collected by the authors through a modest expansion of the collection of LE models considered by Schlunegger et al. (2020) and Gloege et al. (2021).We are particularly interested in evaluating confidence levels for emergent changes between the time intervals 1985-1989 and 2014-2018, where the ensemble mean changes (signal) and natural variability in the changes (noise) provide a means to identify the signal-to-noise ratio (SNR) and thereby characterize the degree of emergence or detectability.
For the attribution of decadal changes in pCO 2 seasonal amplitude, we have applied the CESM2-LE, as this was deemed to have better correspondence with climatological pCO 2 variability over biome scales than the other models (this will be addressed in Figure 12).The application of CESM2-LE will allow us to address the degree of confidence that we have in distinguishing the impact of the invasion flux of anthropogenic carbon from climate-driven perturbations.CESM2-LE is also applied to identify the degree to which thermal and nonthermal drivers of pCO 2 seasonality changes are emergent.

CMIP6 Models
This study also includes analyses of 11 CMIP6 ESMs that are evaluated for their agreement with the pCO 2 products over biome scales, with the models listed in Table 6.As ESMs are used for future projections of the ocean carbon sink, it is of interest to assess their fidelity in reproducing the seasonal cycle and its underlying processes.The CMIP6 models are not intended to correspond to coupled versions of the GOBMs considered here, but our goal is rather to identify similarities or discrepancies between these CMIP6 models, the GOBMs, and observational pCO 2 products.For the CMIP6 models we have opted to focus on a slightly different period, namely the five years (2010-2014) at the end of the historical component of CMIP6 simulations, rather than the 2014-2018 interval considered for the GOBMs.

Overview of Multi-Decadal Changes in CO 2 Flux and pCO 2 Seasonality
The evolution of air-sea CO 2 flux seasonality, considered as averages over the six biomes, is shown in Figure 1 as time series of the annual minimum and annual maximum of the monthly sea-air CO 2 fluxes for the pCO 2 products and GOBMs.Here, we have chosen to show the maximum and minimum of the fluxes to represent the full range of the seasonal cycle.One should keep in mind that the across-product ensemble spread also includes the impacts of wind speed on gas exchange that are not treated the same way across products (the details for the individual products can be found in the references in Table 2).
A prominent feature in Figure 1 is that the amplitude of the seasonal cycle in CO 2 fluxes is consistently larger for the GOBMs than it is for the pCO 2 products, across all biomes (see also Figure S1 in Supporting Information S1).From the RECCAP2 GOBMs alone it is not possible to determine the degree to which these internal disagreements reflect structural differences between the GOBMs and the degree to which it represents differences in forcing of the GOBMs.It is also worth noting that both the GOBMs and pCO 2 products exhibit a non-negligible degree of natural variability (interannual to decadal), for both the maxima and minima.
A decadal trend toward increased seasonal amplitude of biome-integrated CO 2 fluxes is evident over the subtropical biomes, as well as for NP-SPSS, whereas for NA-SPSS and SH-SS biomes this is less evident (especially in the pCO 2 products, see also Figure S1 in Supporting Information S1).The increase in the amplitude, where it exists, can be accounted for mainly by a decrease of the annual minima (i.e., an increase of the uptake flux), whereas the annual flux maxima remained comparably constant over time.To better constrain whether there are discernible decreases in the spread amongst the pCO 2 products over 1985-2018 as pCO 2 observations increase, the biome-averaged pCO 2 for seasonal maximum and seasonal minimum are shown in Figure S2 in Supporting Information S1, where although there is a decrease for both NA-SSS, SH-STPS, and SH-SS between 1985 and 2018, for the other biomes the changes are rather small.The broad-scale view of the evolution of CO 2 fluxes in Large Ensemble Simulations Used in the Analysis in Figure 12 Figure 1 serves to motivate much of the analysis that follows, with emphasis devoted to pCO 2 rather than CO 2 fluxes as a means to facilitate attribution and mechanistic understanding.Nevertheless, we will return to CO 2 fluxes in the discussion of decadal trends in seasonal cycle amplitude.

Climatological Seasonal Cycle of pCO 2 and Multi Decadal Changes
The climatological monthly pCO 2 seasonal cycle over 2014-2018 as a zonal mean is shown in Figure 2.For the pCO 2 products (Figure 2a), pCO 2 is highest in the subtropics in summer and highest in the Northern Hemisphere subpolar region and in the Southern Ocean in winter.The GOBMs (Figure 2e) agree with the pCO 2 products in their seasonal phasing in the subtropics but are in disagreement in the Northern Hemisphere subpolar regions and the Southern Ocean.As a measure of the internal consistency across the nine pCO 2 products and the 11 GOBMs for this 2014-2018 climatology, the standard deviation is calculated month-by-month across the pCO 2 products (Figure 2b) and across the GOBMs (Figure 2f).From this we learn that the agreement across the pCO 2 products is robust (<10% disagreement relative to the annual cycle amplitude), whereas internal disagreements are much larger across the GOBMs (larger than 20% of the annual cycle amplitude at some locations), with the GOBMs diverging the most north of 30°N and south of 30°S.For the Northern Hemisphere subpolar regions, the disagreement in the GOBMs is largest in summer, and over the Southern Ocean it is largest in summer and in late winter.We also consider the multi-decadal changes in the climatological monthly pCO 2 seasonal cycle between 1985-1989 and 2014-2018 for the pCO 2 products (Figure 2c) and for the GOBMs (Figure 2g).For both the pCO 2 products and the GOBMs (Figures 2c and 2g, respectively), the meridional patterns of increases in seasonal amplitude to first order reflect an amplification of the climatological seasonal cycle in the Northern Hemisphere subtropical and subpolar biomes.However, this is less clear in the Southern Hemisphere subtropics and over the Southern Ocean.

Spatial Structure in the Seasonal Amplitude Change
To understand the spatial structure of the amplitude change and the coherence between pCO 2 products and between GOBMs, maps of changes in the amplitude of the pCO 2 seasonal cycle between 1985-1989 and 2014-2018 are shown for the pCO 2 products in Figure 3a and for the GOBMs in Figure 3b.We use two definitions for the amplitude assessment: winter-minus-summer (JFM minus JAS for the Northern Hemisphere and JAS minus JFM for the Southern Hemisphere) in Figures 3a and 3b and climatological monthly minimum-minus-maximum in Figures 3c and 3d.For both the GOBMs and the pCO 2 products, regions where the multi-decadal trends are not significant (defined using a one standard deviation threshold) are stippled.From this we can see that the Northern Hemisphere has more significant changes than the Southern Hemisphere.Nevertheless, for both hemispheres for GOBMs and pCO 2 products the broadest regions of significant changes are in the subtropics, although there is also a significant trend in the eastern subpolar North Pacific.
The two different definitions of seasonal amplitude (winter-minus-summer and minimum-minus-maximum) result in very similar large-scale patterns.While for the minimum-minus-maximum definition negative changes always indicate an increase in seasonal amplitude, the winter-minus-summer amplitude is negative over the subtropical domain and positive in high latitudes, such that negative (positive) values for the multi-decadal amplitude change in Figures 3a and 3b over the subtropics (high latitudes) also indicate an amplitude increase.
Phasing differences in seasonality are treated differently with the two definitions, since winter-minus-summer refers to fixed time intervals while the minimum and maximum can occur at any time throughout the year.While For this analysis, the nine pCO 2 products listed in Table 2 and the 11 GOBMs in Table 3 (excluding OCIM) are used.
The sensitivity to the choice of 5-year versus 10-year versus 15-year averaging for characterizing decadal changes is displayed in Figure S3 in Supporting Information S1.
the winter-minus-summer definition highlights differences in phasing between pCO 2 products and GOBMs, it conceals some of the changes in the seasonal amplitude, specifically if the minimum or maximum of the seasonal cycle occurs outside of winter or summer months (as seen in Figure 3d in the high latitude North Atlantic and North Pacific, where the discrepancy between GOBMs and pCO 2 products is much larger for the winter-minus-summer definition).
For the remainder of our analysis we will use the winter-minus-summer definition of seasonality due to its simplicity when it comes to analysis of multiple models and products, and also for consistency with previous work (e.g., Landschützer et al., 2018), keeping in mind that the specifics of this definition (the choice of JFM and JAS for defining seasonal amplitude) are sensitive to phasing differences between GOBMs and pCO 2 products.The sensitivity of the choice of 5-year versus 10-year versus 15-year climatologies for the case of the pCO 2 products (Figure S3 in Supporting Information S1) indicates that the results are robust to the length of the climatologies, the main difference being that using 10-year or 15-year climatologies shortens the effective decadal-change interval and weakens the signal of interest.
The fields shown in Figures 3a and 3c reveal important information about uncertainty in detecting multi-decadal changes in pCO 2 seasonality using pCO 2 products.As all of the products are based upon the same un-gridded SOCAT data product, differences in multi-decadal changes shown in Figures 3a and 3c are due to mapping differences rather than differences in the basis data itself.As such, the stippled regions indicate where mapping uncertainty is large relative to resolved decadal changes.We interpret mapping uncertainty to be a reducible variety of uncertainty, in that with improved data coverage and expansion of the ocean observing system the products should be expected to converge.
For the GOBMs (Figures 3b and 3d) the regions of significance are very similar to what is shown for the pCO 2 products, but the underlying reasons for uncertainty are of course fundamentally different.For the GOBMs, regions where trends are not significant reflect the combined effect of structural differences between models and differences in GOBM forcing (as was described in Methods Section 2.5, some of the models were forced with CORE, while others were forced with JRA-55).The uncertainty shown in Figures 3b and 3d for the GOBMs should also be considered as a reducible variety of uncertainty.
Nevertheless, it should also be noted that for both pCO 2 products and GOBMs that it is the subtropics where the thermal component of pCO 2 seasonality has the weakest opposition from the nonthermal component that there is emergence.Thus, both mapping uncertainty and presumably structural uncertainty with GOBMs are mostly problematic in regions where the nonthermal drivers become important.This will be discussed further in subsequent analyses of this study.

Biome Aggregated Climatologies
Climatological monthly anomalies in pCO 2 and its primary drivers, SST and surface DIC, relative to their annual means are displayed in Figure 4 as (area-weighted) averages across the RECCAP2 biome regions.A first comparison between the full collection of pCO 2 products (blue) and all GOBMs (green) is shown in the first column of Figure 4 by biome region, where continuous lines and associated shading represent the mean and one standard deviation, respectively.Our first objective is to identify discrepancies between the GOBMs and the pCO 2 products, with such discrepancies to be addressed subsequently in Section 3.2, where we consider attribution.
In broad terms, the seasonal phasing agreement between GOBMs and pCO 2 products is best in subtropical biomes (NH-STSS, NH-STPS, and SH-STPS) and worst in the Northern Hemisphere subpolar and the Southern Ocean biomes (NA-SPSS, NP-SPSS, and SH-SS).In the NA-SPSS and NP-SPSS biomes, GOBMs consistently represent seasonal pCO 2 minima that occur earlier than in the pCO 2 products, and that are weaker for NA-SPSS.For NA-SPSS and NP-SPSS, after the spring minimum in pCO 2 in the GOBMs there is an increase toward a maximum in August, whereas the pCO 2 products remain low throughout the summer before increasing in autumn and winter toward a maximum in January/February.Analogous to what is found in the Northern Hemisphere subpolar biomes, there are also discrepancies in phasing for the GOBMs relative to the pCO 2 products for SH-SS.For the higher latitude regions where phasing in models diverges greatly from that seen in the pCO 2 products, it is widely understood that the early summer minimum in pCO 2 (as seen in pCO 2 products) is directly related to biological drawdown of DIC (Gregor & Gruber, 2021;Lauderdale et al., 2016).It is a known bias of coupled models that primary production is often not realistically sustained through summer after a pronounced spring bloom.This has been shown in the North Atlantic (Goris et al., 2018;Person et al., 2018), and more broadly for the high latitudes of both hemispheres (Cabré et al., 2016;Mongwe et al., 2018;Nevison et al., 2016), and the same limitation appears to be pertinent in our study. 10.1029/2023GB007798 13 of 34 In the subtropical biomes, on the other hand, the GOBMs consistently simulate a pCO 2 seasonal cycle that is larger than what is seen in the pCO 2 products (Table 5).For NH-STSS the amplitude averaged across the GOBMs is 80% larger than for that averaged across the pCO 2 products, for NH-STPS it is 26% larger, and for SH-STPS it is 34% larger.This occurs with phasing that is largely consistent between the GOBMs and pCO 2 products.As seen in Figure 4 (second column), this discrepancy in the amplitude of pCO 2 variations does not reflect biases in the seasonal amplitude of SST variations, as the thermodynamic boundary conditions for forcing the GOBMs result in SSTs that have relatively strong fidelity to observations when averaged over biome scales.Rather, the discrepancy in the amplitude of the pCO 2 seasonal amplitude indicates that surface DIC (Figure 4, second column) as the main driver of the nonthermal component of pCO 2 seasonality is not providing sufficient counterbalance to the realistic thermal component of seasonal forcing (i.e., the seasonal amplitude seems to be too small in GOBMs).For the subtropical biomes it is worth noting that the OCIM model (Figure 4, thin line) has an even larger pCO 2 seasonal amplitude than the GOBMs.The OCIM model has a realistic SST seasonal evolution (not shown) but simulates only circulation-and surface flux-driven DIC variations without including a contribution of the biological drawdown of DIC.The OCIM pCO 2 and DIC seasonal cycle amplitude discrepancies relative to the GOBMs provide an estimate of the contribution of the biological pump to DIC and pCO 2 seasonality.
The fact that the DIC seasonal amplitude is weaker for the GOBMs relative to the pCO 2 products means that the nonthermal component of pCO 2 seasonality is smaller.Over the subtropical regions, this will serve to enhance the thermal dominance of seasonal pCO 2 variations and in the subpolar and circumpolar regions this will reduce the nonthermal dominance of seasonal pCO 2 variations.This provides a means to alter the phasing of pCO 2 variations, although the degree to which this occurs will be modulated by the background mean state of DIC and TA (Fassbender et al., 2018).In all biomes, the DIC seasonal amplitude for GOBMs lies between that of the abiotic OCIM model and that of the pCO 2 /TA products, suggesting that for the GOBMs there may be as-yet undetermined shortcomings in the representation of the seasonal processes that increase DIC (entrainment) and decrease DIC (net biological consumption of DIC).
Comparing the decadal changes (Figure 4, right column) with the climatologies over 2014-2018 (left column), we see again that the decadal changes in pCO 2 to first order amplify the climatological seasonal cycle in the Northern Hemisphere as well as for SH-STSS, consistent with what was identified in the Hovmöller diagrams in Figure 2. Importantly, the pCO 2 products show significantly less internal disagreement for the climatological seasonal cycle (left column of Figure 4 and Table 5) than the GOBMs.Yet, for the decadal changes (right column of Figure 4), the spread across the GOBMs and across the pCO 2 products is similar (one standard deviation for the biome averaged seasonal amplitude changes are roughly between 1 and 10 ppm, Table 5).While this is of the same order of magnitude as for their seasonal amplitude for the pCO 2 products, GOBMs show significantly less internal disagreement for the decadal change than for the seasonal amplitude.This indicates that the processes responsible for the increase in seasonal pCO 2 amplitude are relatively well represented in GOBMs We hypothesize here that the pertinent process in the GOBMs is the invasion flux of C ant dominating pCO 2 seasonality modulations associated with climate-driven changes in the state of the ocean, a hypothesis to which we will return in the next section on attribution.The best agreement for decadal changes in the pCO 2 seasonal amplitude occurs in the subtropical biomes, consistent with Figure 2.For NH-STSS, the GOBMs show a mean net increase in the seasonal amplitude of 10.7 μatm and the pCO 2 products of 9.8 μatm, with the difference in their increase being less than 10%.For NH-STPS, the mean net increase in the seasonal amplitude for the GOBMs and pCO 2 products are the same at 6.3 μatm.For SH-STPS the increase for GOBMs is 2.7 μatm, which is approximately half of the mean pCO 2 product increase of 5.3 μatm.It can also be seen in Table 5 that the relatively consistent ranges of change over the subtropical biomes between GOBMs and pCO 2 products are approximately within each other's uncertainty range.

Attribution
Having evaluated the seasonal cycle of pCO 2 in terms of its phasing and amplitude, as well as related discrepancies between the pCO 2 products and the GOBMs, we now turn our attention to attribution of the drivers.We have already inferred from the biome-aggregated analyses in Figure 4 that DIC should be explored as having a potential first-order role in modulating the differences between pCO 2 products and GOBMs.

Spatial Patterns of the Seasonal Cycle in Surface DIC Concentrations
The spatial patterns for the surface DIC seasonal cycle amplitude (winter-minus-summer, climatology over 2014-2018) for the pCO 2. /TA products is shown in Figure 5a (average across JMAMLR, CMEMS-LSCE-FFNN,  and OceanSODA-ETHZ, as well as for the GOBMs in Figure 5b.We see a good agreement between the amplitude of the surface DIC seasonal cycle from pCO 2 /TA products and the GOBMs for the overall patterns with maximum values seen in the Northern Hemisphere subpolar domains.For the pCO 2 /TA products we also see a zonally oriented local surface DIC maximum across 40°-45°S that can be distinguished from a local minimum over the Antarctic Circumpolar Current (ACC) region between 50° and 60°S, with this feature being more weakly present in the GOBMs.The differences in the seasonal amplitude of DIC (GOBMs minus pCO 2 /TA products) in Figure 5c indicates a negative bias of the GOBMs that is relatively consistent over large scales, and that is particularly pronounced at northern high latitudes.
The zonal average of the amplitude of the seasonal cycle (winter-minus-summer) in DIC is shown in Figure 5d.
Additionally the amplitude of the surface ocean DIC seasonal variability for the two three-dimensional DIC climatologies, namely MOBO-DIC and NNGv2LDEO are also shown.The MOBO-DIC product tends to have a seasonal amplitude that is somewhat larger than the surface DIC products over most latitudes, except over the Southern Ocean and over 30°-45°N.NNGv2LDEO is also larger than the GOBMs over most latitudes, also with an exception to the north of 50°N but not over the Southern Ocean.
The same collection of products considered as zonal means in Figure 5d are shown as biome averages in Figure S4 in Supporting Information S1, along with the OCIM product.This further illustrates that there is a discrepancy between the weaker DIC seasonal amplitude of GOBMs relative to not only the pCO 2 /TA products but also to MOBO-DIC and NNGv2LDEO, although the seasonal DIC amplitude of MOBO-DIC is smaller than that of all other observation-based products (and thus closer to the GOBMs) in the high latitude biomes (NA-SPSS, NP-SPSS, and SH-SS).Interestingly, for both large permanently stratified biomes (NH-STPS and SH-STPS) the late summer seasonal minimum in DIC concentrations for the MOBO-DIC and NNGv2LDEO products precedes that of the pCO 2 /TA products.

pCO 2 Seasonality Differences Between Observational Products and GOBMs
To systematically attribute discrepancies in the pCO 2 seasonal cycle between GOBMs and pCO 2 products to their drivers, we begin with an analysis of the relative importance of the thermal and nonthermal drivers of pCO 2 seasonal variability following the framework of Fassbender et al. ( 2022) (see Methods Section 2.6.2 as well as the final section of the Supporting Information S1 for a description).For this approach, in contrast to the widely applied deconvolution method of Takahashi et al. (1993Takahashi et al. ( , 2022)), monthly pCO 2 deviations in a given year are considered relative to pCO 2 values calculated from annual mean values of SST, SSS, DIC, TA, PO 4 and SiO 4 .This approach is chosen instead of considering pCO 2 deviations from the annual mean of pCO 2 , as working with the annual mean can obscure seasonally asymmetric pCO 2 responses to its drivers.
Within this analysis we consider for each grid point the absolute value of the differences between the winter and summer seasonal averages over 2014-2018.The difference between the thermal and nonthermal pCO 2 amplitudes is shown for the observational products (Figure 6a) and for the GOBMs (Figure 6b).We have opted to show differences between thermal and nonthermal amplitudes, represented in units of μatm, rather than the dimensionless ratio of these two quantities, as this provides some context for the magnitude of the difference between the ranges of pCO 2 extremes for the thermal and nonthermal components.For both panels, positive values (red) indicate that the thermal pCO 2 component dominates the seasonal cycle amplitude and negative values (blue) indicate that the nonthermal component dominates for this climatological diagnostic.
We see consistency between the pCO 2 products (Figure 6a) and the GOBMs (Figure 6b) in that the thermal pCO 2 component is dominant over the subtropics, with this dominance being even stronger for the GOBMs.Over the northern subpolar gyres and the Southern Hemisphere circumpolar regions there is a pronounced disagreement between the pCO 2 products and the GOBMs.Specifically, for the pCO 2 products the nonthermal pCO 2 component is dominant in these biomes, while the GOBMs exhibit a relatively strong thermal component in these regions.These discrepancies between the pCO 2 products and GOBMs are also shown for each GOBM and pCO 2 product separately in Figure 6c, revealing the relatively strong agreement between the pCO 2 products and the disparity between the GOBMs that was also in evidence in Figure 4 for pCO 2 itself averaged over biome scales.
As we have identified in Figures 4 and 5, there is a systematic discrepancy between GOBMs and pCO 2 products, with GOBMs indicating a weaker climatological seasonal cycle in DIC variations, which is consistent with a more pronounced dominance of the thermal component of pCO 2 seasonality highlighted in Figure 6.Next, we evaluate whether the nonthermal pCO 2 component discrepancy between GOBMs and pCO 2 products can be explained by the smaller DIC seasonal variability in the GOBMs.
We consider the drivers of the climatological pCO 2 seasonal cycle during 2014-2018 over the six aggregated biomes in Figure 7.This is done through the application of Equation 1, where the details of the Taylor Series decomposition terms used in Equation 1 are provided through Equations S1-S9 of the Supporting Information S1, and the coefficients by biome are given in Table S1 in Supporting Information S1.The left column of Figure 7 shows pCO 2 climatologies for the pCO 2 products (blue) and the GOBMs (green).This serves as a reminder that over the NA-SPSS, NP-SPSS and SH-SS domains, there are fundamental discrepancies between the GOBMs and the pCO 2 products in their representation of both the phase and amplitude of pCO 2 seasonality, while in NH-STSS, NH-STPS, and SH-STPS regions the amplitude of the seasonal cycle of pCO 2 is consistently larger in GOBMs.In order to understand the cause of these discrepancies, a decomposition of the month-to-month changes in pCO 2 into its main drivers (SST, SSS, nDIC, and nTA) is presented in the second and third column of Figure 7.
For the pCO 2 products (Figure 7, second column) the SST term dominates over the nDIC term in the subtropics, while the nDIC contribution is dominant in the subpolar biomes of the Northern Hemisphere and high latitude regions of the Southern Hemisphere.For each of these six biomes, we see that the GOBMs have a smaller nDIC contribution (third column of Figure 7) when compared to the pCO 2 products (second column of Figure 7).
It is worth noting for the observational pCO 2 products (second column in Figure 7) that although the effect of the spring bloom in NA-SPSS is marked by a short nDIC-driven minimum in the pCO 2 tendency in April, for NP-SPSS the nDIC-driven pCO 2 tendency is sustained from April through July (emphasis on the blue nDIC tendency line in each case).For the GOBMs (third column in Figure 7) there is no sustained dominance of the nDIC contribution relative to SST over most of the seasonal cycle for NA-SPSS and NP-SPSS.For both of these biomes the GOBMs simulate weak phytoplankton blooms in late spring as indicated by the nDIC tendency term for both cases, with the discrepancies in amplitude being larger than discrepancies in phase.This results in an SST dominance of the pCO 2 tendency through the summer months, thereby driving pCO 2 to peak 6-7 months earlier than that of the pCO 2 products for each of the two biomes.This represents the impact of the weaker representation of the seasonal amplitude of surface nDIC variations for the GOBMs.
In the fourth column of Figure 7, the difference between the GOBM contributions (third column) and the pCO 2 product contributions (second column) allows us to identify nDIC discrepancies as the underlying term that leads to the principal different tendencies in pCO 2 seasonal amplitude over biome scales.In other words, over these aggregated scales the most important driver that causes the disagreement between the pCO 2 seasonal cycle of the GOBMs and the pCO 2 products lies in the processes that control nDIC and equivalently for biome scales DIC seasonality.It is important to emphasize that the GOBMs provide a distinct advantage over coupled ESMs in studying disagreements between modeled and observation-based pCO 2 seasonality, because their climatological seasonal evolution of SST is effectively imposed from observations.This allowed us to identify and understand the potential magnitude of pCO 2 seasonality biases that are imposed by a seeming systematic bias with GOBMs to underestimate the importance of DIC seasonality.It is worth recalling here that the DIC from the pCO 2 /TA products is more consistent with the independent products (MOBO-DIC and NNGv2LDEO) than with the GOBMs (Figure 5d), lending credence to the discrepancies highlighted in Figure 7.
A candidate mechanism that could potentially account for why the GOBMs have smaller seasonal surface DIC variability than the DIC products is in the modeled seasonal variability of biological drawdown of surface DIC.However, in addressing this question we are faced with the caveat that seasonal variability in biological drawdown of surface DIC over global scales is poorly constrained by observations (Lee, 2001;Lee et al., 2002).To shed light on this question, we consider in Figure 8 the relationship between climatological (2014-2018) surface salinity-normalized DIC (nDIC) variations and mixed layer depth (MLD) averaged over the six biomes.The MLD product is derived from Argo observations, described in Section 2.6.1 of the Methods section.Figure 8 shows that seasonal MLD variations averaged across the GOBMs are of the same or larger amplitude when compared to the observational product.Systematic discrepancies toward deeper mixing in the models exists in the permanently stratified biomes, and to some extent in the NP-SPSS and NH-STSS biomes.However, during summer months when the MLD reaches a minimum, strong nDIC reductions are consistently more pronounced in the observations than in the GOBMs in all biomes.This suggests that the surface summer DIC drawdown, likely associated with biological drawdown of DIC, is underestimated in GOBMs.This could reflect a spurious exhaustion of nutrients by early summer in GOBMs.Investigating these questions further is not possible with the tools available through RECCAP2.However, given their importance for pCO 2 , it is strongly advised that process understanding of the potential shortcomings of models in this regard be further investigated.

Attribution for Decadal Increases in pCO 2 Seasonality
Here we adopt the framework of Fassbender et al. (2022) to address both mechanistically and quantitatively the drivers of enhanced pCO 2 seasonal amplitude.Specifically, this method allows one to deconvolve the relative importance of these drivers, namely the transient invasion flux of anthropogenic carbon (denoted as the C ant contribution) and other drivers (denoted collectively as the climate contribution).As the attribution methodol- Calculations are performed first individually on each product or GOBM, and subsequently averaged to arrive at a mean respectively for the observation-based products and GOBMs.The differences in the drivers between the pCO 2 products and the GOBMs are shown in the right column-note the differences in scale.The drivers are identified via Taylor Series decomposition (Equations S1-S9 in Supporting Information S1) using PYCO2SYS to calculate the coefficients for DIC and TA, and using constant values for SST and SSS (dlnpCO 2 /dT = 0.0423°C −1 and dlnpCO 2 /dlnS = 1) as in Sarmiento and Gruber (2006).For pCO 2 itself (left column) the monthly dots are positioned mid-month, mid-way between tick marks.For tendencies (columns two and three) dots are located at the end of the month, so on tick marks.SUM in columns two and three represents the sum of the four underlying tendencies.ogy of Fassbender et al. (2022) requires not only surface pCO 2 but also surface DIC and TA concentrations (in addition to SST, SSS, PO 4 , and SiO 4 ), we focus here on the three RECCAP2 pCO 2 /TA products that include the necessary quantities, specifically JMAMLR, OceanSODA-ETHZ, and CMEMS-LSCE-FFNN in Table 2, as well as the GOBMs listed in Table 3.For the observation-based products, monthly climatological PO 4 and SiO 4 fields were taken from World Ocean Atlas 2018 (Garcia et al., 2019).
We begin the analysis for this section on decadal changes in pCO 2 and CO 2 flux seasonality with a descriptive assessment of seasonal asymmetries in the changes in seasonal amplitude between 1985-1989 and 2014-2018 in Figure 9. pCO 2 anomalies are shown relative to the annual mean value computed from annual means of SST, SSS, DIC, TA, and nutrients (pCO 2 AM ) for both observationally based products and GOBMs, following the method described in Fassbender et al. (2022).Calculations are performed separately for each year and then averaged over the time intervals 1985-1989 and 2014-2018.This provides an estimate of the annual mean if pCO 2 were to respond linearly to its individual drivers over the seasonal cycle, making it possible to identify seasonally asymmetric changes in pCO 2 that are obscured when computing a mean of monthly pCO 2 values over an annual cycle.
We first consider the case of the subtropics as approximately 15°-40°N and 15°-40°S, where summer anomalies of pCO 2 for 1985-1989 (Figure 9a) and 2014-2018 (Figure 9b) are larger for the GOBMs than for the pCO 2 products, reflecting discrepancies in the amplitude of DIC seasonality between the GOBMs and observation-based products discussed previously (Figure 5d).Importantly, between 25° and 40°N and 25°-40°S the summer anomalies are larger than the winter anomalies over both time intervals.For the difference between the 2014-2018 and 1985-1989 intervals (Figure 9c), we see that over the Northern Hemisphere subtropical latitudes the summer changes are larger than the winter changes.In other words, notable asymmetry can be identified over the Northern Hemisphere subtropics with larger increases in summer highs than decreases in winter lows.North of 45°N there are relatively symmetric changes in summer and winter for the observation-based products and the GOBMs, but the changes for the two product classes are of reverse sign (see next paragraph) and larger for the pCO 2 products.The changes for the Southern Hemisphere are less pronounced, with some degree of asymmetry over 10°-25°S where for the observation-based products the winter lows decrease more than the summer highs increase.There is no discernible asymmetry over the Southern Ocean.
The disagreement between the GOBMs and pCO 2 products, both for seasonal pCO 2 anomalies (Figures 9a and 9b) and their changes (Figure 9c), is largest in the northern subpolar latitudes (45°-60°N).For both the 1985-1989  1985-1989 period, 2014-2018 period, and the difference between periods (Δ) are shown in separate columns.Seasonal anomalies of surface pCO 2 are calculated relative to an annual mean computed from annual means of temperature, salinity, dissolved inorganic carbon, total alkalinity, and nutrients for each year (pCO 2 AM , see Equation 4).Different magnitudes in the multi-decadal changes for each season (panel c) reflect asymmetries in the pCO 2 seasonal cycle change.Ensemble mean GOBM changes are shown in green hues (11 models listed in Table 3, excluding OCIM).Ensemble mean pCO 2 product (Prod.)changes are shown in blue hues (same pCO 2 products as in Figure 6).Seasons are defined as JFM and JAS.
and 2014-2018 climatologies, the pCO 2 products show relatively strong positive pCO 2 anomalies in winter, and relatively strong negative pCO 2 anomalies in summer, reflecting the interplay between winter convective mixing and warm-season export production in this region where nonthermal pCO 2 effects are dominant.In contrast, over the same region the anomalies for pCO 2 in GOBMs are relatively small in both winter and summer, reflecting the near balance between thermal and nonthermal pCO 2 cycles.Although it is not explicitly shown here, it is important to keep in mind that Figures 2 and 4 revealed that there are discrepancies in the phenology or phasing of seasonal variations between the GOBMs and pCO 2 products, such that the results with the GOBMs are sensitive to the specific definition of the seasonal cycle amplitude with pCO 2 products showing indeed the largest seasonal change in winter and summer for 45°-60°N, while GOBMs show the largest change in early spring and late summer/early autumn (Figures 2 and 4).For the Southern Hemisphere region between 45° and 60°S, the GOBMs and pCO 2 products are in much better agreement than they are over the Northern Hemisphere subpolar regions over both 1985-1989 and 2014-2018, with relatively modest decadal changes for the northern subpolar regions where nonthermal pCO 2 drivers are larger than thermal pCO 2 drivers for the pCO 2 products.
The corresponding asymmetries in zonally integrated seasonal CO 2 fluxes, shown in Figures 9d-9f, are quite different from what was seen for surface pCO 2 anomalies.For both time periods 1985-1989 and 2014-2018, seasonal CO 2 fluxes are larger over the subtropics in winter (ingassing) than in summer (outgassing), which indicates that seasonal wind speed asymmetries more than compensate the summer-skewed pCO 2 seasonal asymmetry over the same latitude range (Fassbender et al., 2022).This is consistent with the winter CO 2 fluxes showing local extrema for ingassing near 40°N and 40°S where westerly winds exhibit meridional maxima.This feature is also found in the 1985-1989 to 2014-2018 decadal change.For the Northern Hemisphere subpolar regions spanning 45°-60°N, the seasonal amplification in CO 2 fluxes between 1985-1989 and 2014-2018 shows enhanced ingassing during summer for the pCO 2 products, whereas the GOBMs show enhanced winter uptake over these latitudes.For the circumpolar regions of the Southern Hemisphere spanning 45°-60°S, both the GOBMs and the pCO 2 products broadly exhibit enhanced ingassing in summer (JFM) and winter (JAS) seasons relative to the 1985-1989 period.
We next turn our attention in Figure 10 to deconvolve the role of C ant invasion and climate perturbations (e.g., changes to the ocean physical or biogeochemical state caused by anthropogenic climate change) in driving multi-decadal changes in pCO 2 seasonality between 1985-1989 and 2014-2018.The changes in the subtropics are qualitatively similar in the pCO 2 products and GOBMs (top and model rows of Figure 10), with the changes in the total seasonal cycle amplitude (Figure 10a) dominated by the C ant invasion impact on the thermal component with a small climate-driven increase in the thermal pCO 2 component (Figure 10b).However, for the subpolar regions of the Northern Hemisphere and the Southern Ocean, the pCO 2 products show that the nonthermal pCO 2 component changes driven by climate perturbations (Figure 10c) are sufficiently large to dominate the total response in these regions (Figure 10a), a feature that is not seen (or considerably weaker) for the GOBMs.C ant -induced declines in carbonate buffering in the high latitudes also work to slightly amplify the total nonthermal pCO 2 seasonal cycle changes in both the pCO 2 products and GOBMs (red lines in right column).
For the total change in the GOBM pCO 2 values (Figure 10d), the impact of climate perturbations is relatively minor so that the change in pCO 2 seasonality is largely dominated by the invasion flux of C ant .For the thermal component of pCO 2 changes in the GOBMs, both the subtropical and subpolar latitudes of the Northern Hemisphere indicate dominance of the effect of transient C ant invasion, with the Southern Ocean posing one region where this is neutralized by climate perturbations.For the nonthermal component of pCO 2 changes in the GOBMs, the impact of C ant invasion is relatively minor, so that the changes are dominated by climate perturbations, which are also small.In summary, for the GOBMs, increases in the total pCO 2 seasonal cycle amplitude are dominated by the direct influence of C ant on the thermal pCO 2 component seasonal cycle amplitude through the enhanced temperature sensitivity associated with elevated ocean pCO 2 values.

Estimating Irreducible Uncertainty in Seasonal Cycle Amplitude Changes
We now turn our attention to the emergence of decadal changes in pCO 2 seasonality, with applications at the gridpoint, zonal mean, and biome scales.

Natural Variability Uncertainty in Detecting Forced Changes in pCO 2 Seasonal Amplitude
As was emphasized in the study of Landschützer et al. (2018), natural variability in the seasonal amplitude of pCO 2 variations can obscure detection of anthropogenic trends, an effect that can be termed natural variability uncertainty.For sea surface pCO 2 , the signal (Figure 11a) and the signal-to-noise ratio (SNR) (Figure 11b) are calculated for a single LE, to illustrate patterns of emergence.We choose here the CESM2 model (Danabasoglu et al., 2020) for which the phasing of the pCO 2 seasonal cycle over the Northern Hemisphere subpolar biomes has the best correspondence with the observation-based products out of the LE simulations described in Table 4 or the CMIP6 models shown in Table 6.
According to the CESM2-LE, the pCO 2 seasonality increases are most emergent (|SNR|>1) (Figure 11b) at the grid point level over much of the subtropical Northern Hemisphere, as well as over a band spanning the subtropical convergence regions of the Southern Hemisphere.Large parts of the North Atlantic subpolar gyre exhibit a decrease in the amplitude of seasonal pCO 2 variations.This analysis is indicative of pCO 2 seasonality changes over the Southern Ocean being non-emergent above the noise level of natural variability at the grid point level over our three-decade interval.In particular, pCO 2 seasonality changes in the Antarctic Circumpolar Current region spanning 50°-60°S are non-emergent, along with some parts of the Southern Hemisphere subtropics.
We next consider the emergence characteristics for the underlying thermal and nonthermal drivers of pCO 2 seasonality as represented by the CESM2-LE.We see that the forced decadal change (ensemble mean) for the thermal component (Figure 11c) and nonthermal component (Figure 11e) are both considerably larger than the decadal change for the full pCO 2 field itself (Figure 11a).For the resultant signal-to-noise ratio, we see that the patterns for the thermal component (Figure 11d) and the nonthermal component (Figure 11f) are spatially offset from each other, and when considered together help to account for the relatively restricted regions over which there is local grid point emergence, predominantly over the thermally dominated subtropical domain, for the pCO 2 seasonal cycle.
From these analyses, we can see that within the subtropics, the emergent regions for total pCO 2 SNR (|SNR|>1 in Figure 11a) correspond to the regions where the ensemble mean nonthermal component (Figure 11e) is small.

Interpretation of Results of Attribution Analysis
Having visualized the patterns of emergence for pCO 2 seasonality changes in Figure 11 we redirect out attention to the last row of Figure 10 (panels 10g-10i), where the same set of 50 large ensemble members with the CESM2-LE  3, and the CESM2 large ensemble, respectively.Shading reflects the standard deviation across the pCO 2 products, GOBMs, and ensemble members.Seasons are defined as JFM and JAS.
are applied to interpret the zonal mean changes in the pCO 2 products (first row of Figure 10).The shading for the LE results represents natural variability uncertainty (one standard deviation about the 50-ensemble-member mean), while the shading for the pCO 2 products (Figures 10a-10c) can be interpreted to represent mapping uncertainty.These forms of uncertainty are distinct but can be used together to assess the degree to which the C ant and climate-driven signals for the pCO 2 products are detectable.We consider the means of the pCO 2 products in Figures 10a-10c to be "best estimates" of the decadal changes.Superimposing the envelope of natural variability from the LE analysis (Figures 10g-10i) onto the pCO 2 product means would suggest that the C ant and climate driven signals are emergent when they and their associated natural variability are separated from zero for their decadal change.
Under this framework, we interpret the C ant component of the pCO 2 product changes to be detectable for all regions and pCO 2 components (pCO 2 , pCO 2 T , and pCO 2 NT ).What is more surprising and compelling with the pCO 2 products is the fact that the climate-driven changes over the northern subpolar domains, with much of this signal found in the nonthermal component, cannot be explained by natural variability according to the CESM2-LE results.Further attribution of this climate-driven signal (whether it reflects processes related to mixing or biology) is beyond the scope of this study, but at the very least this underscores the importance of For the full pCO 2 seasonal cycle, (i) the signal is defined as the ensemble mean change, and (ii) the noise as the ensemble standard deviation, and (iii) the SNR indicates the confidence levels for emergence.This is done analogously for the thermal component of seasonality changes for (iv) the signal, (v) the noise, and (vi) the confidence intervals for emergence.Likewise for the nonthermal component, (vii) the signal, (viii) the noise, and (ix) the confidence intervals for emergence.Absolute magnitudes are chosen to emphasize the strength of the signal, so that the thermal and nonthermal fields in the left column do not necessarily sum to the total pCO 2 signal.The decomposition of the thermal and nonthermal components of the pCO 2 seasonal cycle is described in Section 2.6.2.
sustained seasonally resolving observations in high latitude regions to better constrain and understand such behavior.
Taken together, the attribution analysis in Figure 10 and the analyses of pCO 2 and CO 2 flux seasonal asymmetries in Figure 9 indicate that, over the broad areal extent of the subtropics, the pCO 2 seasonal amplitude increases with an asymmetric larger increase in summer are in large part due to the C ant impact on the thermal driver of seasonal amplitude.Through the effect of stronger wind speeds in winter than summer, the seasonal CO 2 flux asymmetry becomes even larger and of opposite sign over the subtropics, with large increases in winter ingassing.Over the broad expanse of the subtropics, the agreement between the pCO 2 products and GOBMs, together with the clear emergence indicated by the LE, gives us reasonable confidence for both the driver (C ant ) and the asymmetry of pCO 2 seasonality changes and elevated winter CO 2 uptake.However, for the high northern latitudes, pCO 2 products and GOBMs disagree on the drivers of the changing pCO 2 seasonality.Here, the observation-based products indicate that a rather large climate forcing is responsible for pCO 2 seasonal cycle changes.

Emergence of Forced Trends With Biome Aggregation
We next address the degree to which biome aggregation improves the emergence of pCO 2 seasonality trends between 1985-1989 and 2014-2018.To this aim, we use multiple LE simulations with five different ESMs (detailed in Table 4) and apply them to interpret the decadal changes in pCO 2 seasonal variations represented by the pCO 2 products.However, before considering multi-decadal trends we consider first in Figure 12a an analysis for the six aggregated biomes of the simulated seasonal amplitude of pCO 2 variations for a climatology over 2014-2018, both for pCO 2 products and LEs, with uncertainty ranges indicated for each case.For the observation-based pCO 2 products, the uncertainty shown represents the uncertainty across the nine pCO 2 products themselves, which is reflecting the spread across methods (mapping uncertainty) in representing gridded pCO 2 .For an individual LE, the uncertainty represents the range of natural variability represented by the collection of individual ensemble members for that specific LE.We see that for all six aggregated biomes, at most two of the five LEs have a climatological seasonal amplitude (defined as the absolute value of JFM-JAS) where the first quartile range (25%-75%, solid bar) overlaps with the first-to-third quartile range of the pCO 2 products.Over four of the six biomes (NP-SPSS, NH-STPS, SH-STPS, and SH-SS) CESM2-LE overlaps with the pCO 2 products, thereby indicating correspondence.We also see that inter-model differences in seasonal amplitude (ensemble median) are much larger than the ensemble spread (natural variability) for the individual ESMs.Overall, the analysis in Figure 12a indicates that the representation of the climatological seasonal amplitude of pCO 2 in LEs is partly inconsistent with the pCO 2 products, mirroring what was found for GOBMs in the first part of the Results section of this study.
On the other hand, the decadal changes in the amplitude of seasonal pCO 2 variations between 1985-1989and 2014-2018 (Figure 12b (Figure 12b) reveal a much better agreement (relative to Figure 12a) between the median of the observational products and the ensemble spread of the individual LEs.Most importantly, with one exception over the three subtropical biomes (NH-STSS, NH-STPS, SH-STPS) and for each LE within these regions, the range of outcomes for modeled pCO 2 seasonality changes (the range indicated by the whiskers) does not intersect with a zero trend or go negative.Thus in the subtropics the biome-aggregated changes in the seasonal amplitude of pCO 2 are predominantly emergent, with SH-SS also being largely emergent.For the NP-SPSS biome, the median decadal change calculated for each LE indicates an increase in pCO 2 seasonal amplitude, although one of the LEs is not emergent while two other LEs are near the one standard deviation threshold for emergence.This is indicative of the NP-SPSS being moderately emergent under aggregation.
For the NA-SPSS biome the decadal changes in pCO 2 seasonal amplitude are quite distinct from the other regions, in that all LEs are at least moderately emergent, but that they have divergent signs in their simulated decadal changes.Three of the LEs (CESM1, CESM2, and ESM2M) exhibit a decrease in pCO 2 seasonal amplitude, while two of the LEs (CanESM2 and CanESM5) exhibit an increase.Under the caveat that we have imposed a somewhat restrictive definition of the seasonal cycle that is optimally suited for the pCO 2 products (see Figures 2  and 4), the at least moderately emergent decrease in pCO 2 seasonal amplitude for three LEs is distinctive for   13 the NA-SPSS biomes.The study of Goris et al. (2023) suggests that such behavior may be expected from the disparate responses of the Atlantic Meridional Overturning Circulation (AMOC) for the LEs.However, further exploration of this point should be conducted while taking into consideration the sensitivity of the models to the definition of the seasonal cycle.It is interesting that one of the pCO 2 products does show zero net change between 1985-1989 and 2014-2018.The fact that the distribution of multi-decadal trends for individual LEs overlaps with the 25%-75% (bar) range for the observation-based pCO 2 products gives us a degree of confidence that the biome-aggregated changes can  1985-1989 and 2014-2018 (right) based on pCO 2 products (blue) and five large ensemble experiments.For the large ensemble experiments, the number of ensemble members used for each model is shown in parentheses in the legend.The seasonal cycle amplitude is calculated as the absolute value of biome-averaged winter-minus-summer differences (JFM-JAS and JAS-JFM for the Northern and Southern Hemisphere, respectively).The boxes show the values between the first quartile (25%) and the third quartile (75%) with a line at the median (50%), with the whiskers indicating the minimum and maximum values.Note that for the pCO 2 products, the uncertainty range reflects uncertainty between the available nine observational pCO 2 products, whereas the uncertainty range for the large ensemble simulations indicates natural climate variability as simulated by each model.For each ensemble member and data product, the change of the amplitude of the seasonal cycle for each region was measured by (i) detrending the 1985-2018 time series by subtracting a quadratic polynomial fit, (ii) calculating the seasonal cycle amplitude as the winter-minus-summer pCO 2 values for each year, (iii) averaging the five yearly values for each period, (iv) subtracting the earlier period average from the latter, and (v) calculating the spatial mean with latitude weighting.The amplitude of the seasonal cycle in 1985-1989 was obtained in step (iii).For each biome the sequencing of the colors is identical to the sequencing in the legend.
be considered emergent or detectable over this 30-year period.It should be noted, however, that the Northern Hemisphere biomes tend to be faced with a triple-whammy of elevated uncertainty, given the relatively high degree of disaccord between pCO 2 products (mapping uncertainty), the relatively elevated level of natural variability uncertainty, and the offset of the medians of the different LEs (model structural uncertainty).The LEs nevertheless bring great value to the parsing of uncertainty in that they allow us to distinguish between internal variability uncertainty and model structural uncertainty.Qualitatively they are of similar amplitude for the decadal changes (Figure 12b), but the model structural uncertainty is much larger for the 2014-2018 climatology (Figure 12a).
The Southern Ocean (SH-SS) region is characterized by relatively weak trends (Figure 12b), yet due to its relatively weak natural variability is also predominantly emergent at the biome scale despite being non-emergent at the gridpoint scale (Figure 11b).The reason for there being both weak natural variability and a weak decadal trend over the Southern Ocean is due to the combined effects of the tendency of models to have a relatively weak dominance of thermal over nonthermal seasonality drivers over the Southern Ocean, and that the amplitude of seasonal SST variability over the Southern Ocean is relatively weak.Taken together, this limited the leverage typically afforded to the invasion flux of C ant to boost pCO 2 seasonality amplitude between the 1980s and 2010s over the Southern Ocean.

Discussion
The unprecedented collection of pCO 2 products and model simulations available through the RECCAP2 project has provided an opportunity for not only synthesizing but also understanding the changing seasonal cycle in sea surface pCO 2 .The use of new observation-based pCO 2 /TA products that include DIC and TA in addition to pCO 2 has made it possible to address the following scientific questions: (a) what are the underlying causes of discrepancies in pCO 2 seasonal variability between GOBMs and observation-based products?(b) what are the drivers of decadal-timescales trends toward increasing pCO 2 seasonality?, and (c) are there asymmetries in the growth of pCO 2 and CO 2 flux seasonal variations over multi-decadal timescales that could give rise to a feedback mechanism consistent with the mechanism identified in the modeling study of Fassbender et al. (2022)?
One of the main findings of this study is that there are systematic discrepancies in the climatological seasonal amplitude of surface DIC variations between GOBMs and observation-based pCO 2 /TA products over very large spatial scales (Figure 5c) that are primarily responsible for compromising the climatological seasonal cycle of pCO 2 .More specifically, simulated DIC seasonal variability is systematically smaller for GOBMs than for DIC seasonal variability in the pCO 2 /TA products.The climatological seasonal amplitude of DIC variations can be impacted by both the seasonal amplitude of net biological removal of DIC from the surface layer and the seasonal process of entrainment associated with destratification.Although it is beyond the scope of this study to parse these processes, important information can be inferred from previous studies.Models have previously been shown to prematurely exhaust limiting nutrients in summer (Goris et al., 2018), and such behavior could impact the seasonal amplitude of surface DIC variations.
From the biome-averaged time series shown in Figure 4, we have seen that there are some important discrepancies in phasing of the climatological seasonal cycle for pCO 2 between GOBMs and pCO 2 products over biome scales.We have opted here to focus largely on a definition of seasonal cycle amplitude as the difference between winter and summer, defined as JFM-JAS (JAS-JFM) for the Northern (Southern) Hemisphere, as this was deemed to capture the main characteristics of pCO 2 seasonality over large scales, based on the pCO 2 products, as shown in the Hovmöller diagrams in Figure 2.Although our main story of seasonality changes in pCO 2 is not strongly sensitive to this definition for the pCO 2 products or the GOBMs (Figure 3), the presence of phasing differences between pCO 2 products and GOBMs in high latitudes influences the seasonal amplitude calculated according to this definition.We recommend further investigation of such phasing differences to better understand the underlying reasons for this discrepancy.
We emphasized earlier that for our purposes of a synthesis, the GOBMs potentially offer an important advantage over coupled ESMs for our seasonality study and associated attribution analysis.This is due to the fact that the SSTs of the GOBMs are constrained to correspond closely to observations through the imposition of a negative feedback in surface heat fluxes by virtue of using observed surface air temperatures (Large & Yeager, 2009).However, there is no analogous data-constraint applied for the DIC seasonal evolution, other than an indirect and much weaker effective restoring of sea surface pCO 2 to atmospheric pCO 2 over a 1-year timescale (an order of magnitude less than for SST).As the SSTs of the GOBMs effectively track the observed seasonal cycle with fidelity, our attention in this study was drawn to the dominant role of DIC seasonality discrepancies in distorting pCO 2 seasonality for the GOBMs relative to the observation-based products.
Nevertheless, as a complementary exercise to the evaluation of GOBMs in Figure 4, we consider here an analogous analysis of a collection of the 11 CMIP6 ESMs listed in Table 6.The ESMs are evaluated for their climatological seasonal cycle in pCO 2 , SST, and DIC, as well as for decadal changes in pCO 2 .For this case, we have considered changes between 1985-1989 and 2010-2014 (an interval 4 years shorter than in Figure 4) to be consistent with the fact that the historical period in CMIP6 ends in 2014.In general terms there is no fundamental difference between both classes of models.Despite the absence of observational constraints, the simulated SST seasonal cycle in ESMs is in good agreement with observations (dotted lines in central column of Figure 13), although the model spread is significantly larger as expected.The DIC seasonality of the ESMs is very similar to that of the GOBMs (weaker in models than in observation-based pCO 2 products) and the overall pCO 2 changes are very consistent between fully coupled ESMs and GOBMs (right column of Figures 4 and 13).This serves to further underscore our main finding that there is a pervasive discrepancy between ocean carbon cycle models and pCO 2 products, not limited to the RECCAP2 GOBMs, with a seasonal variability in models that is too weak relative to pCO 2 products.Interpretation of projected future seasonal changes in ESMs should therefore take into account the systematic model uncertainties highlighted in this study.As we have seen in Figure 12, a future extension of this analysis for coupled models could benefit from the inclusion of LEs, which offer a means to deconvolve for models the effects of structural uncertainty from natural variability uncertainty.
Our attribution analysis indicates that in the subtropics the GOBMs and pCO 2 products are qualitatively consistent in both the degree of seasonal asymmetry in decadal pCO 2 changes and CO 2 flux changes (Figure 9), as well as for the mechanistic drivers (Figure 10).This reflects the fact that both the GOBMs and pCO 2 products are dominated by the thermal pCO 2 component in the subtropics, albeit with the degree of dominance of the thermal pCO 2 component being exaggerated by the GOBMs.The fact that the Northern Hemisphere subtropics have a larger response than the Southern Hemisphere subtropics is expected given that the seasonal cycle in SST is larger in the north than in the south, and thereby the thermal component of seasonality is larger in the north than in the south.Thus, C ant invasion will thereby preferentially boost the amplitude of pCO 2 seasonality in the subtropics of the Northern Hemisphere relative to the Southern Hemisphere, due to the higher amplitude SST variability to the north.
An important implication of our attribution analysis is that the decadal changes in seasonal pCO 2 and CO 2 flux amplitudes are qualitatively consistent with the existence of a negative climate feedback over the historical period spanning 1985-2018.Consistent with what was found by Fassbender et al. (2022) using future projections with an LE simulation, the seasonal asymmetry in CO 2 fluxes identified here for the historical period is by inference in large part due to seasonal asymmetries in wind speeds and thereby piston velocities, rather than asymmetries in pCO 2 itself.It is our hope that the analysis provided here will motivate future work to quantify the strength of this feedback over the historical period, as doing so is beyond the scope of this study.
Two studies that evaluated the early generation of prognostic ocean biogeochemistry models for their representation of pCO 2 seasonality against observations were McKinley et al. (2006) for the case of the North Pacific and Tjiputra et al. (2012) for the case of the North Atlantic.Although an explicit comparison of the models applied by McKinley et al. (2006) and Tjiputra et al. (2012) with the RECCAP2 GOBMs is beyond the scope of this study, it is important to ask, "how far have we come?".The message from this study is that important discrepancies in pCO 2 seasonality between models and observation-based products persist approximately two decades into the era of prognostic biogeochemical models.The new suite of pCO 2 /TA products that provide time-varying surface DIC products constitute invaluable tools for identifying and attributing potential model biases over global scales.Moving forward, we recommend that efforts be devoted to understanding how complexity in physical and biogeochemical process representation in models is reflected in pCO 2 seasonality.For the case of the biogeochemical components of models, one opportunity would be to evaluate with the same circulation state a climatology of pCO 2 seasonality using different complexity of biogeochemical models, analogous to what was explored by Galbraith et al. (2015) with the BLING and TOPAZ biogeochemical models, or by Karakus et al. (2022).Moreover, it could be beneficial to link the simulated strength of key circulation metrics (e.g., Atlantic Meridional Overturning Circulation for the North Atlantic) to the pCO 2 seasonality of the GOBMs.In a study with several coupled climate models, Goris et al. (2023) were able to show that the representation of a model's pCO 2 seasonality in the North Atlantic is strongly linked to the simulated strength of the Gulf Stream volume transport.In order to address the question of how missing physical processes may impact biogeochemistry and pCO 2 seasonality, it would be beneficial to explore independently the impact of mesoscale turbulence with ocean eddies (Lévy et al., 1998) and shear-induced turbulence associated with near-inertial oscillations, ocean swells, and wind stirring (Jochum et al., 2013;Keerthi et al., 2021;Rodgers et al., 2014) through perturbation studies with individual models.Such efforts would be beneficial in identifying the processes that contribute to the discrepancies identified in this RECCAP2 synthesis effort.

Conclusions
In this study, we have presented as a contribution to the RECCAP2 effort a synthesis of how the seasonal cycle in sea surface pCO 2 and air-sea CO 2 fluxes has changed over the last three decades .Our main findings include: For a 2014-2018 climatology, we identified a systematic discrepancy in the seasonal amplitude of surface DIC variations between the GOBMs, which are too weak, and the surface pCO 2 /TA (DIC) products, such that the pCO 2 seasonal cycle in subtropical biomes is spuriously large, and both the amplitude and phase of seasonal pCO 2 variations are spuriously modulated in subpolar and Southern Ocean biomes.Similar biases are found in ESMs as submitted to CMIP6.When averaged over aggregated biomes, both pCO 2 products and GOBMs represent consistent increases in the seasonal amplitude of pCO 2 and CO 2 fluxes between 1985-1989 and 2014-2018.
Decadal increases in pCO 2 seasonal amplitude in subtropical biomes are shown to be dominated process-wise by the invasion flux of anthropogenic carbon (C ant ) for both the pCO 2 products and GOBMs.For pCO 2 seasonality changes over subpolar and circumpolar biomes, GOBMs see this change as being dominated by C ant invasion, but for pCO 2 products it is dominated rather by modulations of the climate state, which can be associated with changes in SST seasonality or in perturbations occurring in mixing or biological processes.Our analysis of LE simulations indicates that the climate-driven signal found in the subpolar biomes of the Northern Hemisphere (bold blue line in Figure 10a) is likely too large to be explained by natural variability, as estimated from the amplitude of natural variability in the LE (width of blue shading in Figure 10g).Considered together, the subtropical biomes exhibit decadal increases in CO 2 flux seasonality that are larger during winter than summer, consistent with the mechanism identified by Fassbender et al. (2022) whereby increased seasonal cycles may be promoting a negative feedback in the climate system.
We have also identified and quantified two distinct sources of uncertainty that can impact the detection of anthropogenic trends in pCO 2 seasonal amplitude.The first type is mapping uncertainty (Figures 3a and 3c), which is in essence a reducible form of uncertainty.The second type is natural variability uncertainty (Figures 11 and 12), which we identified with the LE simulations, with this being an irreducible form of uncertainty.For both cases, we identified that there are large regions of the global domain where the forced anthropogenic signal does not rise above the one standard deviation (noise) level at the gridpoint scale.Even for the idealized case of a perfectly resolving observational system where mapping error disappears, one would still need to quantify the degree of confidence one has that a local measured decadal change represents forced changes in the presence of natural variability.We tend in the subtropics to find unanimous agreement at the gridpoint scale across pCO 2 products and GOBMs on the dominant mechanism of enhanced pCO 2 seasonality, namely the C ant -driven boost to the thermally component of pCO 2 seasonality.
Thus our synthesis provides evidence that the seasonal pCO 2 amplitude has significantly increased globally, and for many regions, this increase has already emerged from the noise with biome aggregation.This highlights the importance of including seasonal changes when considering long-term acidification trends, as we will likely cross critical thresholds in some seasons, before such thresholds are crossed in the annual mean (Burger et al., 2020(Burger et al., , 2022;;Henson et al., 2017;Kwiatkowski & Orr, 2018;McNeil & Sasse, 2016).For similar reasons, we recommend that modeling resources be devoted to identifying missing processes and implementing representations or parameterizations of the pertinent missing processes in ocean biogeochemical models.
Our conclusions rest on the growing availability of pCO 2 products that include associated time-varying DIC and TA products, with this representing a significant advance for scientific inquiry into the mechanisms controlling the ocean uptake of anthropogenic carbon.The fact that the three pCO 2 /TA products (considered here as surface DIC products) that exist to date (JMAMLR, OceanSODA-ETHZ, and CMEMS-LSCE-FFNN) indicate consistent patterns of seasonal variability over large scales, which are also consistent with depth-resolving three-dimensional monthly DIC climatologies (MOBO-DIC and NNGv2LDEO) over large scales, should boost community confidence in the scientific utility of these products.We therefore strongly recommend to the pCO 2 product providers to provide DIC and TA along with sea surface pCO 2 and air-sea CO 2 fluxes.In parallel, we recommend a rigorous skill-assessment of the mapping methods that transform from SOCAT measurements to globally gridded pCO 2 products, for example, by using LEs with ESMs as a synthetic testbed similar to the work presented by Gloege et al. (2021).
+ SP-STPS + Southern Indian Ocean Southern Hemisphere seasonally stratified (subpolar and subtropical combined) SH-SS SO-STSS + SO-SPSS Note.The full names are given in the left column, the abbreviated names used throughout this study in the central column, and the original RECCAP2 biomes from which they are constructed are listed in the right column.Our six aggregated biomes are shown in the central panel of Figure 1.As has been noted in the Methods section, spatial and temporal data coverage of data products and models might further constrain the area of the aggregated biomes used for our analysis.

Figure 2 .
Figure 2. Hovmöller diagrams of zonally averaged climatological (2014-2018) pCO 2 anomalies for pCO 2 products (a) and GOBMs (e), and corresponding standard deviations (b, f); multi-decadal change in pCO 2 climatology between 1985-1989 and 2014-2018 for pCO 2 products (c) and GOBMs (g) and corresponding standard deviation (d, h).For the ensemble mean, the seasonal maxima (black triangles) and minima (white triangles) are highlighted in panels (a) and (e).Anomalies were calculated, gridpoint by gridpoint, by first fitting a quadratic polynomial through the full monthly time series over the period 1985-2018 and removing this trend from the data.Then, mean pCO 2 values over1985-1989 and 2014-2018  were calculated for each month as a deviation from the detrended mean pCO 2 over the whole 5-year period.

Figure 3 .
Figure 3. Multi-decadal changes in the pCO 2 seasonal amplitude between the two 5-year periods 1985-1989 and 2014-2018, for (a, c) pCO 2 products and (b, d) GOBMs.For panels (a) and (b), seasonal amplitude is calculated as a winter-minussummer difference, where winter and summer are defined as JFM and JAS (JAS and JFM) for the Northern (Southern) Hemisphere, respectively.For panels (c) and (d) seasonal amplitude is calculated as minimum-minus-maximum.Dotted regions indicate where the spread of the ensemble members (STD = 1 level) is larger than the ensemble mean multi-decadal change.For this analysis, the nine pCO 2 products listed in Table2and the 11 GOBMs in Table3(excluding OCIM) are used.The sensitivity to the choice of 5-year versus 10-year versus 15-year averaging for characterizing decadal changes is displayed in FigureS3in Supporting Information S1.

Figure 4 .
Figure 4. Climatological seasonal cycle anomalies for 2014-2018 for pCO 2 (left column), SST (center column, dashed) and DIC (center column, solid), and change in pCO 2 seasonal cycle anomalies between1985-1989 and 2014-2018 (right column)  over the six aggregated biomes for GOBMs (green), pCO 2 products (blue), and for OCIM (brown, dashed).The biome names as defined in Table1and Figure1are given in each panel title.Ensemble means for GOBMs and pCO 2 products are shown as thick lines while the shading indicates the standard deviation around the mean.Note the different scales of the left and right columns, and that the Southern Hemisphere sub-plots begin with the month of July.

Figure 5 .
Figure 5. Climatological seasonal amplitudes in surface DIC concentrations over 2014-2018 considered as winter-minus-summer (JFM-JAS for the Northern Hemisphere and JAS-JFM for the Southern Hemisphere) for (a) the average of three pCO 2 /TA (DIC) products (JMAMLR, CMEMS-LSCE-FFNN, and OceanSODA-ETHZ), (b) the average across 11 GOBMs, and (c) the difference between (b) and (a), where negative (positive) values indicate that GOBMs under-represent (over-represent) the climatological seasonal amplitude of DIC variability.(d) shows latitudinal averages for GOBMs (green) and data products (blue), where bold lines represent product averages, and light lines indicate individual products.Also shown in (d) are MOBO-DIC (red dashed) and NNGv2LDEO (red solid).The MOBO-DIC climatology corresponds to the period 2004 through 2017 and NNGv2LDEO to 1995.

Figure 6 .
Figure 6.Ensemble mean difference between the thermal and nonthermal pCO 2 component of its seasonal cycle amplitude (ΔA) for (a) the three pCO 2 products that include associated DIC and TA products (JMAMLR, OceanSODA-ETHZ, and CMEMS-LSCE-FFNN and (b) 11 GOBMs.(c) Zonal mean of ΔA values for GOBMs and pCO 2 products with bold lines showing the ensemble mean values.The olive-colored line shows the abiotic OCIM model, which is not included in the model ensemble mean.pCO 2 component seasonal cycle amplitudes were computed as the absolute value of the difference between the winter-minus-summer seasonal averages over 2014-2018.For each case the separation into thermal and nonthermal components follows the approach of Fassbender et al. (2022), where a carbon system calculator (CO2SYS) is used to quantify the thermal component, and thereby the nonthermal component via difference from the full pCO 2 annual cycle.Seasons are defined as JFM and JAS.

Figure 7 .
Figure7.Decomposition of the seasonality drivers for pCO 2 products and GOBMs (one standard deviation shading for left column) for all six biomes, representing the transformed Taylor Series decomposition in Equation 1.The 2014-2018 climatological seasonal cycle of pCO 2 (left column) is decomposed into month-tomonth changes due to two physical (SST (orange) and SSS (purple)) and two biogeochemical (salinity-normalized DIC (nDIC, blue) and salinity-normalized TA (nAlk, yellow)) drivers (2nd and 3rd column).The observational products are the three products that include DIC and TA, namely JMAMLR, OceanSODA-ETHZ, and CMEMS-LSCE-FFNN.Calculations are performed first individually on each product or GOBM, and subsequently averaged to arrive at a mean respectively for the observation-based products and GOBMs.The differences in the drivers between the pCO 2 products and the GOBMs are shown in the right column-note the differences in scale.The drivers are identified via Taylor Series decomposition (Equations S1-S9 in Supporting Information S1) using PYCO2SYS to calculate the coefficients for DIC and TA, and using constant values for SST and SSS (dlnpCO 2 /dT = 0.0423°C −1 and dlnpCO 2 /dlnS = 1) as inSarmiento and Gruber (2006).For pCO 2 itself (left column) the monthly dots are positioned mid-month, mid-way between tick marks.For tendencies (columns two and three) dots are located at the end of the month, so on tick marks.SUM in columns two and three represents the sum of the four underlying tendencies.

Figure 8 .
Figure 8. Seasonal evolution of the relationship between MLD and anomalies of nDIC over the six biomes for a climatology over 2014-2018 for observation-based products (blue) and GOBMs (green).The horizontal axis represents MLD, and the vertical axis represents seasonal anomalies in nDIC concentrations relative to the annual mean.The stars indicate January, and the subsequent months of March, June, September, and December are marked with 3, 6, 9, and 12, respectively.The pCO 2 products that include DIC and TA, namely JMAMLR, OceanSODA-ETHZ, and CMEMS-LSCE-FFNN, are considered along with the full suite of 11 GOBMs.Observationally based MLDs have been derived for this study from the gridded Argo-derived temperature and salinity product ofRoemmich and Gilson (2009) using a density threshold criterion(Holte & Talley, 2009) using the years 2014-2018, for consistency with the nDIC fields.

Figure 9 .
Figure 9. Mean summer (S) and winter (W) (a-c) surface pCO 2 anomalies and (d-f) CO 2 sea-to-air flux.Results for the1985-1989 period, 2014-2018 period, and  the difference between periods (Δ) are shown in separate columns.Seasonal anomalies of surface pCO 2 are calculated relative to an annual mean computed from annual means of temperature, salinity, dissolved inorganic carbon, total alkalinity, and nutrients for each year (pCO 2 AM , see Equation4).Different magnitudes in the multi-decadal changes for each season (panel c) reflect asymmetries in the pCO 2 seasonal cycle change.Ensemble mean GOBM changes are shown in green hues (11 models listed in Table3, excluding OCIM).Ensemble mean pCO 2 product (Prod.)changes are shown in blue hues (same pCO 2 products as in Figure6).Seasons are defined as JFM and JAS.

Figure 10 .
Figure 10.Mean change in the seasonal cycle amplitude (winter-minus-summer) of pCO 2 and its thermal (T) and nonthermal (NT) components from the period 1985-1989 to the period 2014-2018.Changes are broken down into the total change and the portion of the change attributed to the accumulation of anthropogenic carbon (C ant ) versus climate (Clim: changes to the ocean physical or biogeochemical state caused by anthropogenic climate change).The top, middle, and bottom rows show the ensemble mean results for the three pCO 2 products, the eleven GOBMs listed in Table3, and the CESM2 large ensemble, respectively.Shading reflects the standard deviation across the pCO 2 products, GOBMs, and ensemble members.Seasons are defined as JFM and JAS.

Figure 11 .
Figure11.Multi-decadal changes in pCO 2 seasonality between1985-1989 and 2014-2018 (left column)  and signal-to-noise ratio (SNR) maps for pCO 2 seasonality changes between1985-1989 and 2014-2018 (right column).For each case, the seasonal amplitude is defined as the absolute value of the difference between JFM and JAS averages of pCO 2 in the model, such that negative values indicate decreasing amplitude.(a) the full decadal change in pCO 2 seasonal amplitude, (b) the SNR for the full pCO 2 signal, (c) the change in the thermal component of pCO 2 seasonal amplitude, (d) the SNR for the thermal component change, (e) the change in the nonthermal component, and (f) the SNR for the nonthermal component.All fields are from the CESM2 large ensemble (50 members).For each ensemble member, the change in the amplitude of the annual cycle was measured by (i) detrending the 1985-2018 time series by subtracting a quadratic polynomial fit, (ii) calculating the annual cycle strength as the absolute value of the JFM-JAS difference in pCO 2 for each year, (iii) averaging the five yearly values for each period, and (iv) subtracting the average over the earlier period from the average over the later period.For the full pCO 2 seasonal cycle, (i) the signal is defined as the ensemble mean change, and (ii) the noise as the ensemble standard deviation, and (iii) the SNR indicates the confidence levels for emergence.This is done analogously for the thermal component of seasonality changes for (iv) the signal, (v) the noise, and (vi) the confidence intervals for emergence.Likewise for the nonthermal component, (vii) the signal, (viii) the noise, and (ix) the confidence intervals for emergence.Absolute magnitudes are chosen to emphasize the strength of the signal, so that the thermal and nonthermal fields in the left column do not necessarily sum to the total pCO 2 signal.The decomposition of the thermal and nonthermal components of the pCO 2 seasonal cycle is described in Section 2.6.2.

Figure 12 .
Figure12.Box plots of the absolute amplitude of the sea surface pCO 2 seasonal cycle for the period2014-2018 (left)  and the change in the amplitude of the pCO 2 seasonal cycle between1985-1989 and 2014-2018 (right)  based on pCO 2 products (blue) and five large ensemble experiments.For the large ensemble experiments, the number of ensemble members used for each model is shown in parentheses in the legend.The seasonal cycle amplitude is calculated as the absolute value of biome-averaged winter-minus-summer differences (JFM-JAS and JAS-JFM for the Northern and Southern Hemisphere, respectively).The boxes show the values between the first quartile (25%) and the third quartile (75%) with a line at the median (50%), with the whiskers indicating the minimum and maximum values.Note that for the pCO 2 products, the uncertainty range reflects uncertainty between the available nine observational pCO 2 products, whereas the uncertainty range for the large ensemble simulations indicates natural climate variability as simulated by each model.For each ensemble member and data product, the change of the amplitude of the seasonal cycle for each region was measured by (i) detrending the 1985-2018 time series by subtracting a quadratic polynomial fit, (ii) calculating the seasonal cycle amplitude as the winter-minus-summer pCO 2 values for each year, (iii) averaging the five yearly values for each period, (iv) subtracting the earlier period average from the latter, and (v) calculating the spatial mean with latitude weighting.The amplitude of the seasonal cycle in 1985-1989 was obtained in step (iii).For each biome the sequencing of the colors is identical to the sequencing in the legend.

Figure 13 .
Figure 13.Same as Figure 4 but with CMIP6 models as listed in Table6.Note that the climatologies for pCO 2 (left column) and both DIC and SST (middle column) are taken over the years 2010-2014 (last 5 years of the CMIP6 historical simulations) and that the change in pCO 2 seasonality (right column) is calculated as the difference between2010-2014 and 1985-1989.

Table 1 The
Six Aggregated Biomes Used in This StudyFigure1.Central panel showing biomes used for this study (see Table1for details): North Atlantic subpolar seasonally stratified (NA-SPSS), North Pacific subpolar seasonally stratified (NP-SPSS), Northern Hemisphere subtropical seasonally stratified (NH-STSS), Northern Hemisphere subtropical permanently stratified (NH-STPS), Southern Hemisphere subtropical permanently stratified (SH-STPS), and Southern Hemisphere seasonally stratified, incorporating the subpolar and subtropical components (SH-SS).Surrounding panels show time series plots of biome-integrated sea-air CO 2 fluxes (annual maximum and minimum values irrespective of the month of occurrence).pCO 2 products (blue) and GOBMs (green) are shown for both the ensemble mean (bold) and for one standard deviation (shaded).Positive (negative) values indicate outgassing (ingassing) of CO 2 .The regions in white in the central panel are not included in this analysis, in each panel, winter is designated by W, and summer by S, in order to distinguish the seasonal phasing between the subtropical and subpolar/Southern Ocean biomes.

Table 2
Note.The products are described more completely in the individal papers referenced in the table.Note that the CMEMS-LSCE-FFNN, JMAMLR, and OceanSODA-ETHZ products include associated surface DIC and TA fields that are used in this study.
. All of these products are built on surface ocean pCO 2 measurements within the Surface Ocean CO 2 Atlas

Table 2 pCO
2 Data Products Used for This Study Note.Seasons are defined as the 3-month periods JFM and JAS.Change represents the difference between1985-1989 and 2014-2018, Units are μatm.

Table 6
CMIP6 Models Evaluated Against pCO 2 Products in Figure