Drivers of Air‐Sea CO2 Flux Seasonality and its Long‐Term Changes in the NASA‐GISS Model CMIP6 Submission

Climate change will affect both the mean state and seasonality of marine physical and biogeochemical properties, with important implications for the oceanic sink of atmospheric CO2. Here, we investigate the seasonal cycle of the air‐sea exchange of CO2 and pCO2,sw (surface seawater pCO2) and their long term changes using the CMIP6 submission of the NASA‐GISS modelE (GISS‐E2.1‐G). In comparison to the CMIP5 submission (GISS‐E2‐R), we find that on the global scale, the seasonal cycles of the CO2 flux and NPP have improved, while the seasonal cycles of dissolved inorganic carbon (DIC), alkalinity, and macronutrients have deteriorated. Moreover, for all ocean biogeochemistry fields, changes in skill between E2.1‐G and E2‐R display large regional variability. For E2.1‐G, we find similar modeled and observed CO2 flux seasonal cycles in the subtropical gyres, where seasonal anomalies of pCO2,sw and the flux are temperature‐driven, and the Southern Ocean, where anomalies are DIC‐driven. Biases in these seasonal cycles are largest in the subpolar and equatorial regions, driven by a combination of biases in temperature, DIC, alkalinity, and wind speed. When comparing the historical simulation to a simulation with an idealized increase in atmospheric pCO2, we find that the seasonal amplitudes of the CO2 flux and pCO2,sw generally increase. These changes are produced by increases in the sensitivity of pCO2,sw to its respective drivers. These findings are consistent with the notion that the seasonality of pCO2,sw is expected to increase due to the increase of atmospheric pCO2, with changes in the seasonality of temperature, DIC, and alkalinity having secondary influences.

Previous modeling efforts have shown that in response to increasing atmospheric and oceanic pCO 2 , the seasonality (i.e., the seasonal amplitude, or the difference between the maximum and minimum values within a seasonal cycle) of sea surface pCO 2 will also increase (Gallego et al., 2018;Gorgues et al., 2010;Rodgers et al., 2008). These modeling studies have been supported by a recent paper showing that the seasonal cycle of pCO 2 has increased in response to the increase in annual mean pCO 2 between 1985 and 2014 (Landschützer et al., 2018). The seasonal amplitude of surface ocean pCO 2 is primarily controlled by temperature, alkalinity, and dissolved inorganic carbon (DIC; Sarmiento & Gruber, 2006;Takahashi et al., 2002). When the seasonal cycle of surface ocean pCO 2 is examined, it is often split into its thermally driven and non-thermally (DIC + alkalinity) driven components (Fay & McKinley, 2017;Landschützer et al., 2014Landschützer et al., , 2018Sarmiento & Gruber, 2006). These components often act in opposing directions. For example, in the subtropical gyres of the ocean, the summertime increase in thermally driven pCO 2 is somewhat counteracted by the decrease in pCO 2 due to increased stratification bringing less DIC to the surface, although the thermal effect dominates (Fay & McKinley, 2017;Takahashi et al., 2002). In the subpolar regions, the opposite situation occurs, with the wintertime mixing driving an increase in non-thermal pCO 2 that is somewhat, but not completely, counteracted by the decrease in temperature driving a decrease in thermal pCO 2 . Calculations using only the carbonate system and assuming equilibrium between atmospheric and ocean pCO 2 show that under a scenario of increasing atmospheric CO 2 concentration, there is amplification of both thermal and non-thermal seasonal forcing of pCO 2 (Riebesell et al., 2009). Observational evidence suggests that drivers of increases in the seasonality of pCO 2 are regionally dependent. In the subtropical gyres, an increase in the seasonality of pCO 2 is driven by an increase in the seasonality of its thermally-driven component. On the other hand, in the subpolar regions and Southern ocean, the increase in the seasonality of pCO 2 is driven by the increase in the seasonality of its non-thermally driven component (Landschützer et al., 2018). However, the previous generation of Coupled Model Intercomparison Project version 5 (CMIP5) models showed that the increase in the seasonality of pCO 2 is driven more by temperature than DIC in all regions except the Southern Ocean (Gallego et al., 2018). These previous studies motivate two important questions regarding the new generation of CMIP6 models: (1) how well do these models capture the seasonality of the CO 2 flux and surface ocean pCO 2 , and (2) what drives future changes in the seasonality of the CO 2 flux and surface ocean pCO 2 in these models? Addressing these questions rigorously requires a detailed investigation into the seasonal cycle of ocean carbon uptake, as well as its biases. Thus, the primary goal of this study is to document the NASA-GISS CMIP6 ocean carbon cycle simulation, with a focus on understanding the biases and changes in the seasonal cycle of ocean carbon uptake and pCO 2 .
In this study, we compare the modeled and observed time-averaged seasonal cycles of the CO 2 flux and ΔpCO 2 (the difference between surface ocean and atmospheric pCO 2 ), and determine how the seasonality of these CO 2 fields changes after 70 years of linearly increasing atmospheric CO 2 by 1% annually, at the time of a doubling of atmospheric CO 2 from the pre-industrial level. First, we provide an overview of the simulated and observed CO 2 flux, and also compare model and observed surface properties that can influence the bias in the flux (Section 3.1.1). We also briefly evaluate how these model biases have changed from the CMIP5 submission of the NASA-GISS model E2-R (Section 3.1.2). In sections comparing the CMIP5 to CMIP6 submissions, we refer to the CMIP5 incarnation of the model as E2-R and CMIP6 incarnation as E2.1-G. We then compare simulated (by E2.1-G) and observed seasonal cycles of the CO 2 flux, ΔpCO 2 , and sea surface properties that can drive these fields (Section 3.2.1). The averaging period for the model is near the end of the historical period, and is always aligned to the temporal coverage of every data set against which the model is compared. As with the annual climatologies, the ability of model E2.1-G to capture these seasonal cycles is compared to that of E2-R (Section 3.2.2).
In our analysis, we tease apart the drivers of the seasonal cycles of the CO 2 flux and surface ocean pCO 2 in both model E2.1-G and observations (Sections 3.2.3-3.2.4). We then difference the model means and seasonal amplitudes between the end of the historical period (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) and from the simulation in which atmospheric CO 2 increases by 1% annually after reaching twice the pre-industrial level ; Section 3.3.1). We refer to these differences as long term changes in ΔpCO 2 and the CO 2 flux. Finally, the major drivers of changes in the seasonal cycle of the CO 2 flux are assessed in both model and observations LERNER ET AL.

Model
The characteristics of the physical ocean component of the GISS GCM (GISS-E2.1-G) are described more fully in Kelley et al. (2020). Since the CMIP5 incarnation of GISS-E2-R (Romanou et al., 2013;Schmidt et al., 2014), key updates have included finer upper-ocean layering, improved implementation of mesoscale eddy transport, three-dimensional variation of mesoscale diffusivity, a higher-order advection scheme, vertical mixing driven by tidal dissipation, increased ventilation of marginal seas, and a scheme for accelerated tracer transport. As in E2-R, mass-conserving numerics and surface-flux formulations ensure consistent treatment of biogeochemical (BGC) constituents under riverine inputs, sea ice formation/melt, evaporation-precipitation, and all other model processes (Ito et al., 2020). The horizontal resolution is 1 × 1.25° in latitude and longitude respectively, and the vertical resolution is 40 quasi-constant pressure layers. Ocean only simulations of passive tracer uptake (CFC) showed good agreement with observations (Romanou et al., 2017).
The ocean carbon cycle module is an update of the one used in E2-R (Romanou et al., 2013(Romanou et al., , 2014 and which originated from the NASA Ocean Biogeochemistry Model (NOBM; Gregg & Casey 2007). It includes four phytoplankton species (diatoms, chlorophytes, cyanobacteria, and coccolithophore), four nutrient species (nitrate, silicate, ammonia, and iron), three detrital pools (nitrate/carbon, silicate, and iron) and one heterotroph species. Carbon cycling is represented through dissolved organic (DOC) and DIC and interacts with atmospheric CO 2 through gas exchange parameterization, following the CMIP6 protocol (Orr et al., 2017). Light profiles from the atmospheric radiation module are propagated underwater into the ocean and spectrally decomposed to 33 wavebands that are used to compute growth of the phytoplankton groups, sinking profiles, as well as changes in the vertical distribution of ocean temperatures due to biologically mediated absorption and scattering in the water column (Gregg & Conkright, 2002). Latto and Romanou (2018) showed that ocean carbon states estimated from the GISS E2.1-G model with the E2-R carbon cycle agreed well with observations in most regions, and that study led to the current improvements in the carbon cycle simulations.
Latest improvements to the ocean carbon cycle model include the introduction of exponential profiles for the diatom sinking and the detritus settling. A prognostic alkalinity scheme has been implemented following OCMIP2 parameterization in order to better simulate carbonate chemistry and the oceanic carbonate pump. Atmospheric dust deposition is now interactive and consistent with the model climate, and thus the ocean iron cycle is forced from the historical annual cycle climatology extracted from the GISS E2.1-G online dust simulations which include eight externally mixed minerals (illite, kaolinite, smectite, carbonates, quartz, feldspar, iron oxides, and gypsum) plus internal mixtures between seven minerals and iron oxides . The masses of free and structural iron and their fractions of total iron have been evaluated using measurements from Izaña Observatory (Pérez García-pamdo et al., 2016). Riverine delivery of BGC constituents is also recently implemented and coupled to the prognostic river runoff calculated in GISS E2.1-G as part of the simulated global hydrological cycle. The concentrations of DOC, DIC, nitrate, silicate, and Fe at all the major and many minor river mouths are obtained from an annual climatology (Da Cunha et al., 2007) and they modulate the BGC characteristics of the freshwater outflow into the ocean at these sites. The riverine contribution to alkalinity is neglected in the model.

Observations
We use a suite of observationally based climatologies to evaluate output from the historical simulation, each of which are interpolated onto the ocean model grid prior to comparison. The monthly climatologies used in this study include those for ΔpCO 2 , the CO 2 flux, DIC, alkalinity, macronutrients (nitrate and silicate), NPP, temperature, salinity, mixed layer depth, and surface wind speed. The pCO 2 and the air-sea CO 2 flux are evaluated against the climatologies of Landschützer et al. (2014), which use a neural-networking approach to map ΔpCO 2 and CO 2 flux observations from the Surface Ocean Carbon Atlas, version 2 to a gridded product with a 1° × 1° horizontal resolution. DIC and alkalinity are obtained from the surface ocean climatology of Takahashi et al. (2014). In this climatology, DIC is derived by using an inorganic carbon chemistry model that uses as inputs climatological distributions of pCO 2 , total alkalinity, salinity, and temperature. These authors obtained distributions of salinity and temperature from the World Ocean Atlas, 2009 climatology, and sea surface pCO 2 using the methodology of Takahashi et al. (2009) updated to use 2005 as the reference year. The methodology of Takahashi et al. (2009) interpolates sea surface pCO 2 from the LDEO database onto a 4° × 5° ocean grid using an interpolation method based on a 2-d advection-diffusion equation. They further derived alkalinity using a linear relationship between alkalinity and salinity from 24 different oceanic regions, and tested their salinity-derived alkalinity and model-derived DIC values against measurements from the GLODAP database (Key et al., 2004), the CARINA program (Key et al., 2010), and measurements from the LDEO database. They found that the calculated values were consistent with the measured values within their respective measurement uncertainties. The data provided in this climatology are referenced to year 2005. Note that, while more recent climatologies of DIC and alkalinity are available from GLODAPv2 (Key et al., 2015;Lauvset et al., 2016), to our knowledge, these datasets are only available as annual climatologies, not monthly climatologies which were necessary for this study.
Macronutrients are evaluated against climatologies of nitrate and silicate from the World Ocean Atlas, 2013 (Garcia et al., 2014). NPP is obtained from the Carbon-based production model, version 2 (Westberry et al., 2008). The CbPMv2 relates NPP to the product of phytoplankton carbon concentration, as estimated by remote sensing retrievals of particulate backscattering coefficients, and phytoplankton growth rates, as estimated from carbon-to-chlorophyll ratios. Here, we use NPP estimates provided by the Ocean Productivity repository: http://www.science.oregonstate.edu/ocean.productivity/index.php (Westberry et al., 2008), calculated based on optical properties retrieved from SeaWiFS measurements between 1998 and 2007.
Temperature and salinity are evaluated against the Roemmich-Gilson ARGO climatology derived from ARGO floats deployed between 2004 and 2012 (Roemmich & Gilson, 2009). The mixed layer depth is also evaluated against an ARGO-based climatology, in which vertical density differences obtained from ARGO floats deployed between 2000 and 2016 are used to determine the depth of the mixed layer (Holte et al., 2017). Finally, surface wind speed is obtained from the second Modern-Era Retrospective analysis for Research and Applications (MERRA2), a NASA atmospheric reanalysis that begins in 1980, and includes monthly gridded values with a horizontal resolution of 0.625° × 0.5° (Bosilovich et al., 2015). In comparing to the model, we take the average of the reanalysis fields each month between 2004 and 2012.
To evaluate model skill in capturing the CO 2 flux, pCO 2 , and other sea surface properties that impact these CO 2 fields, we calculated the normalized root-mean square error and bias for the global surface ocean. The normalized root mean square is: NRMSE

Experiments
We analyze results from NASA-GISS-E2.1-G simulations for the historical period, forced with observed forcing, and an idealized future forcing scenario in which the atmospheric concentration of CO 2 , starting from the pre-industrial CO 2 concentration in 1850, increases by 1% annually until it reaches double the pre-industrial value (around year 70 of the simulation). For convenience, we refer to the latter scenario as the 1% simulation. For the historical period, model output is averaged over the same period that is represented by the corresponding data set (Table 1). For example, when comparing to the climatology of Landschützer et al. (2014), which includes the monthly averaged CO 2 flux and ΔpCO 2 between 1998 and 2011, we take the model's monthly average CO 2 flux and pCO 2 fields over the same period. In this way, we attempt to limit model-data bias that may be due to differences in the timing of sampling of model output versus timing of sampling of observations. However, we note that when the averaging period is shorter than 20 years, some of the discrepancy between the model and the observed change might be attributed to internal variability. For the idealized 1% simulation, we take a 20 years average with the year in which the atmospheric CO 2 concentration reaches 2× its pre-industrial value as the central value. Model spin ups were carried out for 700 years off an equilibrated simulation of the pre-industrial ocean-atmosphere-ice system with prognostic CO 2 and about 500 years of accelerated DIC, and then the carbon cycle was again simulated online for another 100 years, for a total of ∼1,300 years.
In addition to the time-averaging, the model and observational fields are also either regionally weighted-averaged (for ΔpCO 2 , mixed-layer depth [MLD], DIC, alkalinity, nitrate, and silicate), or regionally integrated (for CO 2 flux and NPP). We ensure that for each property, only grid-cells that contain both simulated and observed values are included in the regional average or regional integration. For the weighted averages, properties with concentration units are weighted by the mass of water in the model grid cell (ΔpCO 2 , DIC, alkalinity, nitrate, and silicate), and MLD is weighted by the surface area of the model grid cell. The regionally integrated fields (CO 2 flux and NPP) are further normalized by the number of grid cells in the region for which both modeled and observed values exist, so that regionally-averaged carbon fluxes (in Pg C/yr) are calculated. The regions considered are the subpolar North Atlantic (45°N-60°N and 75°W-0°E), the subtropical North Atlantic (15°N-45°N and 75°W-0°E), the equatorial Atlantic (15°S-15°N and 75°W-0°E), the subtropical South Atlantic (60°S-45°S and 75°W-15°E), the subpolar North Pacific (45°N-60°N and 140°E-65°W), the subtropical North Pacific (140°E-110°W), the equatorial Pacific (15°S-15°N and 150°E-75°W), the subtropical South Pacific (45°S-15°S and 150°E-75°W), and the Southern Ocean (60°S-45°S and 180°E -180°W). While not encompassing the entire ocean, the regions are meant to represent either major basins which are important carbon sinks (e.g., the North Atlantic and the Southern Ocean), or areas of the largest biases in the CO 2 flux (e.g., the subpolar regions and the equatorial Pacific). So that our examination of the seasonality of the flux and pCO 2 is not influenced by seasonal sea ice, we only consider latitudes equatorward of 60°N and 60°S. We also remove all remaining grid cells that include seasonal sea ice from the regional averages. Our motivation for being restricted to ice-free waters is two-fold. First, in our analysis of the drivers (Sections 3.2.3-3.2.4), we do not consider the influence of changes in sea ice coverage on the air-sea flux, being mainly interested in the BGC and thermal drivers of the seasonality of the CO 2 flux and ΔpCO 2 . Second, the model's diagnostic pCO 2,sw output was erroneously scaled by a factor of (1-(% sea ice cover)). This error caused the model's diagnostic for pCO 2,sw to be much lower than expected in regions of seasonal sea ice.
Finally, when discussing seasonality in specific regions, "boreal" will be used to specify northern hemisphere and equatorial winter (December, January, and February), spring (March, April, and May), summer (June, July, and August), and fall (September, October, and November). Likewise, "austral" will be used to LERNER ET AL.  specify southern hemisphere winter (June, July, and August), spring (September, October, and November), summer (December, January, and February), and fall (March, April, and May). Seasons will not be specified as "boreal" or "austral" when referring to regions in both hemispheres.

Comparison to Observations
We find that pCO 2,sw shows better agreement with observations in the subtropical regions and Southern Ocean than in the equatorial and subpolar regions ( Figure 1a). In the subtropical and Southern Ocean regions, the model shows a difference that amounts to ∼±5% of the observed pCO 2,sw , while the difference is ∼±25% of the observed pCO 2,sw in the equatorial Pacific, subpolar North Atlantic, and subpolar North Pacific. In the subtropical regions, the positive pCO 2 bias is consistent with the positive temperature biases, since the higher model temperature decreases the solubility of pCO 2 in seawater. The bias in pCO 2,sw is also consistent with the negative NPP bias, which causes the model to draw down less DIC than in observations and hence increases the DIC available to be converted to pCO 2 .
In the subpolar gyres and Southern Ocean, model pCO 2,sw is generally less than observed (by ∼25% in the subpolar regions, and ∼5% in the Southern Ocean; Figure 1a). In the subpolar North Pacific, the negative temperature bias partly contributes to the negative bias in pCO 2 , with higher solubility driving pCO 2 LERNER ET AL.
10.1029/2019MS002028 6 of 33 downward in the model compared to observations (Figure 1f). The large bias in the mixed layer depth in the subpolar North Atlantic (Figure 1d) largely contributes to the pCO 2 bias. This is because while there is a positive bias in alkalinity and DIC throughout the water column in this region, the model only underestimates the vertical gradient in DIC ( Figures S1 and S2). Thus the increase in surface DIC in the model due to wintertime mixing is less than that in observations, leading to a more negative difference between DIC and alkalinity in the model than in observations. Since DIC and alkalinity have opposing effects on pCO 2,sw , with DIC increasing pCO 2,sw and alkalinity decreasing pCO 2,sw , the smaller difference between DIC and alkalinity in the model leads to negative bias in pCO 2,sw in this region. While the Southern Ocean does not exhibit the negative temperature bias in the northern subpolar regions, it does exhibit a positive NPP bias ( Figure 1e) and negative DIC bias, suggesting the larger than observed productivity is leading to enhanced DIC drawdown and reduced pCO 2,sw .
Like the subtropical regions, the equatorial Atlantic exhibits a positive pCO 2,sw bias (∼20%) that appears associated with a positive temperature bias (Figure 1a and 1f). On the other hand, in the equatorial Pacific, there is a negative bias in pCO 2,sw (∼25%) that appears associated with a positive bias in NPP ( Figure 1e) and temperature ( Figure 1f). In particular, the region of large pCO 2,sw bias in the model off the coast of Peru is associated with the region in the equatorial Pacific where the NPP and sea surface temperature (SST) biases are largest (Figures 1e and 1f), suggesting that the reason for the discrepancy between the observed and model pCO 2,sw in this region is a combination of the model's overestimation of productivity and underestimation of DIC upwelled from waters below the thermocline.
Turning to the CO 2 flux, the largest biases are in the equatorial Pacific, where the model has weak outgassing compared to observations, and in the subpolar North Atlantic and the Southern Ocean, where CO 2 uptake in the model is too large compared to observations. The model also underestimates the CO 2 flux in the subtropical gyres (North and South Atlantic, North and South Pacific) and overestimates the outgassing in the equatorial Atlantic. Additionally, there are a few local discrepancies between the direction of CO 2 gas exchange. Along the Peruvian margin, observations show outgassing of CO 2 , while the model simulates strong uptake and, in the Southern Ocean south of 60°S, observations show either no net flux or slight outgassing of CO 2 , while the model shows strong uptake.
The biases in pCO 2,sw should at least in part control the biases in the CO 2 flux, as the air-sea exchange of CO 2 is partly a function of the air-sea gradient in pCO 2 . The relationship between the CO 2 flux (F CO2 ) and ΔpCO 2 follows: where k (T, Ws), α (T, S), T, S, Ws, and sice are the piston velocity, the solubility of CO 2 , temperature, salinity, the surface wind speed, and fractional sea ice coverage, respectively. Biases in any one or a combination of the above factors, along with the air-sea gradient in pCO 2 , may contribute to the biases in the CO 2 flux. To determine the extent to which the biases in the CO 2 flux are related to the biases in pCO 2,sw , we provide scatterplots of the two CO 2 fields in each region ( Figure S3). We find a significant positive relationship between the biases in pCO 2,sw and the biases in the CO 2 flux in most regions. The one exception is the equatorial Atlantic, where biases in wind speed may be just as important in controlling the flux (see Section 3.2.4). The close relationship between the CO 2 flux and pCO 2,sw biases also explains the reversal in sign of the CO 2 flux near the Peruvian margin and Southern Ocean compared to observations. In both of these regions, the positive biases in NPP increase the drawdown of DIC, decreasing pCO 2,sw . Furthermore, underway measurements of pCO 2,sw suggest that high pCO 2,sw values are associate with cold, nutrient-rich, and oxygen depleted upwelled waters (Boutin et al., 1999;Copin-Montégut & Raimbault, 1994). This is because upwelled waters, particularly from low-oxygen waters that have been subject to intensive remineralization, are enriched in DIC, and their vertical transport to the surface increases pCO 2,sw (Feely et al., 2006;Sutton et al., 2014;Takahashi et al., 2009). In the Southern Ocean, wintertime upwelling in the Antarctic Divergence zone has a similar effect; bringing DIC-enriched waters to the surface to increase pCO 2,sw (Gruber et al., 2019;Lovenduski et al., 2007). Since the positive biases in temperature in both regions suggest reduced upwelling along the Peruvian margin and Southern Ocean, the reduction of both processes may serve to further increase the negative pCO 2,sw bias in these regions. These biases induce air-sea pCO 2 gradients that are positive (atmospheric pCO 2 greater than ocean pCO 2 ) in the model, whereas in observations they are negative. Notice that in the Southern ocean, the observed CO 2 flux is near-zero, so that it would only take a small bias in pCO 2,sw for the model's Southern Ocean to become a region of atmospheric CO 2 uptake. Table 2 lists the NRMSE and bias for each of the fields examined in this study. For pCO 2,sw , we find that the root mean square difference is only ∼10% of the observed annual mean value. The bias is similarly small, being ∼10 μatm. On the other hand, we find a large difference between model and observed fluxes: the root mean square difference between simulated and observed fluxes is ∼ 3× larger than the mean observed CO 2 flux ( Table 2). The bias in the CO 2 flux is much smaller, ∼−0.1 g C/m 2 /yr, which amounts to ∼2% of the annual mean observed flux. These differences suggest that even though regional discrepancies in the CO 2 flux are large, the model does a better job at capturing the global CO 2 flux.
The normalized root mean square differences are larger for the CO 2 flux than the remaining BGC fields, and are much larger for the MLD, NPP, nitrate, and silicate than for temperature, salinity, alkalinity, DIC, and wind speed. The biases are small for all fields, except for (1) wind speed, where the bias approaches 20% of the observed annual mean, and (2) for nitrate, silicate, and NPP, where the biases are at least 50% of their respective observed annual means.

Changes from CMIP5
Table 2 also lists the global NRMSE and bias for fields from E2-R, the NASA-GISS CMIP5 submission (this model is described in further detail in Romanou et al. (2013)). The NRMSE and bias in pCO 2,sw have increased by ∼33% and ∼22% from E2-R to E2.1-G, respectively. While the NRMSE for the CO 2 flux has increased by ∼1.5, the bias has been reduced by >90%. Thus model skill in reproducing pCO 2,sw has deteriorated overall, whereas for the CO 2 flux, the model has a reduced global bias at the expense of increased regional biases. These increased regional biases are most pronounced in the subpolar regions and Southern Ocean, where E2.1-G has stronger CO 2 sinks, in the equatorial Atlantic, where it has a stronger CO 2 source, and in the equatorial Pacific (Figures 1b, 1c, 2a, and 2b), where it has a weaker CO 2 source. Moreover, the NRMSE and bias in DIC, alkalinity, and macronutrient concentrations in E2.1-G show large increases, from 50% to 300% of these respective metrics in E2-R. On the other hand, NPP shows slight improvement in E2.1-G, with a decrease in the NRMSE of ∼6% and bias by ∼10%. The main improvements in NPP are in the Southern Ocean and equatorial Pacific, which show higher productivity compared to E2-R and similar productivity compared to (though slightly higher than) observations (Figures 1e and 2d).

Seasonality during Observational Periods
In this section, we provide an overview of the time-averaged seasonal cycles of the CO 2 flux and ΔpCO 2 in each of the regions examined in this study. We also discuss the biases in both of these seasonal cycles and use LERNER ET AL. DIC, dissolved inorganic carbon; MLD, mixed-layer depth; NPP, net primary production.

Table 2 Annual Mean Values, Normalized Root Mean Square Errors (NRMSE), and Global Bias Between Model Surface Fields and Observations Discussed in This Study
an analysis based on first-order Taylor series expansions to attribute these biases to one or a combination of biases in the fields that impact the CO 2 flux and ΔpCO 2 . To this end, we also show the seasonal cycles of the properties presented in Figure 1, including the MLD, sea surface temperature, sea surface salinity, NPP, sea surface DIC, sea surface alkalinity, and surface wind speed. While not directly causing changes in ΔpCO 2 , the roles of nitrate and silicate may also be important in seasonally and regionally limiting NPP, and thus drawing down sea surface DIC and pCO 2,sw . Hence we also present these macronutrients in Figures S4-S6.

Subtropical Gyres
The simulated seasonal cycle of ΔpCO 2 and the CO 2 flux shows good agreement with the observations (Figure 3a), with wintertime ΔpCO 2 and CO 2 flux minima and summertime ΔpCO 2 and CO 2 flux maxima. The ΔpCO 2 and CO 2 flux seasonal cycles are out-of-phase with the seasonal cycle of mixed layer depth, wind speed, and DIC (Figures 3c and 3i), and in-phase with sea surface temperature (Figures 3e and 3g). The seasonal cycles of alkalinity and salinity are small in both the model and the observations. Contrary to a large bias in the annual mean (discussed in the previous section), the NPP seasonality is realistic in the northern hemisphere basins of the Atlantic and the Pacific but it is out of phase from the observations in the respective southern hemisphere basins.

Equatorial Regions
In the equatorial Atlantic and Pacific, the simulated seasonality of ΔpCO 2 ( Figure 4b) is more pronounced than in observations. In both of these regions, the model ΔpCO 2 minima occur ∼4 months earlier than in observations. The month of maximum CO 2 uptake in the model also occurs ∼4 months earlier than in observations in both regions. Unlike for ΔpCO 2 in the equatorial Pacific, the CO 2 flux in this region shows weaker seasonality in the model than in observations. The Positive values represent outgassing. The other panels show biases between the E2-R and the observed annual climatologies of (c) pCO 2,sw (μatm) and (d) net primary production (NPP, g C/m 2 /yr). NPP in the model fails to capture the observed seasonal variations in both equatorial regions ( Figure 4d). In the equatorial Pacific, the NPP maximum is ∼6 months out of phase with observations, while in the equatorial Atlantic its seasonal amplitude is much lower than observed.

Subpolar Regions and Southern Ocean
In both the subpolar North Atlantic and Pacific, the seasonal cycles of both ΔpCO 2 and the flux are quite different than in observations. Observations show ΔpCO 2 maxima in late boreal winter, while the model maxima are shifted by ∼6 months compared to the observations. The model also produces ΔpCO 2 minima in boreal fall in both regions that are not observed. The corresponding observed CO 2 flux shows minima twice per year, in May and in October, and maxima in February and August. While the August CO 2 flux maximum is reproduced by the model in both regions, the model cannot reproduce the February maximum or May minimum in either region. The October minimum in the CO 2 flux is captured by the model in the subpolar North Pacific, but not the subpolar North Atlantic.
In the subpolar North Atlantic, the model shows a decrease in ΔpCO 2 from late boreal summer to boreal fall/winter, whereas in observations ΔpCO 2 continues to increase through boreal fall and winter. While the model and observations show similar changes in temperature during the period (Figure 5e), the model shows a smaller change in DIC, by about ∼30 μmol/kg, while observations show a ∼50 μmol/kg change. The larger increase in observed DIC causes ΔpCO 2 to increase in observations, whereas in the model ΔpCO 2 slightly decreases due to an increase in solubility that overcompensates the increase in DIC. Interestingly, the model shows a smaller DIC change from boreal summer to winter despite having a much larger increase in the MLD. The difference in boreal summer to winter DIC change is explored further in Section 4.1. The model decrease in ΔpCO 2 in boreal winter, combined with the decrease in wind speed, results in stronger uptake in boreal winter versus summer ( Figure 5a). In the subpolar North Pacific the large discrepancy LERNER ET AL.
10.1029/2019MS002028 10 of 33 2), and solid lines are the simulated values from the historical run averaged over the period of the respective observations. In the legend, "stNatl" is the subtropical North Atlantic, "stNPac" is the subtropical North Pacific, "stSPac" is the subtropical South Atlantic, and "stSPac" is the subtropical South Pacific.  Figure 3, but for the equatorial regions. In the legend "EqAtl" is the Equatorial Atlantic, and "EqPac" is the equatorial Pacific.
Figure 5. Same as Figure 3, but for the subpolar regions and Southern Ocean. In the legend, "spNAtl" is the subpolar North Atlantic, "spNPac" is subpolar North Pacific, "SOc" is the Southern Ocean.
between modeled and observed boreal wintertime ΔpCO 2 ( Figure 5b) coincides with the slightly larger than observed mixed layer depth ( Figure 5c) and the period when the model positive alkalinity bias is largest (Figures 5g and 5h). DIC and alkalinity have opposite effects on ΔpCO 2 ; while increasing DIC increases ΔpCO 2 , increasing alkalinity decreases ΔpCO 2 . The overestimation of alkalinity that is brought to the surface leads to lower ΔpCO 2 and stronger uptake than in observations.
In the Southern Ocean, ΔpCO 2 and the CO 2 flux show similar seasonal patterns and biases. The air-sea flux of CO 2 indicates maximum uptake in late austral summer/early fall (Figures 5a and 5b) and minimum in late austral winter. The timings of the peaks in the seasonal cycles of both CO 2 fields are shifted by about a month in the model compared to the observations, while the amplitudes of these seasonal cycles are also larger in the model compared to those in observations. There is also a shift in the timing of the peak in primary production, but unlike for the CO 2 fields, model peak NPP occurs earlier, rather than later, in the year compared to observations ( Figure 5d). The seasonal cycle of the remaining fields is reproduced well by the model, including the mixed layer depth, sea surface temperature, DIC, and wind speed (Figures 5c, 5e, 5g, and 5i). Both model and observations show almost no seasonal variations in alkalinity (Figure 5h).

Changes in Seasonality from CMIP5
In this section, we present the seasonality of the CO 2 flux and pCO 2,sw in E2-R and compare these seasonal cycles to those from observations and E2.1-G. The seasonal amplitudes and timing of seasonal extrema of the CO 2 fields are presented in Figures 6c-6f, while the seasonality of other fields are shown in Figures S8-S9. For ease of comparison of annual mean to seasonal cycle properties, we also show the annual means for each region in Figures 6a-6b, and Figure S7. For pCO 2,sw , the (N)RMSE and biases of the seasonal cycle (seasonal amplitude and timing of seasonal extrema) in E2.1-G are nearly the same as for E2-R (Table 3; for the timing of extrema, we use the RMSE instead of the NRMSE, since timing for all fields are in units of months). One exception is the bias in the timing of seasonal pCO 2,sw extrema, which has been reduced by ∼10% in E2.1-G from E2-R. The seasonal cycle of the CO 2 flux, on the other hand, has generally improved in E2.1-G compared to that of E2-R. The (N)RMSE and bias of the seasonal amplitude of the flux are smaller in E2.1-G than E2-R by ∼33% and ∼20%, respectively. The RMSE and bias for the timing of the seasonal CO 2 flux extrema have increased, but by only 1%-2%. Overall, Table 3 indicates that on a global scale the model's skill in capturing the seasonal cycle of pCO 2,sw remains largely unchanged between E2-R and E2.1-G, while its skill in capturing the seasonal cycle of the CO 2 flux has improved. The largest improvements (change in bias from E2-R to E2.1-G > 10%) in the CO 2 flux seasonal amplitude are in the subpolar regions, subtropical North Pacific, subtropical South Pacific, and Southern Ocean, while model skill in capturing the seasonal amplitude of the flux has decreased in the equatorial regions and subtropical South Atlantic. The bias in the seasonal extrema timings have been reduced from E2-R to E2.1-G in the subtropical North Atlantic, subtropical North Pacific, subtropical South Atlantic, and Southern Ocean while it has increased in the subpolar regions and subtropical South Pacific. The seasonal amplitude of pCO 2,sw also shows an improved fit to observations in the subpolar regions, while biases in the seasonal amplitudes of this field have increased in the equatorial Atlantic and subtropical South Atlantic (Figure 6b). The largest improvements in the seasonal extrema timings of pCO 2,sw are in the subpolar North Atlantic and Southern Ocean, while the timings have deteriorated in the subpolar North Pacific.
Interestingly, while the NRMSEs and biases of annual mean DIC and alkalinity are larger in E2.1-G than E2-R, the NRMSEs of their seasonal amplitudes are slightly smaller (by < 10%; Table 3). The timings of their seasonal extrema also show a better fit to observations in E2.1-G. However, this improvement is small, amounting to at most a 12% reduction in the bias of the timing for DIC (less than this for the other seasonal extrema timing metrics for DIC and alkalinity). Moreover, the biases in the seasonal amplitudes of both DIC and alkalinity show dramatic increases; for DIC the bias is inflated by a factor of ∼35 (from a very small bias to one that is ∼20% of the mean seasonal amplitude of DIC; Figure S8), while for alkalinity, the bias increases by a factor of ∼3. Nitrate and silicate show an overall worse fit in their seasonal cycles to observations in E2.1-G versus E2-R. The biases and RMSEs of the seasonal extrema, as well as NRMSEs of the seasonal amplitudes, for nitrate and silicate are largely unchanged between E2-R and E2.1-G. However, the biases of the seasonal amplitudes have increased by ∼21% for nitrate and ∼12% for silicate. Thus, for macronutrients, DIC, and alkalinity, the seasonal amplitude biases show changes (increases) that are >10%, with the remaining metrics showing comparatively small changes. Overall, then, Table 3 indicates a general decline in model skill for capturing the seasonal cycles of these four fields. Finally, E2.1-G displays only about half of the bias in NPP as that of E2-R, suggesting a pronounced improvement in the seasonal cycle of this field. However, the improvements for the other seasonal cycle metrics for NPP are more marginal (<10%). Figure S8 shows that this improvement in the seasonal amplitude bias masks important regional variability in changes to model skill; while the seasonal amplitude of NPP has been improved in some regions (subpolar North Atlantic, subtropical North Atlantic, equatorial Atlantic, and equatorial Pacific), it has deteriorated in others (subtropical South Atlantic and Southern Ocean).

Drivers of Seasonal pCO 2
The seasonality of ΔpCO 2 varies widely between regions both in the observations and in the model. Over the course of the seasonal cycle, the magnitude (by over a factor of 4 in some regions) as well as the sign of ΔpCO 2 can change so that a region may switch from an atmospheric CO 2 sink to source. pCO 2,sw is a function of temperature, salinity, DIC, and alkalinity. The seasonality of any of these drivers can potentially impact pCO 2,sw seasonality. In this section, we determine the main drivers of the seasonality (or monthly anomaly from the annual mean) of pCO 2,sw , which itself is critically important in determining the air-sea CO 2 flux. To determine the main drivers of surface seawater pCO 2,sw anomalies, we construct a linear, first order Taylor series expansion for the monthly anomalies in pCO 2 (Doney et al., 2009;Gallego et al., 2018;Landschützer et al., 2018;Lovenduski et al., 2007;Mogollón & Calil, 2018): where i is an index for months (i = JAN,FEB…DEC), and H.O.T. stands for higher-order terms neglected in this analysis. The δ terms represent deviations from the annual mean, that is, are obtained by differencing the (climatological) monthly (X i ) and annual mean fields ( X ), so that δX i =  X X . The partial derivatives are obtained by (i) calculating pCO 2,sw offline, (ii) adding a perturbation to each driver that amounts to 0.1% of their annual mean surface value, and recalculating the pCO 2,sw offline (Mogollón & Calil, 2018), and (iii) taking the difference between perturbed and unperturbed pCO 2,sw . Both the model derived and observed drivers are used to assess the right-hand side terms in (Equation 2 DIC, dissolved inorganic carbon; NPP, net primary production.

Table 3 Normalized Root Mean Square Errors and Global Biases Between Model and Observed Seasonal Amplitudes (s.amp.) and Timings of Seasonal Extrema
cycle of pCO 2,sw are those associated with the largest right-hand side terms. Similar methods based on Taylor series expansions have been used to assess drivers in the variability of pCO 2,sw and the CO 2 flux (Doney et al., 2009;Latto & Romanou, 2018;Lovenduski et al., 2007;Mogollón & Calil, 2018). Offline calculations of pCO 2,sw are performed using the OMIP-BGC carbonate chemistry protocols, which estimate pCO 2,sw from * 2 CO (the sum of carbonic acid and dissolved CO 2 ) and the solubility of CO 2 by (i) calculating the solubility of CO 2 as a function of temperature and salinity (as for the CO 2 flux), and (ii) calculating * 2 CO from pH and DIC, where pH is obtained by solving the pH-alkalinity equation with the SolveSaphe algorithm (Munhoven, 2013;Orr et al., 2017).
In each region, we find that temperature and DIC play the dominant role in driving pCO 2,sw . While the temperature and DIC effects often oppose each other, they are seldom of the same strength such that alkalinity partially influences pCO 2,sw . The relative importance of these three drivers depends on the region and time of year, and varies between the model and observations.
In the subtropical regions (Figures 7b, 7d, 7f, and 7h), temperature drives the wintertime minimum and summertime maximum of pCO 2,sw in both model and observations. However, the role of the DIC-driven anomalies, which act in opposite direction to the temperature-driven pCO 2,sw anomalies, is underestimated in the model in the subtropical North Atlantic, South Atlantic, and South Pacific, leading to stronger seasonality in pCO 2,sw in the model compared to observations in these regions (Figure 3b). On the other hand, in the subtropical North Pacific, the contributions of DIC and temperature to the pCO 2,sw anomalies agree with observations, although alkalinity contributes slightly more to the boreal summertime maximum and wintertime minimum in ΔpCO 2 in the model than in observations in this region.
LERNER ET AL.  In the subpolar regions (Figures 7a and 7e) and the Southern Ocean (Figure 7i), DIC plays a much larger role in driving the seasonal cycle of pCO 2,sw in both the model and in observations, being commensurate to the role of temperature in the subpolar regions (Figures 7a and 7e) and playing a larger role than temperature in the Southern Ocean (Figure 7i). However, the model can often underestimate the DIC-driven pCO 2,sw anomalies, resulting in differences between the model and observed seasonality of pCO 2,sw . For example, in late boreal summer in the subpolar North Atlantic, DIC drives negative pCO 2,sw anomalies in observations, while in the model temperature drives positive pCO 2,sw anomalies that overcompensate for the negative DIC driven pCO 2,sw anomalies, so that the model shows a boreal summertime maximum in pCO 2,sw and ΔpCO 2 (Figure 5b). In boreal winter, both the model and observations exhibit a DIC-driven pCO 2,sw maximum, but this DIC-driven maximum is underestimated by the model.
Similarly, in the subpolar North Pacific in boreal summer the DIC driven minimum in pCO 2,sw in the model is compensated by temperature and alkalinity effects on pCO 2,sw , which drive a maximum in model pCO 2,sw (Figure 5b). In observations, however, the alkalinity driven pCO 2,sw anomaly is small in boreal summer, and the DIC driven anomaly is greater than the temperature driven anomaly, so that overall the observed pCO 2,sw is lower on average during boreal summer. In the Southern Ocean, both the model and observations show DIC driven minima in pCO 2,sw in austral summer and maxima pCO 2,sw in austral winter, although the model DIC driven component in both seasons is greater in magnitude than observed. As elaborated in Section 4.1, we speculate that the bias in the DIC-driven component of the pCO 2,sw anomalies is due to overestimation of NPP in austral summer and overestimation of mixing of DIC from subsurface waters in austral winter.
Driver contributions in the equatorial regions are more complex. The dominant component of the seasonal cycle of pCO 2,sw in these regions varies with time of year and differs between model and observations. In the equatorial Atlantic, temperature and DIC both drive a pCO 2,sw maximum in early boreal spring (maximum delta pCO 2 in Figure 4b), which is partially compensated by the effect of alkalinity driving pCO 2,sw downward during the same season (Figure 7c). In August, temperature drives a pCO 2,sw minimum in the model in late boreal summer, whereas in observations, the temperature effect is roughly compensated by the effect of alkalinity, so there is neither an observed pCO 2,sw minimum nor maximum in August (Figure 4b). Instead, there is an observed pCO 2,sw minimum in November, driven by negative anomalies in observed temperature and DIC during this period. In the model, however, there is a positive temperature anomaly during this period which, combined with a positive alkalinity anomaly, drives a local pCO 2,sw maximum in late boreal fall. In the equatorial Pacific, the temperature, DIC, and alkalinity driven pCO 2,sw anomalies are generally smaller in magnitude in the model than in observations (Figure 7g). In boreal spring, temperature drives a pCO 2,sw maximum in both the model and observations. In early boreal fall, however, there is a mainly temperature driven minimum in pCO 2,sw in the model, while in observations boreal fall minima in alkalinity and maxima in DIC drive positive pCO 2,sw anomalies. Importantly, in the equatorial Pacific, the drivers of the observed seasonal cycle of pCO 2,sw must be viewed speculatively. The climatology of Takahashi et al. (2014) does not include DIC and alkalinity in the equatorial Pacific between 8°S and 8°N, due to the strong interannual variability in the measurements driven by El Ninõ and La Ninã events. Thus, the averages computed in Figure 7g are biased to drivers of pCO 2,sw outside of this equatorial band (in the model and in observations, since model averages were computed only using grid cells where observations are available). For example, the decrease in observed alkalinity in boreal fall, or "spike" in DIC in October, may not accurately reflect the seasonality of DIC and alkalinity of the entire equatorial Pacific, but only in the latitude bands outside of 8°S and 8°N. Given this observational uncertainty, it is unclear how well our model captures the seasonality of pCO 2,sw in this region.

Drivers of the Seasonal CO 2 Flux
As in the previous section, we evaluate the main drivers of the seasonality of the CO 2 flux using a Taylor series analysis. In this section, we only consider temperature and salinity effects imparted on the flux through the piston velocity and solubility, as their effects on ΔpCO 2 were detailed in the previous section. Similarly, we only consider the direct effect of wind speed on the piston velocity, as this analysis cannot account for wind-driven transport processes that also impact the air-sea exchange of CO 2 . Both the piston velocity and solubility follow the protocols set for OMIP-BGC (Orr et al., 2017), and the calculation of the piston velocity uses the same quadratic gas transfer formulation as that used by Landschützer et al. (2014). The effect of sea ice is not considered, since in this analysis, we exclude grid cells where seasonal sea ice is present.
To determine the dominant drivers of the CO 2 flux during the seasonal cycle, we employ a Taylor series expansion of the monthly anomalies in the flux: The procedure to estimate the right-hand side terms is largely the same as for pCO 2,sw . For the CO 2 flux anomalies, the partial derivatives are obtained by (i) calculating the CO 2 flux offline using the model's CO 2 flux routine, which follows the protocols for OMIP-BGC (Orr et al., 2017); (ii) adding a perturbation to each driver that amounts to 0.1% of their annual mean surface value, and recalculating the CO 2 flux offline, and (iii) taking the difference between perturbed and unperturbed CO 2 flux. Figure 8 shows the monthly anomalies in the CO 2 flux (δF CO2 as computed from each term in Equation 3).
The largest contribution to the flux anomaly reveals the dominant mechanism that determines the flux anomaly each month. In all regions, monthly anomalies in the CO 2 flux are driven mainly by wind speed and ΔpCO 2 (Figure 8). Specifically, in the subtropical gyres (Figures 8b, 8d, 8f, and 8h)  Ocean ( Figure 8i) the anomalies are nearly entirely driven by anomalies in ΔpCO 2 . In the North and South Atlantic (Figures 8b and 8d) and Pacific (Figures 8f and 8h), in both model and observations, the summertime maximum and wintertime minimum in the CO 2 flux are driven by the maximum and minimum, respectively, in ΔpCO 2 . However, in observations, the ΔpCO 2 driven anomaly in winter is smaller in magnitude than that in the model, which at least partly explains why the seasonality in the observed CO 2 flux is less than that in the model (Figure 3a). Similarly, the austral summer minimum and winter maximum CO 2 flux in the Southern Ocean are driven by ΔpCO 2 in both the model and observations, although the wind speed appears to have a larger role in driving the seasonal variability of the flux in the model than in observations (Figure 8i).
In the subpolar and equatorial regions, the CO 2 flux anomalies are driven by both ΔpCO 2 or wind speed depending on time of year, although the model and the observations do not always agree. In the subpolar North Atlantic in the model, the (late) boreal summertime maximum CO 2 flux (Figure 5a) is driven mostly by wind speed (Figure 8a). However, the contribution from the ΔpCO 2 driven flux anomalies in summer are not negligible. In boreal winter months, the simulated flux of CO 2 is more negative than in observations, due to the wind driven CO 2 flux anomalies being much larger in magnitude in the model than in observations. In the subpolar North Pacific, the model does not capture well the drivers of the monthly flux anomalies ( Figure 8e). Here, the simulated seasonal cycle of the CO 2 flux (maximum during boreal summer, minimum during boreal fall) is driven primarily by wind speed, and is opposite in phase to that of the observations, where the boreal summertime minimum and wintertime maximum flux are driven primarily by ΔpCO 2 (Figure 8e).
In the equatorial Atlantic (Figure 8d), the simulated flux has a realistic dependence to wind speed and ΔpCO 2 . The boreal springtime maximum CO 2 flux (Figures 4a and 4b) is driven by a boreal springtime maximum in ΔpCO 2 , but this effect is more pronounced in the model than in the observations (Figure 8c). However, in boreal summer, the model shows a CO 2 flux minimum driven by a ΔpCO 2 minimum, whereas observations show the flux to be higher than average (but not at a maximum) in boreal summer due a combination of above average ΔpCO 2 and above average wind speed.
Finally, in the equatorial Pacific, the model flux seasonal cycle is somewhat shifted compared to observations (Figures 4a and 4b) and is dominated by ΔpCO 2 (Figure 8g), while in observations, the flux seasonal cycle is controlled by the seasonal cycle of the wind speed. The discrepancy in the wind speed effect can be explained by considering that the CO 2 flux scales linearly with ΔpCO 2 and with the square of the wind speed (Equation 1). Thus the sensitivity (i.e., partial derivative) of the flux with respect to wind speed scales linearly with both ΔpCO 2 and the wind speed. Because in the equatorial Pacific both the wind speed and the magnitude of ΔpCO 2 are systematically underestimated by the model, the sensitivity of the CO 2 flux to wind speed is smaller in the model than in observations, leading to the model underestimating the effect of wind speed on the flux in this region.

Changes in Annual Means and Seasonal Amplitudes
In this section, we compare the change in model means and seasonal amplitudes of ΔpCO 2 and the CO 2 flux between the end of the historical period (1990-2014) and the run with 1% annual increase of atmospheric CO 2 after the pre-industrial era. For the latter run, we analyze results averaged over a 20 year period around the time of doubling of atmospheric CO 2 (i.e., simulation years 60-80). As in Section 3.2.2, we define the seasonal amplitude as the maximum minus the minimum of a field's seasonal cycle in a certain period. The changes in model means and seasonal amplitudes of the CO 2 fields during each period are examined in each of the nine regions ( Figure 9). We provide the changes in the model mean and seasonal amplitudes of all the surface properties in Figures 3-5 in Figure S10. The full seasonal cycle of the CO 2 fields, as well as the other aforementioned surface properties examined are provided in Figures S11-S13.
All regions show a decrease in ΔpCO 2 , with the largest change in the equatorial Pacific, and the smallest change in the subtropical North Pacific (Figure 9b). The changes in the seasonal amplitude of ΔpCO 2 are larger than the mean changes everywhere except the tropical regions (Figure 9b). The regions with the larg-est changes to the seasonal amplitude of ΔpCO 2 are the subpolar North Atlantic, subtropical North Atlantic, and Southern Ocean. In the subpolar and subtropical North Atlantic (Figures S11a and S12a), the change in the seasonal amplitude is driven by an increase in the boreal summertime maximum and a decrease in the boreal winter minimum, while in the Southern Ocean it is driven by a decrease in the austral summer minimum ( Figure S13a). Consistent with the increase in ΔpCO 2 , all regions show a tendency toward stronger uptake (Figure 9a), with the largest changes occurring in the subtropical North Atlantic, subpolar North Atlantic, equatorial Pacific, and Southern Ocean. Again, changes in the seasonal amplitude of the CO 2 flux are much larger than changes to the annual mean value in all except the equatorial regions (Figure 9a).

Drivers of Long Term Change in Seasonality
Since changes in the seasonal amplitude are found to be larger than changes in the mean state, we proceed to examine the drivers of the changes in the seasonality of pCO 2,sw and the CO 2 flux. To this end, we expand the Taylor series of the change in pCO 2,sw between the 1% and the historical (hist) simulations. The procedure for this Taylor series expansion follows previous studies on drivers of changes in the seasonality of the pCO 2,sw (Gallego et al., 2018;Landschützer et al., 2018). For example, for a driver X for the pCO 2,sw flux, we calculate the difference as: where δ represents the monthly change in the flux or its drivers and Δ represents the change between the two simulations. Equation 4 describes the difference in the monthly anomaly of the pCO 2,sw , as driven by property X, between the two experiments.
Continuing with our example using the pCO 2,sw anomalies, the change in the anomalies driven by any given variable X can be further decomposed following a second Taylor series expansion: Here, the first term on the right-hand side represents the change in the monthly anomaly of pCO 2,sw driven by X due to the change in the sensitivity of pCO 2,sw to X, and the second term on the right hand side represents the change in pCO 2,sw due to the change in the monthly anomaly in X. This analysis is applied to each driver of pCO 2,sw (T, S, DIC, and Alk) and to each driver of the CO 2 flux (T, S, Ws and ΔpCO 2 ). This decomposition enables us to identify the drivers responsible for the changes in the seasonality of ΔpCO 2 and the CO 2 flux, as well as whether changes in the sensitivity to the drivers or the drivers themselves influence the change in the seasonality of the CO 2 fields.
LERNER ET AL.

10.1029/2019MS002028
19 of 33 We find regional variability in the main driver of the changes in the pCO 2,sw anomalies ( Figure 10). In the subtropical regions (Figures 10b, 10d, 10f, and 10h), the main driver of the changes in the pCO 2,sw anomalies is temperature, which tends to increase the pCO 2,sw anomalies in the summer and decrease the anomalies occurring in winter. The largest contribution to the change in the temperature driven pCO 2,sw anomalies is from the change in the sensitivity of pCO 2,sw to temperature, as opposed to the change in temperature anomalies themselves.
In the subpolar North Pacific and the Southern Ocean, the main driver of the changes in the pCO 2,sw anomalies is DIC (Figures 10a, 10e, and 10i), with the largest increase in the anomalies occurring in the winter for each basin and the largest decrease in the anomalies occurring in the summer. The largest contribution to the change in the DIC-driven pCO 2,sw anomalies in these two regions is from the change in the sensitivity of pCO 2,sw to DIC, as opposed to the change in DIC itself. In the subpolar North Atlantic, the changes in the DIC driven pCO 2,sw anomalies are also larger than the changes in pCO 2,sw anomalies due to the other drivers ( Figure 10a). However, the combined change in the pCO 2,sw anomalies due to alkalinity and temperature are larger than the DIC driven changes in the pCO 2,sw anomalies, so that in this region the anomalies show the maximum increase in boreal summer and the anomalies show the maximum decrease in boreal winter. As in the subtropical regions, in the subpolar regions and Southern Ocean these changes are due to changes in the sensitivity of pCO 2,sw to DIC, alkalinity, and temperature, with changes in temperature, DIC, and alkalinity themselves playing secondary roles (Figures 10a, 10e, and 10i).
LERNER ET AL.
10.1029/2019MS002028 20 of 33 In the equatorial regions, multiple drivers are important in influencing changes in pCO 2,sw . In the equatorial Atlantic in boreal winter, the increase in the DIC driven pCO 2,sw anomalies is compensated by the decrease in the alkalinity driven pCO 2,sw anomalies, whereas in boreal summer, the increase in the alkalinity driven pCO 2,sw anomalies is compensated by the decrease in the DIC and temperature driven pCO 2,sw anomalies ( Figure 10c). In the equatorial Pacific, the temperature driven increase in the pCO 2,sw anomalies in boreal spring and decrease in pCO 2,sw anomalies in boreal summer are only partly offset by a concurrent decrease (increase) in the DIC driven pCO 2,sw anomalies in boreal spring (boreal summer). This results in slightly larger changes in the seasonality of pCO 2,sw in the equatorial Pacific than the equatorial Atlantic, though in both regions changes in the seasonality are small compared to the other regions (Figures 10c and 10g). Like the other ocean regions, changes in seasonality in the equatorial regions are driven by changes in the sensitivity of pCO 2,sw to its drivers, as opposed to changes in the drivers themselves.
Turning to the CO 2 flux, in all regions except the equatorial regions, the changes in the CO 2 flux anomalies are driven by ΔpCO 2 (Figures 11a, 11b, 11d, 11e, 11f, 11h, and 11i). In each of these regions, the greatest contribution to the long term change in the ΔpCO 2 driven CO 2 flux anomalies is from the long term change in the monthly ΔpCO 2 anomalies themselves. The long term change in the sensitivity of the CO 2 flux to ΔpCO 2 , on the other hand, plays a relatively minor role in affecting the long-term change in the ΔpCO 2 driven CO 2 flux anomalies. In fact, we see that the sensitivity of the flux to ΔpCO 2 changes little in the future. Changes in ΔpCO 2 act to decrease the flux in the winter and increase it in the summer in the Atlantic regions and subtropical North and South Pacific. This is in contrast to the subpolar North Pacific and Southern Ocean regions, where changes in the ΔpCO 2 lead to increases in the flux during winter.
LERNER ET AL. In the equatorial regions, the changes in the flux anomalies are driven by both ΔpCO 2 and wind speed (Figures 11c and 11g). However, the drivers of the change to the wind-driven CO 2 flux differ between the two equatorial regions. In the equatorial Pacific during boreal summer and late fall, the change in the sensitivity of the CO 2 flux to wind speed, as opposed to the change in the wind speed anomalies themselves, mainly drives the change in the wind-driven CO 2 flux anomalies. However, during the rest of the year in the equatorial Pacific, and throughout the entire year in the equatorial Atlantic, the changes in the sensitivity to wind speed and wind speed anomalies play commensurate roles in driving the change to the wind-driven CO 2 flux anomalies. The opposite is the case for the ΔpCO 2 driven flux anomalies; the change in the flux anomalies are driven by the change in ΔpCO 2 anomalies. The maximum change in the flux anomalies occurs in boreal summer, with the flux anomalies decreasing due to an increase in the sensitivity to wind speed in June/July and to a decrease in the ΔpCO 2 anomalies in August/September. Note that the changes in the CO 2 flux anomalies are much smaller in the equatorial regions than the other regions, showing that there is little change in the seasonality of the CO 2 flux in the equatorial regions.

Discussion
We find that in the simulated seasonal cycle of pCO 2,sw , temperature tends to be the dominant driver in subtropical latitudes, the contributions of temperature and DIC are of similar magnitude in the northern subpolar regions, and DIC tends to be the dominant driver in the Southern Ocean (Figure 7). Such a latitudinal contrast is not as apparent for the seasonal cycle of the CO 2 flux, which is controlled by ΔpCO 2 in the subtropical gyres, Southern Ocean, and subpolar North Pacific, but is driven nearly equally by wind speed and ΔpCO 2 in the subpolar North Atlantic (Figure 8). In the equatorial regions, both the seasonal cycle of pCO 2,sw and the CO 2 flux have drivers that are highly dependent on the time of year (Figures 7 and 8). We also find that the long term change in the seasonal cycle of pCO 2,sw is controlled by the same drivers that control its seasonal cycle at the end of the historical period (e.g., temperature in the subtropics and DIC in the Southern Ocean). Moreover, we find that the change in the sensitivity of pCO 2,sw to the drivers, as opposed to changes in the seasonal cycles of the drivers themselves, control the long term change of the seasonal cycle of pCO 2,sw . For the CO 2 flux, we find that in all regions except the tropics, the long term change in the seasonal cycle of the CO 2 flux is controlled by changes in the seasonal ΔpCO 2 anomalies. However, these predicted long term changes must be viewed with caution, since the model's seasonal cycles and their drivers show major biases in multiple regions when confronted with observations. Below, we discuss the biases in the seasonal cycles in each of the drivers, focusing on the regions with the largest biases (equatorial Pacific and Atlantic, subpolar North and South Pacific, and Southern Ocean), and highlight processes that may produce these biases in the model. We then discuss changes in mechanisms driving the CO 2 flux from E2-R to E2.1-G, and speculate on changes in the model that may have led to different drivers between the two model versions. Finally, we analyze the mechanisms underlying the long term change of the sensitivity of pCO 2 .

Discrepancies in the Drivers of CO 2 Seasonal Cycles
In the model, the largest biases in the seasonal cycles of pCO 2 and the CO 2 flux occur in the non-subtropical gyre regions, including the equatorial Pacific and Atlantic, the Southern Ocean, and the subpolar North Atlantic and Pacific. However, the biases are driven by different processes depending on the region considered. Table 4 summarizes the main drivers of the seasonal variability in pCO 2,sw and the CO 2 flux in both the model and in observations. The dominant driver is reported in each region in boreal winter, boreal spring, boreal summer, and boreal fall.

Equatorial Regions
In the equatorial Pacific, much of the discrepancy between the model and observed CO 2 flux seasonality can be attributed to differences in the seasonality of ΔpCO 2 itself. Depending on the time of year, the model shows large biases in the temperature-driven, DIC-driven, and alkalinity-driven pCO 2 anomalies. For ex-ample, in late boreal fall, the model's temperature-driven pCO 2 anomalies are much smaller (closer to zero) than in observations (Figure 7g). The smaller than observed anomalies in temperature may be caused by the model underestimating upwelling in the eastern equatorial Pacific. This underestimation appears to be partially caused by the model having a more stratified eastern equatorial Pacific than in observations, with the model having warmer than observed waters above the thermocline and colder than observed waters below the thermocline ( Figure S14). In addition, notice that the maximum wind speed occurs in boreal summer in the model but it is sustained in boreal summer and fall in observations (Figure 4i). This suggests that wind-driven upwelling of cooler waters, which drives down pCO 2,sw , is restricted to a shorter period of time in the model versus observations.
The observed DIC maximum in boreal late summer/early fall (Figure 5g) is consistent with other modeling efforts examining the seasonal variability of pCO 2,sw in the equatorial Pacific (Valsala et al., 2014;Wang et al., 2015), and has been attributed to the eastern equatorial upwelling during this part of the year. That the model cannot reproduce the DIC-driven pCO 2,sw maximum in boreal fall (Figure 7g) may therefore also be a reflection of the model underestimating the contribution of upwelling to surface DIC. However, if upwelling were to drive the observed fall DIC maximum, one would also expect alkalinity to be maximum during this period, since upwelled waters should be enriched in both DIC and alkalinity. Instead, we find that in observations, alkalinity is below average in boreal fall (Figure 4h). Since there is no observed NPP maximum in boreal fall, the observed lower-than average alkalinity may reflect the preferential growth of phytoplankton that are particularly efficient at removing alkalinity from the surface ocean, such as cocco-LERNER ET AL.

Table 4 Dominant Drivers of Anomalies in the CO 2 Flux (Outside Parentheses) and pCO 2,sw (Inside Parentheses) in Each region
lithophores. The observation based seasonal cycles of alkalinity and DIC, however, are highly uncertain, given that their monthly climatologies do not include observations in the central equatorial pacific ).
In our model, the seasonal variability of DIC and alkalinity is small compared to observations in the equatorial Pacific (Figures 4g and 4h). Similarly, the seasonal variability in the contribution of alkalinity and DIC to the monthly pCO 2,sw anomalies in the model is much smaller than observed (Figure 7g). To determine whether the lack of seasonal variability in model DIC and alkalinity is related to other sea surface properties (e.g., NPP), we re-calculate the average seasonal cycles of all properties shown in Figure 4, but only at locations outside of the 8°S to 8°N band, where climatological values of DIC and alkalinity are present ( Figure S15). Aside from DIC and alkalinity, we find reduced seasonal variability in NPP and sea surface temperature compared to observations. Since primary production draws down both DIC and alkalinity, reduced seasonal variability in primary productivity may partly explain the reduced variability of DIC and alkalinity. The reduced variability in SST, on the other hand, may reflect a reduction in the variability of upwelling, which would reduce the contrast in SST during periods of maximum and minimum upwelling. Thus, we speculate that reduced seasonal variability of both vertical transport and NPP contributes to the lack of seasonal variability in DIC and alkalinity in this region.
Aside from reduced seasonal amplitudes, the seasonality of the model's DIC and alkalinity is also out of phase with observations. In the model in boreal summer, upwelling appears to bring DIC and alkalinity to the surface (when the model reaches a temperature minimum; Figure S10e), resulting in DIC-driven and alkalinity-driven pCO 2 maxima and minima, respectively (Figure 7g). This is in contrast to observations, which show the period of maximum upwelling of DIC appears to be in boreal fall rather than summer (as suggested by the temperature minimum and DIC maximum; Figures S15g and S15h). In boreal fall, the model shows maxima and minima in the alkalinity and DIC-driven pCO 2,sw anomalies, respectively (Figure 7). These extrema are consistent with the model having (i) a lack of upwelling in boreal fall and (ii) having an NPP maximum in boreal fall (Figure 4d), which drives down both alkalinity and DIC. In contrast to the model, observations show lower than average NPP in boreal fall. This increase in model NPP, which is not observed, is likely due to (i) the model systematically underestimating nitrate throughout the year, and (ii) having a nitrate maximum associated with the model's late boreal summer upwelling ( Figure S5). The systematic underestimation of nitrate makes some regions of the equatorial Pacific, such as the eastern upwelling region, nitrate limiting, when observations indicate that the eastern upwelling region of the equatorial Pacific should be iron limiting (Feely et al., 2002). Thus, in the model, the pulse of nitrate that occurs in boreal summer ( Figure S5) appears to promote the higher than average productivity in boreal summer and fall that is not observed.
There is a stark contrast in the role of wind speed in determining the seasonality of the CO 2 flux in the model and in observations. Throughout most of the seasonal cycle, wind speed dominates the seasonal cycle of the flux in observations but ΔpCO 2 is the main driver in the model. The role of the model's systematically weaker winds and lower (in magnitude) ΔpCO 2 in lowering the magnitude of the wind driven CO 2 flux anomalies has already been mentioned in Section 3.2.4. However, there are also discrepancies in the seasonal variation of the wind speed between the model and observations. For example, in early boreal fall, the model has weaker than average winds, which drive a negative flux anomaly, whereas observations have stronger than average winds, which drive a positive flux anomaly (Figure 8g). Landschützer et al. (2013) reported that in the equatorial Atlantic, the thermal-and non-thermal-driven pCO 2,sw changes roughly compensated one another, so that the net seasonal change is small. This is the case in our model as well, although there are clear biases in the alkalinity and DIC-driven components of pCO 2,sw . The largest biases in the components of the seasonal cycle of pCO 2 in the equatorial Atlantic occur in the alkalinity and DIC-driven components during boreal summer. In observations, the alkalinity and DIC minima correspond to an NPP maximum in August, suggesting that productivity causes the extrema in the observed alkalinity and DIC-driven pCO 2,sw in late boreal summer. In the model, however, no summertime NPP maximum is apparent, likely because the model has consistently low nitrate in this region throughout the year ( Figure S5). Instead, model DIC and alkalinity-driven monthly pCO 2,sw anomalies reach their minimum and maximum, respectively, earlier in boreal summer, although the cause of these extrema are not apparent from the sea surface properties examined in Figure 4.
There is broad similarity between the CO 2 flux and ΔpCO 2 observed and modeled seasonal cycles in the equatorial Atlantic, including the CO 2 flux and ΔpCO 2 trough in boreal fall and the CO 2 flux and ΔpCO 2 peak in boreal spring (Figures 4a and 4b). These extrema are similar to features found in previous modeling studies (Wang et al., 2015) and in other data sets (Lefèvre et al., 2013;Padin et al., 2010;Parard & Boutin, 2010). However, the drivers of the CO 2 flux are different between model and observations throughout the year (Table 4; Figure 8c). Both wind speed and ΔpCO 2 driven components of the flux are larger in magnitude in the model than in observations, though are generally of the same sign except in late boreal fall. This is consistent with the model having ΔpCO 2 that is systematically larger in magnitude than in observations, so that the sensitivity of the CO 2 flux to wind speed is larger in the model than in observations.

Subpolar North Atlantic and Pacific
In the subpolar North Atlantic and Pacific, both the model and observations show DIC-driven pCO 2,sw anomaly maxima in boreal winter and temperature-driven anomaly maxima in boreal summer (Figure 7; Table 4). In the summer, DIC and temperature drive pCO 2,sw anomaly extrema in the opposite direction to those in boreal winter. These extrema are due to convective mixing bringing cold deep waters enriched in DIC to the surface (Miller et al., 1999), with the subsequent decrease in DIC in spring due to biological draw-down (Landschützer et al., 2014;Takahashi et al., 2002) and increased boreal summertime stratification. Key differences between the model and observations, however, are that the model shows smaller DIC-driven pCO 2,sw extrema than observed in the subpolar North Atlantic, and larger alkalinity-driven pCO 2,sw extrema than observed the subpolar North Pacific. In the subpolar North Atlantic, the model shows less seasonality in surface DIC because the vertical gradient in DIC in the top 1,000 m is much smaller than in observations ( Figure S1). The smaller vertical gradient means that although boreal wintertime mixing in the subpolar North Atlantic is larger than in observations, it does not increase DIC as much as in observations (Figure 5g). In the subpolar North Pacific, seasonality differences in pCO 2,sw are more closely linked to seasonality differences in alkalinity. Like DIC, alkalinity is enriched in deep waters, so that in the model boreal wintertime mixing brings alkalinity to the surface ( Figure 5h) and reduces pCO 2,sw . The model's negative alkalinity-driven pCO 2,sw anomaly ( Figure 7b) appears to be caused by the model having more alkalinity between 100 and 200 m than in observations ( Figure S16). Thus, although the model and observed boreal wintertime mixed layer depth are similar (Figure 5c), the model overestimates the entrainment of alkalinity into surface waters during winter in the subpolar North Pacific.
The alkalinity-driven boreal summertime pCO 2,sw maxima in the model in the subpolar North Pacific coincides with the period of maximum stratification, after the model's spring NPP maxima (Figures 5c and 5d). While increased NPP is responsible for the DIC minima in both the model and observations, only the model shows a productivity driven alkalinity minimum in boreal summer. This is likely because in the model, the processes affecting alkalinity and DIC are similar. Alkalinity calculations follow OCMIP2, which computes alkalinity based on the following assumptions: (i) alkalinity changes proportionally to the rate of change in phosphate (increasing phosphate decreases alkalinity), which itself is computed from the rate of change in nitrate scaled by the Redfield ratio, at each time step, (ii) calcium carbonate production (which decreases alkalinity) is proportional to NPP in the euphotic zone, and (iii) dissolution of calcium carbonate (which increases alkalinity) is proportional to the divergence of the downward flux of organic material (Najjar & Orr, 1999). Since biologically driven changes in alkalinity are largely proportional to biologically driven changes in DIC, it is somewhat expected that increasing NPP simultaneously decreases alkalinity and DIC in the model. However, in the real ocean, alkalinity changes are not only dependent on total productivity or the net divergence of organic material, but rather dependent on productivity of organisms chiefly responsible for alkalinity changes (e.g., particulate organic carbon forming organisms such as coccolithophore), and dissolution of calcium carbonate at depth. Thus, differences in the model and observed alkalinity-driven pCO 2,sw anomalies may partly be due to actual versus modeled biologically driven alkalinity changes, in addition to differences in the transport of deep-ocean alkalinity to surface waters.
Finally, differences between model and observed alkalinity may also be related to salinity biases. For example, in the annual climatologies, alkalinity and salinity biases appear closely associated in the Gulf stream, near Central America, and in the southern tropical Pacific (Figures 1g and 1i). However, in the majority of the subtropical regions, alkalinity, and salinity biases are opposite in sign. Moreover, in the subpolar North Pacific, there is little difference between the model and observed seasonal cycles of salinity, while observations show seasonality in alkalinity that is reduced compared to the model (e.g., no summertime minima; Figures 5f and 5h). Thus, both the annual climatologies and seasonal cycles suggest that processes controlling salinity biases (e.g., differences in surface freshwater fluxes) only partially contribute to alkalinity biases.
In both the subpolar North Atlantic and Pacific, the model flux is shown to be more influenced by wind speed than the observed flux at all times except during boreal spring (Figures 8a and 8e; Table 4). Similar to the equatorial Pacific, the higher effect of wind speed occurs because the model overestimates the magnitude of ΔpCO 2 throughout most of the year (Figure 5b), leading to a higher sensitivity of the CO 2 flux with wind speed. The lower effect of ΔpCO 2 on the CO 2 flux in the model, on the other hand, is caused by the model underestimating wind speed throughout the year (Figure 5i), which lowers the model's sensitivity of the CO 2 flux to ΔpCO 2 .

Southern Ocean
The model's seasonal cycle of pCO 2,sw and the CO 2 flux is somewhat more consistent with observations in the Southern Ocean than in the other non-gyre regions. Minimum uptake and ΔpCO 2 both occur in late austral summer, with the CO 2 flux being ΔpCO 2 (Figure 8i) driven and pCO 2,sw being DIC driven (Figure 7i). The non-thermal component of pCO 2,sw drives the seasonal cycle through mixing of DIC to the surface in austral winter and DIC drawdown by biological utilization in austral summer (Gruber et al., 2019;Landschützer et al., 2014;Mongwe et al., 2018;Pasquer et al., 2015). Mongwe et al. (2018) have shown that the CMIP5 Earth System Models had difficulty capturing the seasonal cycle in the Southern Ocean, though in different ways. North of the Polar Front (where we restrict our analysis), some models overestimated the DIC-driven component of the seasonal cycle of pCO 2,sw (3 of 10), while most models (9 of 10) overestimated the SST-driven component of the seasonal cycle. Our analysis shows that the SST driven anomalies in pCO 2,sw are roughly consistent with observations, while the model generally overestimates the magnitude of the DIC-driven pCO 2,sw anomalies in austral winter (Figure 7i). In austral winter, this overestimation may be due to the model overestimating the mixing of DIC below ∼100, since above 200 m (the average wintertime mixed layer depth in the Southern Ocean) in the Atlantic Sector of the Southern Ocean between 45°N and 56°N, the model overestimates DIC concentrations in subsurface waters ( Figure S17). The apparent overestimation may also be an artifact of the lack of wintertime observations of DIC and pCO 2 in the Southern Ocean Lauvset et al., 2016). In austral summer, the model underestimation of the DIC-driven pCO 2,sw anomalies may be due to the model having greater austral summertime NPP than in observations (Figure 5d).
In our analysis, we have separately examined the biases in the seasonal cycles in each region. However, a bias common to all regions in our model is the underestimation of nitrate. The underestimation of nitrate may partly explain the negative NPP bias in multiple regions, including the subtropical and equatorial regions. Our model currently computes nitrogen fixation supported production based on iron and light availability, and does not include the effects of oxygen inhibition (Dunne & John, 2013). Furthermore, denitrification is not explicitly represented; instead, nitrate is removed in a given grid cell based on (i) the vertical integral of nitrogen fixation, and (ii) the ratio of nitrate in a grid cell to the vertical integral of nitrate. Nitrification is also not explicitly represented; instead, we assume that upon remineralization, ammonia is immediately converted to nitrate. Importantly, these representations of denitrification and nitrification do not account for their impacts on alkalinity. For example, denitrification should add alkalinity to the water column, while nitrification should remove alkalinity (Paulmier et al., 2009), but neither processes is accounted for in our model. Inclusion of a more realistic nitrogen cycle that includes oxygen inhibition for nitrogen fixation and an explicit representation of denitrification and nitrification may improve simulated nitrate values and allow for more realistic limitation regimes, improving DIC and alkalinity, and hence the air-sea exchange CO 2 .

Changes in Drivers of the CO 2 Between CMIP5 and CMIP6
Table 4 also shows the dominant drivers of the CO 2 flux anomalies for E2-R. For most regions, the dominant drivers of the seasonal cycle of the CO 2 flux are the same between model versions. Exceptions are the subpolar North Atlantic in boreal fall, the equatorial Atlantic in boreal spring and fall, the subpolar North Pacific in boreal winter, spring, and summer, and the equatorial Pacific in boreal summer. For the subpolar North Atlantic, although both E2-R and observations show ΔpCO 2 to be the dominant driver of the monthly CO 2 flux anomalies in boreal fall, the ΔpCO 2 -driven anomalies in E2-R are of opposite sign to those of observations. In contrast, the ΔpCO 2 -driven anomalies in E2.1-G are of the same sign as those in observations, although are of smaller magnitude (Figures 8a and S19a). The increase in the ΔpCO 2 -driven anomalies in E2.1-G in boreal fall is due largely to an increase in the DIC-driven pCO 2,sw anomalies during this season (Figures 7a and S18a). While the seasonal amplitude of DIC in observations is better represented by E2-R than E2.1-G, the fall increase in DIC in E2-R lags the increase in observations by ∼1 month. Thus, when DIC-driven pCO 2,sw anomalies are averaged in boreal fall, they are negative in E2-R and positive in observations and (less so) in E2.1-G. In the equatorial Atlantic, ΔpCO 2 dominates the boreal fall and spring anomalies in the CO 2 flux in both models and observations. However, in boreal spring the pCO 2,sw anomalies are driven by DIC in E2-R and temperature in E2.1-G and observations. This appears to be caused by the boreal spring temperature maximum being higher in E2.1-G and observations than in E2-R, thus driving a higher temperature-driven pCO 2,sw maximum (Figures 7c and S18c). In contrast, in boreal fall DIC drives the pCO 2,sw anomalies in E2-R and observations, while for E2.1-G temperature continues to be the driver behind the pCO 2,sw anomalies. While the temperature driven pCO 2,sw anomalies are still better reproduced by E2.1-G than E2-R in boreal fall, the DIC-driven anomalies are better reproduced by E2-R, and in both E2-R and observations they are more negative than in E2.1-G. Notice also that alkalinity is also an important contributor to the seasonal cycle of pCO 2,sw in the equatorial Atlantic in boreal fall, being of similar magnitude to DIC, and the alkalinity-driven pCO 2,sw anomalies are better reproduced by E2-R than E2.1-G.
The subpolar North Pacific shows the largest number of changes in drivers of the CO 2 flux seasonal cycle between E2-R and E2.1-G, as these have changed in all seasons but boreal fall. Moreover, the dominant drivers in E2-R and observations are the same in all seasons, while they are only the same in boreal fall for E2.1-G. In the boreal winter and spring, wind speed is the dominant driver of the CO 2 flux in E2.1-G, and ΔpCO 2 is the dominant driver in E2-R and observations. This is due to a reduction in the seasonal amplitude of DIC in E2.1-G from E2-R (Figures 7e, S8e, and S18e), reducing its impact on the pCO 2,sw anomalies in E2.1-G. Note also that the effect of temperature on pCO 2,sw seasonality is underestimated in E2.1-G (with a comparable bias to that of DIC), and comparable to observations in E2-R (Figures 7e and S18e). Finally, in the equatorial Pacific, ΔpCO 2 is the dominant driver of the CO 2 flux anomalies in boreal summer in E2.1-G, while wind speed is the dominant driver in observations and E2-R. While the wind speed effect on the CO 2 flux anomalies is larger in observations that in E2-R during boreal summer, the discrepancy is increased in E2.1-G to the point that the anomaly is opposite in sign (Figures 8g and S19g). The decrease in the wind-driven CO 2 flux anomalies are largely driven by the reduction in both the seasonal amplitude in wind speed and annual mean pCO 2,sw (Figures 6b and S8) in this region from E2-R to E2.1-G.
While in most regions and periods the drivers of the CO 2 flux have remained the same between E2-R and E2.1-G, in the regions and seasons in which they are different, E2-R generally shows a better match to observations. Changes in the model from E2.1-G to E2-R that altered the drivers remain a subject of future investigation. Some candidate model changes that may have contributed to differences between E2-R and E2.1-G drivers of the CO 2 flux seasonality include: (i) the implementation of exponential detrital sinking profiles, which may have changed the fraction of fixed carbon exported below the euphotic zone and hence surface DIC and nutrient concentrations, (ii) the inclusion of prognostic alkalinity, which appears particularly important in controlling pCO 2,sw in the equatorial Atlantic, and (iii) changes in the parameterizations of vertical mixing and mesoscale eddies (Kelley et al., 2020). The latter may have impacted vertical and lateral transport of heat and DIC, contributing to the reduction temperature-driven and DIC-driven pCO 2 anomalies in the subpolar North Atlantic and Pacific.

Importance of pCO 2,sw in Driving Seasonality Changes
A somewhat expected result is that, between the 1% and the historical simulation in E2.1-G, changes in the monthly anomalies in the ΔpCO 2 gradient drive changes in the CO 2 flux in nearly all regions ( Figure 11). On the other hand, long-term changes in pCO 2,sw anomalies can be roughly divided by regions where the change is controlled by temperature and by DIC (Figure 10). Our findings are consistent with the currently observed changes in the seasonal cycle of pCO 2,sw , which show temperature dominated changes in the subtropical gyres and DIC + alkalinity driven changes in the subpolar regions and Southern Ocean (Landschützer et al., 2018).
We further find that it is the change in the sensitivity of pCO 2,sw to temperature or DIC (i.e., the change in the partial derivatives of Equation 2), as opposed to the change in temperature or DIC anomalies, that drives the change in the pCO 2,sw anomalies. This finding can be explained if we consider simple decompositions of the partial derivatives with respect to temperature (Takahashi et al., 1993) and DIC (Gallego et al., 2018): pCO DIC , respectively, to pCO 2,sw . γ T is a factor that has been shown experimentally to be independent of temperature (Takahashi et al., 1993). Hence, the increase in the sensitivity of pCO 2,sw to temperature between the end of the historical and the 1% simulation at the time of doubling of atmospheric CO 2 is expected, given that increasing atmospheric CO 2 also increases pCO 2,sw and inflates the right-hand side of Equation 6a.
On the other hand, γ DIC is a measure of the ocean's capacity to buffer changes in pCO 2,sw due to the accumulation of atmospheric CO 2 (Egleston et al., 2010;Hauck & Völker, 2015). γ DIC is related to the Revelle factor (R), the ratio of the relative change in pCO 2,sw to the relative change in DIC    where DIC is the annually averaged surface ocean DIC concentration. To determine the contributions of γ DIC and pCO 2,sw to the long term change of   2,sw pCO DIC requires an expression for R that can be evaluated at the end of the historical period and at the time of doubling of atmospheric CO 2 . In seawater, the Revelle factor may be approximated as (Sarmiento & Gruber, 2006)  (Figure 12). The differences between the two contributions are largest in the equatorial and subtropical regions, where the contribution of pCO 2,sw is much larger than that of γ DIC , and are very small in the high latitude regions, specifically the Southern Ocean and subpolar North Pacific. We also find that γ DIC decreases in each region, which in turn contributes to an increase   2,sw pCO DIC following Equation 6b. In addition, the decrease in γ DIC must be due to an increase in R, since the annually averaged DIC increases in each region ( Figure S7g), and increasing DIC increases γ DIC . Other studies have found that future changes in γ DIC impact the sensitivity of pCO 2,sw to DIC (Fassbender et al., 2018;Hauck & Völker, 2015). Using the MITgcm coupled to REcoM-2, Hauck and Völker (2015) found that between 2011 and 2,100, the decrease in γ DIC significantly contributed to the increase in the seasonal drawdown of pCO 2,sw in the Southern Ocean, with the change in biological productivity playing a relatively minor role. Fassbender et al. (2018), while not explicitly estimating γ DIC , used the ESM2M run under the RCP8.5 concentration pathway to investigate changes in the Revelle factor due to anthropogenic carbon invasion. They found that at sites in the Kuroshio extension, the subtropical North Atlantic, Irminger Sea, and south of the Antarctic Polar Front Ocean, R increased between 1,860 and 2,100, with the largest increases occurring at high latitudes (the Irminger Sea and Southern Ocean). They also emphasized that the effect of this increase in R is to magnify the sensitivity of pCO 2 to DIC. The result obtained in this study, that a decrease in γ DIC (increase in R) results in an increase in   2,sw pCO DIC , with the largest γ DIC -driven magnifications being in the high-latitudes (specifically the Southern Ocean and subpolar North Pacific; Figure 12), is consistent with these previous studies.
When examining changes in the seasonality of pCO 2 in seven globally-coupled CMIP5 models between averaging periods 2006-2026 and 2080-2,100, Gallego et al. (2018) found that the effects of changes in pCO DIC from changes in γ T and γ DIC , respectively, were small compared to the contribution from the changes in the annual mean pCO 2,sw in most regions. Since Gallego et al. (2018) do not provide their estimates of the change in γ DIC , we cannot compare them to the changes found in our study. Qualitatively, their results are consistent with those in our study, since the effect of the change in γ DIC is smaller than that of the change in pCO 2,sw in each region. Overall, the drivers of the change in seasonality in our model are consistent with those in the CMIP5 ensemble. Landschützer et al. (2018) compared drivers of the change in the seasonality of the non-thermal component of pCO 2,sw averaged between 1985 and 1989 and between 2010 and 2014, based on integrating SOCATv4 pCO 2,sw measurements into a neural-network. In their analysis, they decomposed the change in the monthly anomalies in non-thermal pCO 2,sw into changes driven by the change in the Revelle Factor (R = γ DIC DIC) and driven by changes in pCO 2,sw . They found that the Revelle-factor driven changes were about one-half of the pCO 2,sw driven changes, and that the change in the annual mean pCO 2,sw explained the majority of the and pCO 2,sw (blue bars) to the long term changes in   2,sw pCO DIC (μatm/ μmol) in each region. Changes are calculated as differences between quantities averaged over the historical and future climate periods described in Section 2.3. change in the seasonality of pCO 2,sw in the ocean, with the exception of the temperate latitudes in the southern hemisphere. They also found that at high latitudes, particularly near 60°N/S, the Revelle-factor driven changes are comparable with the pCO 2,sw driven changes (see their Figure S8). Thus drivers of changes in the seasonality of pCO 2,sw obtained from the neural-network are generally consistent to those found here.
Two notes of caution regarding our findings of long-term changes in the drivers of the seasonal cycles of pCO 2,sw and the CO 2 flux should be considered. First, in regions where the biases for both CO 2 seasonal cycles are largest, such as the subpolar North Atlantic and equatorial Pacific, the response of the drivers to changes in atmospheric forcing in our model may under or overestimate the magnitude of the response in the real ocean. Second, it is unclear how robust the long term changes in drivers found in our model are. A rigorous measure of robustness (e.g., the internal variability of changes in the drivers) would require an ensemble of 1% simulations which are not available. Thus, while finding some consistency between our results and previous studies is encouraging, further investigation is required to constrain the uncertainty (both structural and due to internal variability) of the response of the seasonal cycle of the CO 2 flux and pCO 2,sw , as well as their drivers, to increases in atmospheric pCO 2 . However, given the magnitude of the increase in atmospheric pCO 2 in the 1% simulation, we believe that this internal variability is overpowered by our forcing signal.

Conclusion
In this study, we have examined the seasonal cycles of ΔpCO 2,sw and the CO 2 flux in nine different oceanic regions in the NASA-GISS modelE GCM (GISS-E2.1-G). We compared them to the seasonal cycles of an observation-based climatology as well as those in the NASA-GISS submission to CMIP5 (model E2-R), and identified the drivers of both the observed and modeled seasonal cycles using an analysis based on first-order Taylor Series expansions. We then examined the change in the seasonal cycle between two simulations: a historical simulation averaged between 1995 and 2014, and a simulation in which atmospheric CO 2 reaches twice its pre-industrial value after increasing by 1% per year since 1850. Finally, we identified the drivers of the change in the seasonal cycles, elucidating the mechanisms by which the drivers change the seasonal cycles of the CO 2 flux and pCO 2,sw between the historical and 1% simulations.
We found general improvement in the seasonal cycle of the CO 2 flux in E2.1-G compared to E2-R. However, when viewed at a regional level, changes in model skill in capturing the CO 2 flux seasonal cycle were much more heterogeneous. Five (three) out of the nine regions analyzed here showed a smaller (larger) bias in the seasonal amplitude, while four (three) of these regions showed a smaller (larger) bias in the timing of seasonal extrema. The seasonal cycle of pCO 2,sw shows only marginal changes when viewing global skill metrics (Table 3). Regionally, the same number (three) of regions showed improvement and deterioration in the seasonal amplitude of pCO 2,sw , while the timing of seasonal extrema for pCO 2,sw have improved (deteriorated) in 2 (1) regions. In addition, while the seasonal cycle of NPP shows a general improvement, the seasonal cycles of DIC, alkalinity, and macronutrients show a general reduction in model skill in E2.1-G compared to E2-R. For E2.1-G, we found that the model seasonal cycles of the CO 2 flux and pCO 2,sw showed similar phasing to the observed seasonal cycles, though with generally increased seasonal amplitudes, in the subtropical regions and Southern Ocean. In the subpolar and equatorial regions, the model was often out of phase with observations (e.g., CO 2 flux in the equatorial Pacific and ΔpCO 2 in the subpolar North Pacific), or did not exhibit the same number of extrema as observations (e.g., CO 2 flux in the subpolar North Atlantic). In most regions, we find that temperature and DIC play the dominant roles in driving pCO 2,sw in both model and observations. However, the effects of DIC and temperature often oppose each other, such that in a few regions (e.g., Southern Ocean), alkalinity determines the seasonal cycle of pCO 2,sw .
In these regions of greatest disagreement, a combination of differences in model and observed DIC, temperature, and alkalinity-driven pCO 2 anomalies was responsible for the differences between the model and observed seasonal cycle of pCO 2,sw . In the subpolar regions, the largest differences were in the DIC and alkalinity components of the seasonal cycle of pCO 2,sw , while in the equatorial regions, the differences in the temperature, DIC, and alkalinity components of the seasonal cycle of pCO 2,sw were of similar magnitude, but varied depending on the time of year. These differences reflect a combination of model biases in NPP, transport of DIC, alkalinity, and nutrients to the surface ocean, and are also perhaps influenced by the lack of adequate observational data coverage in winter compared to summer. For the equatorial Pacific, they may also be influenced by the lack of climatological data for DIC and alkalinity. Future improvements to the model, including a particulate inorganic carbon tracer that would remove the dependence of alkalinity on total productivity, as well as explicit representations of denitrification and nitrification, should decrease these biases in the seasonal cycles.
When considering changes between the historical simulation and the 1% simulation, we found that in all regions, the ocean becomes a stronger atmospheric CO 2 sink, and the seasonality of ΔpCO 2 and the CO 2 flux increases in all regions except the subpolar North Pacific. However, the drivers of this stronger sink were regionally dependent. The effects of temperature are most important in the subtropical regions, while the effects of DIC are most important in the subpolar North Atlantic and Southern Ocean. DIC, alkalinity, and temperature all have effects that are similar in magnitude on changes in ΔpCO 2 seasonality in the equatorial regions. However, a common thread in all regions is that changes in the seasonality of ΔpCO 2 are driven by changes in the sensitivity of pCO 2,sw to changes in DIC, alkalinity, and temperature, as opposed to changes in the variables themselves. This finding is consistent with previous studies that show that the sensitivity of pCO 2,sw to its drivers increase (i) as the seawater pCO 2,sw concentration increases, and (ii) as the buffer capacity of seawater decreases. These findings do not account for natural variability in the CO 2 flux, pCO 2,sw , or their drivers, since we only analyze a single historical and a single 1% simulation, and not an ensemble. Finally, we clarify that the aim of this paper is to document the advances and biases for only a single member of the CMIP6 ensemble. The documentation presented here should be useful for assessments of future GISS model development and for investigations across the climate modeling community seeking to interrogate the multi-model spread of the CO 2 flux and its seasonality in the CMIP6 ensemble. It should also be helpful for investigators attempting to understand processes that need to be better parameterized, and hence better constrained by observations, in order to reduce model bias and increase model skill in future projections of the CO 2 flux. To these ends, future studies should examine the response of the seasonality of surface ocean pCO 2,sw and the CO 2 flux to climate change across the suite of CMIP6 models.

Data Availability Statement
Computational resources were provided by NASA Center for Climate Simulation High-End Computing Program through the NASA Center for Climate Simulation (NCCS) at Goddard Space Flight Center. All data sets used in this study are available from the NASA GISS climate group and through the CMIP6 repository (https://esgf-node.llnl.gov/projects/cmip6/). All methods, scripts, and codes will be readily available through the GISS climate group repository. Access will be unrestricted at https://data.giss.nasa.gov/.