A framework to evaluate and elucidate the driving mechanisms of coastal sea surface pCO 2 seasonality using an ocean general circulation model (MOM6-COBALT)

The temporal variability of the sea surface partial pressure of CO 2 (pCO 2 ) and the underlying processes driving this variability are poorly understood in the coastal ocean. In this study, we tailor an existing method that quantifies the effects of thermal changes, biological activity, ocean circulation and fresh water fluxes to examine seasonal pCO 2 changes in highly-variable 15 coastal environments. We first use the Modular Ocean Model version 6 (MOM6) and biogeochemical module Carbon Ocean Biogeochemistry And Lower Trophics version 2 (COBALTv2) at a half degree resolution to simulate the coastal CO 2 dynamics and evaluate it against pCO 2 from the Surface Ocean CO 2 Atlas database (SOCAT) and from the continuous coastal pCO 2 product generated from SOCAT by a two-step neuronal network interpolation method (coastal-SOM-FFN, Laruelle et al., 2017). The MOM6-COBALT model not only reproduces the observed spatio-temporal variability in pCO 2 but also in sea 20 surface temperature, salinity, nutrients, in most coastal environments except in a few specific regions such as marginal seas. Based on this evaluation, we identify coastal regions of ‘high’ and ‘medium’ agreement between model and coastal-SOM-FFN where the drivers of coastal pCO 2 seasonal changes can be examined with reasonable confidence. Second, we apply our decomposition method in three contrasted coastal regions: an Eastern (East coast of the U.S) and a Western (the Californian Current) boundary current and a polar coastal region (the Norwegian Basin). Results show


Introduction
The ocean plays an important role in offsetting human-induced carbon dioxide (CO2) emissions associated with cement production and fossil fuel combustion (Friedlingstein et al., 2019). Globally, the ocean is a net sink that absorbs roughly one quarter of the anthropogenic CO2 emitted into the atmosphere (-2.5 ± 0.6 Petagram of carbon per year (Pg C yr -1 ) for the  2018 decade, Friedlingstein et al., 2019). The spatio-temporal variability of this oceanic CO2 uptake is relatively well constrained in the open ocean thanks to several method including sea surface CO2 data-derived interpolations (e.g., Landschützer et al., 2014;Rödenbeck et al., 2014Rödenbeck et al., , 2015Takahashi et al., 2002), models and atmospheric inversions (e.g., Gruber et al., 2009Gruber et al., , 2019Keeling and Manning, 2014;Manning and Keeling, 2006), but it is less constrained and understood in the coastal ocean. Nonetheless, in recent decades, significant progress have been made with regard to the quantification and 40 analysis of the spatial distribution of the coastal air-sea CO2 exchange (FCO2) globally and regionally (e.g., Borges et al., 2005;Cai, 2011;Chen et al., 2013;Laruelle et al., 2010, Roobaert et al., 2019. The FCO2 seasonal cycle was also recently analyzed in coastal regions worldwide by Roobaert et al. (2019). This study identified that at the annual timescale, the global coastal ocean acts as an atmospheric CO2 sink (-0.2 ± 0.02 Pg C yr -1 ) with a more intense CO2 uptake occurring in boreal summer because of the disproportionate contribution of high latitude coastal regions in the Northern Hemisphere which cover 45 25 % of the total coastal area and are characterized by an intense CO2 sink in summer. A more in-depth analysis also revealed that the majority of the coastal seasonal FCO2 variations stems from the air-sea gradient in partial pressure of CO2 (pCO2), although changes in wind speed and sea-ice cover can be significant regionally.
Several processes influence the seasonal variations of surface ocean pCO2 and thus, the seasonality in FCO2. These processes 50 include changes in sea surface temperature (SST) tied to air-sea heat fluxes and ocean circulation, changes in sea surface salinity (SSS) associated with evaporation, fresh water fluxes (from land, ice-melt, precipitation and evaporation) and ocean circulation, as well as variations in sea surface alkalinity (ALK) and dissolved inorganic carbon (DIC) tied to biological activity, fresh water fluxes and ocean circulation (Sarmiento and Gruber, 2006). In the open ocean, the respective influence of these processes on the pCO2 variability has been interpreted using changes in SST, SSS, ALK and DIC observed in-situ (e.g., 55 Landschützer et al., 2018;Takahashi et al., 1993) or based on global/regional ocean biogeochemical models relying on a mechanistic, quantitative description of the physical, chemical and biological processes controlling the ocean carbon cycle (e.g., Doney et al., 2009). These investigations reveal that changes in SST (i.e. the thermal effect) is the main driver of the seasonal pCO2 in tropical oceanic regions, while non-thermal components (change associated with DIC, ALK and SSS) dominate at mid-and high-latitude (poleward of 40° N and 40° S, e.g., Landschützer et al., 2018;Takahashi et al., 2002). 60 In the coastal ocean, the processes controlling the pCO2 seasonal dynamics was mostly investigated regionally (e.g., Arruda et al., 2015;Frankignoulle & Borges, 2001;Laruelle et al., 2014;Nakaoka et al., 2006;Shadwick et al., 2010Shadwick et al., , 2011Signorini et al., 2013;Turi et al., 2014;Yasunaka et al., 2016) and only a few observation-based studies attempted to analyze the coastal pCO2 seasonal variability into processes at the global scale (Cao et al., 2020;Chen and Hu, 2019;Laruelle et al., 2017). 65 Regional studies using either observations or model results have covered, e.g., the shelves of the entire Atlantic basin (Laruelle et al., 2014), the West (California Current, Turi et al., 2014) and East (e.g., Shadwick et al., 2010Shadwick et al., , 2011Signorini et al., 2013) coasts of the United States, as well as the South and Southeast Brazilian shelves, Uruguayan and Patagonia shelves and shelves along the southwestern Atlantic ocean (Arruda et al., 2015). In the California Current, the strong upwelling of carbon-rich waters was identified as the main control of the pCO2 seasonality (Turi et al. 2014). On the Patagonia shelf, the thermal effect 70 and biological pumps were found to be the main drivers of the seasonal pCO2 variability with only a small contribution from the ocean circulation (Arruda et al., 2015), while along the East coast of the U.S, seasonal thermal changes play the major role (Shadwick et al., 2010(Shadwick et al., , 2011Laruelle et al., 2015;Signorini et al., 2013). These studies are, however, confined to specific regions and a global picture of the mechanisms driving the coastal pCO2 dynamics is still missing. In addition, the attribution analysis into specific physical and biological processes is incomplete. Indeed, the attribution relies on a linear decomposition 75 linking variations in sea surface ocean pCO2 to seasonal changes in DIC, ALK, SST and SSS (e.g., Signorini et al., 2013Lovenduski et al., 2007;Takahashi et al., 1993;Turi et al., 2014), or on a series of sequential simulations isolating biological and physical terms therefore ignoring how covariations between the different terms dampen or reinforce each other (e.g., Arruda et al., 2015;Turi et al., 2014).

80
In this study, we develop a new framework to elucidate the seasonal pCO2 dynamics of the global coastal ocean. This framework relies on the global Modular Ocean Model version 6 (MOM6, Adcroft et al., 2019) from the NOAA Geophysical Fluid Dynamics Laboratory coupled to the biogeochemical module Carbon Ocean Biogeochemistry And Lower Trophics version 2 (COBALTv2, Stock et al., 2014Stock et al., , 2020. MOM6-COBALT model outputs provide the relevant variables and processes that are required to perform an explicit decomposition of the inorganic carbon dynamics (Liao et al., 2020) in the 85 entire coastal domain. These outputs are then analyzed using a novel approach to attribute seasonal variations in surface ocean pCO2 to changes in biological activity, ocean circulation, SST, air-sea CO2 fluxes and fresh water fluxes (Liao et al., 2020), and which is here enhanced for the coastal ocean. The decomposition method constitutes a significant improvement upon previous studies. First, it accounts for co-variations in biological and physical processes and how their evolution jointly modulates the pCO2 signal. Second, it improves on the traditional linear approaches developed for the open ocean (Sarmiento 90 and Gruber, 2006;Takahashi et al., 1993) and used since then (e.g. Lovenduski et al., 2007), because, as shown later in this study, the linear decomposition introducing significant biases in coastal waters due to the larger range in DIC, ALK, pH and salinity values encountered in the variable coastal environment (Egleston et al., 2010).
the coastal regions where the model best reproduces the observed ocean pCO2 variability and can thus be considered most suitable for a detailed analysis of the drivers of the pCO2 seasonal changes.
-Second, to illustrate the capabilities of our upgraded decomposition framework, we examine the drivers of the pCO2 100 seasonality in three contrasted coastal regions: The East coast of the U.S, the West coast of North America and the Norwegian Basin.

Ocean biogeochemical model description
In this study, we used the ocean model MOM6 and the Sea Ice Simulator version 2 (fourth generation of ocean-ice models 105 called OM4) detailed in Adcroft et al. (2019). The version of OM4 adopted here is OM4p5 which has a nominal horizontal resolution of 0.5° (i.e. with a finer latitudinal resolution of 0.26° in the tropical region). On the vertical, it includes 75 hybrid coordinates with a z* coordinate near the surface (geopotential coordinate allowing free surface undulations) and a modified potential density coordinate below. The vertical spacing increases from 2 m in the upper 20 m (i.e first 10 layers) to larger isopycnal layers below. Layers in z* broadly deepens towards high latitudes (see Adcroft et al., 2019 for details on the grid). 110 This ocean-ice model is coupled to the biogeochemical module COBALT version 2 (COBALTv2), which includes 33 state variables to resolve global-scale cycles of carbon, nitrogen, phosphate, silicate, iron, calcium carbonate, oxygen and lithogenic materials (Stock et al., 2020). Details about the planktonic food web dynamics in COBALT, and global assessments of largescale carbon fluxes through the food web such as net primary production can be found in Stock et al. (2014Stock et al. ( , 2020. The ocean model is forced by the 55-km horizontal resolution Japanese atmospheric reanalysis (JRA55-do) version 1.3 at a 3-hour 115 frequency between 1959 and 2018 (Tsujino et al., 2018), and the atmospheric CO2 concentration data (xCO2) from the Earth System Research Laboratory (Conway et al., 1994;GLOBALVIEW-CO2, 2004). The xCO2 is converted to pCO2 using atmospheric and water vapor pressures by the model. SST, SSS, sea surface nutrients (nitrate, phosphate, silicate) and oxygen were initialized from the World Ocean Atlas version 2013 (Garcia et al., 2013a(Garcia et al., , 2013b. Initial DIC and ALK conditions are taken from GLODAPv2 (Olsen et al., 2016). The initial DIC is corrected for the 120 accumulation of anthropogenic carbon to match the level expected in the first year of simulation (1959) using the data-based estimate of ocean anthropogenic carbon content of Khatiwala et al. (2013). At the end of a 81-year spin-up repeating year 1959, the model has reached a near-equilibrium between atmospheric pCO2 and surface ocean pCO2, with a drift in global airsea CO2 flux < 0.004 Pg C yr -1 over the last 10 years of spin-up. Further details on the configuration, spin-up and simulation can be found in Liao et al. (2020). 125

Observational products and model evaluation
We first evaluate the ability of MOM6-COBALT to reproduce the observed spatial distribution of environmental variables in the coastal domain, namely the SST, SSS and sea surface nutrients (nitrate, phosphate and silicate). The observational SST and SSS fields are from the daily NOAA OI SST V2 (Reynolds et al., 2007) and the daily Hadley center EN4 SSS (Good et al., 2013), respectively. The observed nutrient fields in the sea surface are extracted from the World Ocean Atlas version 2018 130 (Garcia et al., 2019). We also compare the simulated coastal pCO2 directly to "raw", un-interpolated observations extracted from the Surface Ocean CO2 Atlas database (SOCAT), using monthly observations from SOCAT version 6 gridded at the spatial resolution of 0.25 degree (SOCATv6, Bakker et al., 2016). For the evaluation period used in this study (1998 -2015), this database contains 9.8 million pCO2 observations within the coastal domain. All data from SOCATv6 are converted from fugacity of CO2 in water to pCO2 using the formulation of Takahashi et al. (2012). We finally compare the pCO2 simulated by 135 the MOM6-COBALT model to the 0.25° continuous monthly pCO2 fields generated from the SOCAT observations by the two-step neuronal network (SOM-FFN) in coastal regions (Laruelle et al., 2017). The SOM-FFN data product of Laruelle et al. (2017) is thus not "raw" and implies a significant amount of statistical modelling. It is also derived from an earlier version of SOCAT (SOCATv4, Laruelle et al., 2017) than the "raw" one. In what follows, the pCO2 products generated by the model, the statistical interpolation of observations, and the un-interpolated observations will be referred to as MOM6-COBALT, 140 coastal-SOM-FFN and SOCATv6, respectively. All observational and simulated fields are converted from their original spatiotemporal resolution to monthly 0.25° gridded climatologies for the 1998 -2015 period to match the one used by the coastal-SOM-FFN. Cells that are covered by more than 95 % of sea-ice are removed from the comparison since we assume no transfer of our master variable (pCO2) through sea ice. In our analysis, we apply the broad definition of the coastal zone by Laruelle et al. (2017), using a global mask that excludes estuaries and inland water bodies while its outer limit is set 300 km away from 145 the shoreline. This definition leads to a total surface area of 77 million km² which is split into 45 coastal regions using the MARgins and CATchment Segmentation (MARCATS, Laruelle et al., 2013). These 45 regions are grouped into 7 broad classes with similar hydrological and climatic settings (Liu et al., 2010): (1) Eastern and (2) Western Boundary Currents (EBC and WBC respectively), (3) tropical margins, (4) subpolar and (5) polar margins, (6) marginal seas and (7) Indian margins.

150
The model evaluation of all gridded environmental variables including pCO2 is performed for the annual mean and the seasonal cycle both globally and within each of the 45 MARCATS regions. For the seasonal analysis, for each variable, a climatological monthly anomaly is calculated as the difference between the variable x for a given month and its climatological annual mean.
The evaluation of the seasonal amplitude is then performed using the bias between observed and simulated root mean square (RMS) of their monthly anomalies. A positive bias represents a larger simulated seasonal amplitude than derived from the 155 observations. The temporal shift between observed and simulated seasonal cycles is also assessed from the Pearson correlation coefficient (no units) of the regression between monthly times series simulated by MOM6-COBALT and those extracted from the observations. These comparisons not only serve to assess the overall model's performance in reproducing observations but also help identifying potential discrepancies between observed and simulated environmental fields (e.g., SST, SSS) that are used by the two-step neuronal network coastal-SOM-FFN to generate the continuous pCO2 climatology. We use two metrics 160 to evaluate SOCATv6 spatial and temporal coverage. First, we evaluate the spatial coverage at the MARCATS scale by computing the percent surface area sampled by SOCATv6 data for each MARCATS. A 50 % spatial coverage means that SOCATv6 data are available in 50 % of the 0.25° x 0.25° cells included in this specific MARCATS (this metric is used in Fig.   1a). Second, we evaluate the ability of SOCATv6 to capture the seasonality at the grid cell scale by computing the number of months where at least one SOCATv6 pCO2 measurement for each 0.25° x 0.25° grid cells. A 8-months temporal coverage 165 means that 8 out of the 12 months are sampled at least once in this grid cell (this metric is used in Fig. 6a).
Finally, from this global and regional spatio-temporal evaluation, we label the model to coastal-SOM-FFN agreement ('high', 'medium' and 'low') for each MARCATS and identify regions for which our results are the most robust for further in-depth analysis of the processes driving the coastal pCO2 dynamics. The labels of agreement are based on 3 criteria: First, we assess 170 whether the simulated annual mean pCO2 is within 20 μatm of the one extracted from the coastal-SOM-FFN. This threshold of 20 µatm roughly corresponds to the globally averaged pCO2 gradient between the atmosphere and the coastal sea surface (Laruelle et al., 2018). The second and third criteria evaluate the magnitude and phasing of the simulated pCO2 seasonal cycle against the coastal-SOM-FFN, using an absolute bias in the seasonal magnitude < 20 μatm and a Pearson coefficient > 0.5 as threshold. The agreement is considered 'high' when the 3 criteria are fulfilled, 'medium' when criteria 2 and 3 are satisfied 175 and 'low' when only one or zero criteria is met on the seasonality.

Processes controlling seasonal pCO2 variability: a method tailored for coastal regions
pCO2 in surface sea water can be computed from DIC and ALK following Eq. (1) (Sarmiento and Gruber, 2006;Wolf-Gladrow et al., 2007): where K 0 ′ is the aqueous-phase solubility constant of CO2 in water and K 1 ′ and K 2 ′ represent the apparent equilibrium dissociation constants of the carbonate system. Several physical and biogeochemical processes can thus affect pCO2 via changes in DIC, ALK and/or via the K 2 ′ K 0 ′ K 1 ′ term which depends on SST and SSS. To quantify the processes controlling the pCO2 185 variability at the seasonal timescale of interest to this study, we adopt the method of Liao et al. (2020). The method starts from the traditional approach that links variations in sea surface ocean pCO2 to changes in DIC, ALK, SST and SSS using the following linear decomposition Lovenduski et al., 2007;Takahashi et al., 1993;Turi et al., 2014): Where the "∆x" terms represent the seasonal anomaly of x (i.e. the departure from the annual mean) and are coefficients that describe the sensitivity of pCO2 to changes in DIC, ALK, SST and SSS. The coefficients for DIC, SST and SSS are always positive as pCO2 increases with increases in DIC, SST or SSS, while the coefficient for ALK is always negative as pCO2 systematically decreases with increasing ALK. These coefficients are generally estimated using the 195 approach of Sarmiento and Gruber (2006) (see Eq. S1-S4 in Appendix), which has been widely used in the open ocean (Liao et al., 2020;Sarmiento and Gruber, 2006;Takahashi et al., 1993). In this study, we refine the estimation of the coefficients so they can be used for the wide range of DIC/ALK ratios that can be encountered in the coastal waters. This includes conditions when the DIC/ALK ratio is close to 1, such as in regions with significant freshwater discharge like those found near estuarine mouths or on polar shelves subject to sea-ice melting, when pH is around 7.5 (Egleston et al., 2010). In these cases, the 200 traditional approximation method using mean DIC, ALK, SSS and SST fields breaks down (see Eq. (S1-S2) and Figure S1 in the Appendix). To circumvent this important limitation, we computed the coefficients of the pCO2 dependency using a regression approach based on the CO2SYS program (Lewis and Wallace, 1998). At each point in space, pCO2 was computed using the 1998 -2015 average of DIC, ALK, SSS and SST with CO2SYS (method 14 in CO2SYS Matlab program, Millero, 2010). The ∂pCO 2 ∂DIC coefficient was then computed as the slope of the linear regression between pCO2 and DIC obtained by 205 allowing DIC to vary around the local mean DIC value while keeping other tracers (ALK, SST, SSS) constant. The DIC range used to compute the slope was set to the ± 2 standard deviation of the 1998-2015 monthly values at that location with an upper bound at ± 60 µmol kg -1 (see Appendix for further details). The same approach was repeated to compute the coefficients for the pCO2 dependence on ALK, SST and SSS, respectively. Our methodology leads to coefficients that are constant in time but space dependent. In Fig. S1, we compare the coastal pCO2 reconstructed from the traditional decomposition (using the space 210 varying coefficients reported by Sarmiento and Gruber, 2006) with those computed here using the CO2SYS regression. For the global coastal ocean, we find a large bias (global mean rmse of fitting pCO2 anomaly in Eq. (2) = 14.6 µatm), which is especially pronounced at high latitudes. In contrast, the decomposition method based on our methodology reduce drastically the biases (global mean rmse = 2.8 µatm) in coastal regions and allows a more robust reconstruction of the pCO2 variability.

215
We further evaluated how using coefficients that are both time and space varying could reduce the residual biases between our pCO2 decomposition (using space dependent coefficients that are constant in time) and the pCO2 simulated in the model that are found in regions with large freshwater discharge, such as the mouth of the Amazon River or Arctic coastal waters. We compare the pCO2 seasonality simulated by the model to the pCO2 reconstructed by the three methods (space varying coefficients from Sarmiento and Gruber (2006); regression-based space varying coefficients; and regression-based space and 220 time varying coefficients) using a point in the Amazon River plume (points at 310.25°E -1°N, Fig. S1d and S1e). At this location, the use of the regression-based coefficients greatly improves the recovery of the simulated pCO2 compared to using the traditional coefficients of Sarmiento and Gruber (2006), reducing the rmse from 83 µatm to 24 µatm. The use of both space and time dependent regression-based coefficients further reduces this bias, bringing down the rmse from 24 uatm to 18 uatm.
This additional improvement is however marginal, motivating our choice to use the simpler approach of the space dependent 225 only coefficients.
Here we assume that the coefficients are constant in time, and the temporal change in pCO2 (∂ t pCO 2 in µatm month -1 ) can therefore be expressed as a simple function of the temporal changes in DIC (∂ t DIC), ALK (∂ t ALK), SST (∂ t SST) and SSS Temporal changes in DIC, ALK, SST, and SSS (∂ t DIC, ∂ t ALK, ∂ t SST and ∂ t SSS) are controlled by surface heat flux, ocean transport, freshwater fluxes, biological processes, and the air-sea CO2 flux. Using the model results, we further expand the 235 decomposition to quantify the contribution of these physical and biological processes (see details of derivation in Liao et al, where the temporal changes in pCO2 (time tendency called pCO2 change) is on the left-hand side (LHS), and the five terms that control this change in pCO2 are on the right-hand side (RHS) of the equation. Subscripts h and v denote the contribution from horizontal (advection and diffusivity in the meridional and zonal directions) and vertical (vertical advection and diffusivity) transports on SST, SSS, DIC and ALK, bio denotes the DIC and ALK changes induced by biological processes (photosynthesis, respiration, and calcium carbonate dissolution/precipitation, denitrification and nitrification), q denotes the 250 effect of surface heat flux on SST, fw denotes the effect of fresh water fluxes (i.e., precipitation, evaporation, river runoff and sea-ice formation and melting) on SSS, DIC and ALK, and the term CO2 flux denotes the DIC change induced by air-sea CO2 exchange.
Here we examine changes in pCO2 attributed to three oceanic processes that modify the concentration in dissolved species (i.e. 255 DIC, ALK and SSS), namely their transport by oceanic circulation (circ, which include horizontal and vertical transport), the effect of dilution/concentration due to freshwater fluxes (fw) and the effect of biological activity (bio), and isolate the thermal influence tied to SST changes induced by both oceanic transport and air-sea exchange of heat. Finally, the air-sea CO2 exchange (CO2 flux) pushes the surface pCO2 concertation towards its equilibrium with the atmosphere and systematically acts to offset the pCO2 changes associated with the sum of the internal oceanic processes (circ, bio, fw and thermal). In this study, we apply 260 Eq. (4) using averages between the sea surface and the mixed layer depth (MLD), defined here as the depth where the water density is 0.01 kg m -3 denser than the water at the surface (minimum MLD is 5 meters). Positive contributions on the RHS would yield an increase in pCO2 (positive pCO2 response on the LHS). Positive values of the CO2 flux correspond to an ocean CO2 uptake. This method to decompose the pCO2 seasonality into controlling processes in the coastal domain is illustrated in three coastal regions: The East and West coast of North America and in the Norwegian Basin. 265 be analyzed in detail in Section 3.1.3, but before we do so we first perform a data-model evaluation according to the following:

Annual mean state and seasonal cycle model evaluation and identification of coastal regions
We first evaluate the model by comparing simulated fields of SSS, SST, sea surface nutrients to observations globally and regionally (Sect. 3.1.1, Figs. 2 and 3). Second, the ability of the model to capture the coastal pCO2 annual mean and seasonality is assessed against the "raw" SOCATv6 data and the continuous monthly observational-based pCO2 product (coastal-SOM-275 FFN, Laruelle et al., 2017), respectively (Sect. 3.1.2, Figs. 3-6).

Model evaluation for coastal waters environmental variables
MOM6-COBALT captures fairly well the main spatial patterns of key environmental parameters (SST, SSS and sea surface nutrients) in the coastal domain (Fig. 2). The global SST field simulated by the model reproduces the strong large-scale tropical to polar SST gradients, with a global median bias of -0.2 °C (Fig. 2a-c), and biases at the scale of MARCATS regions ranging 280 from 0 °C in the North East Atlantic (M17) to 1.3 °C in the East coast of the U.S (M10, Fig. 3a and Table S1). With a global median bias value of 0.2, the model also correctly reproduces the observed SSS patterns which are mainly regulated by evaporation and freshwater inputs from precipitation, riverine runoff and ice melt, with lower SSS values in polar regions and along the coasts in Southeast Asia and higher SSS values along the coasts of evaporation basins such as in the Arabian or the Mediterranean Sea (Fig. 2d-f). The SSS analysis at the MARCATS scale reveals absolute SSS biases generally less than or 285 close to 1 except for five MARCATS where absolute biases exceed 2. These MARCATS are mainly located in marginal seas (the Baltic Sea, M18, the Black Sea, M21 and the Persian Gulf, M29), but also include one polar region (the Canadian Archipelago, M13) and one tropical region (Tropical West Atlantic, M7, Fig. 3b and Table S1). Similar to SSS, largest modeldata discrepancies for nutrients are mostly found in marginal seas (Fig. 3c-e and Table S1). For instance, the largest PO4 and SiO4 biases are encountered in the Black Sea (M21, absolute biases of 3 and 75 µmol kg -1 , respectively). The Peruvian 290 upwelling (M4), the Bay of Bengal (M31) and the N-E Pacific (M1) also present large biases in NO3 and PO4, respectively (e.g., NO3 bias of 8 µmol kg -1 for M4). The global median nutrients biases are however much smaller, reaching 0.3, -0.2 and -0.4 µmol kg -1 for nitrate (NO3, Fig. 2i), phosphate (PO4, Fig. 2l) and silicate (SiO4, Fig. 2o), respectively, The model-data seasonal evaluation reveals that MOM6-COBALT reproduces the global SST and SSS amplitudes remarkably well (median absolute bias of 0.1 °C and 0.0, respectively, Table S2

Model evaluation for coastal pCO2
The spatial distribution of the annual mean pCO2 simulated by MOM6-COBALT is in good agreement with the observational pCO2 values extracted from the "raw" SOCATv6 database with generally low pCO2 values (blue colors) in temperate and high latitudes and high pCO2 values (yellow and red colors) in tropical and sub-tropical regions (Figs. 4a-c). The model-data pCO2 310 evaluation at the regional scale shows that 33 of the 45 MARCATS present absolute biases lower than 20 µatm (Table S1). pCO2 product which uses a neural network interpolation method to fill data gaps and resolve the spatio-temporal coastal pCO2 variability globally.
Our results show that MOM6-COBALT reproduces the main spatial features of the annual mean pCO2 field captured by the coastal-SOM-FFN product, as revealed by the relatively low globally averaged bias of 2.5 µatm (Figs. 4a and 4d). In both the 320 model and the SOM-FFN product, low coastal pCO2 values are consistently found in temperate and high latitude regions in both hemispheres, while high pCO2 values are largely limited to (sub)tropical regions. Largest discrepancies (Fig. 4e) are found at high latitudes (poleward of 60° N and 60° S, negative bias), along the Eastern Boundary Peruvian and Namibian upwelling systems (high positive bias) and more locally close to the mouth of some large rivers (e.g., the plume of the Amazon or the Rio de la Plata, high negative bias). We note however that these regions are poorly sampled in the SOCATv6 dataset 325 (Fig. 4b) and are thus likely weakly constrained in the coastal-SOM-FFN product (Fig. 4d). At the regional scale, differences in annual mean pCO2 between MOM6-COBALT and coastal-SOM-FFN are lower than 20 µatm in 35 MARCATS (Table S1, Fig. 3f), which partly is a reflection of the low annual mean biases observed in the environmental driver variables in these regions (see Sect. 3.1.1). In EBC, WBC, and subpolar coastal regions, the model tends to overestimate the regional mean pCO2 compared to coastal-SOM-FFN (positive bias), except along the East coast of U.S 330 (M10), in the China and Kuroshio seas (M39) and in the North East Atlantic (M17 , Table S1). In polar regions, the model generally underestimates the mean pCO2 compared to coastal-SOM-FFN, except around the South of Greenland (M15). In Indian, marginal, and tropical coastal regions, no general trend can be identified regarding the sign of the bias, which can be positive or negative.
Quantitatively, the 10 MARCATS with absolute biases > 20 µatm are mainly located in regions for which very limited or no 335 observational data have been compiled in the SOCATv6 database (Table S1)  Our analysis reveals that the seasonal amplitudes simulated by MOM6-COBALT are systematically larger than the ones estimated by the coastal-SOM-FFN product ( Fig. 5a-b, red colors in Fig. 5c and positive biases in Table S2) for all coastal regions belonging to EBC, WBC, Indian and tropical margins. For the majority of polar and subpolar margins and for some 345 marginal seas, the model simulates lower seasonal pCO2 amplitudes (blue colors in Fig. 5c and negative biases in Table S2).
Quantitatively, absolute biases between the modelled and coastal-SOM-FFN amplitudes do not exceed 20 µatm except for marginal seas where larger discrepancies are calculated (6 of the 9 marginal MARCATS, Table S2). The monthly mean pCO2 seasonal cycle simulated by MOM6-COBALT is also well in phase (Pearson correlation coefficients > 0.5) with the one extracted from coastal-SOM-FFN in 34 out of the 45 MARCATS (red colors in Fig. 5d and Table S2). The agreement is 350 especially good in the best monitored MARCATS regions (MARCATS where > 50 % of the area is covered by SOCATv6 observations, Table S1). For instance, in regions with good data coverage such as along the East coast of the U.S (M10), the Norwegian Basin (M16), the Californian Current (M2), the Leeuwin Current (M33), or the Brazilian Current (M6), the Pearson correlation coefficient is higher than 0.9 (Table S2). In contrast, the seasonal pCO2 cycle simulated by MOM6-COBALT substantially diverges from that of the coastal-SOM-FFN in four poorly monitored marginal seas and in a few of the EBCs, 355 Indian margins, subpolar margins, and tropical margins (Pearson correlation coefficient < 0.5, Table S2 and blue colors in Fig.   5d).
The model pCO2 seasonal evaluation against SOCATv6 is only performed in 11 MARCATS namely the Californian Current Zealand (M36). The modeled seasonal cycle is in good agreement with that one derived from SOCATv6 ( Fig. 6b-n, Table S2) with absolute biases < 20 µatm for all of the 11 selected MARCATS and Pearson correlation coefficients close to 0.5 or higher except for the Iberian Upwelling (M19, Pearson value of 0.2) and in the New Zealand shelf (M36, value of 0.3). We did not perform the SOCATv6-model seasonal evaluation for the other MARCATS because the vast majority of grid cells only include data for less than 4 climatological months (Fig. 6a). However, we also evaluated the simulated pCO2 seasonality against 365 SOCATv6 in regions where this evaluation is not possible to be performed at the MARCATS scale. To do so, we selected four sites of smaller spatial extent than MARCATS for which we calculated climatological seasonal pCO2 signals from the SOCATv6 dataset and compared them with the model pCO2. These sites are located off the Antarctic Peninsula, on the Queensland Plateau in NE Australia, in coastal waters of Papua New Guinea and of Terra Nova (see black boxes in Fig. 6a).
In those regions, the absolute biases on the seasonal amplitude between    (Table S1). Together, these 29 MARCTAS represent 65 % of the global coastal ocean surface area. For the 11 MARCATS that are best covered by observations (MARCATS where > 50 % of the surface area is covered 380 by SOCATv6 observations, Table S1), absolute biases for the annual mean are always < 20 µatm for the three product intercomparison, except in the Californian Current (M2), in the Baltic Sea (M18) and along the N-E Pacific (M1). The seasonal MOM6-COBALT against coastal-SOM-FFN evaluation also reveal that 39 of the 45 MARCATS have pCO2 seasonal amplitude biases < 20 µatm and 34 MARCATS have a Pearson correlation coefficient > 0.5 (Table S2).

385
Based on this evaluation, we attribute for each MARCATS a level of confidence on the model to coastal-SOM-FFN agreement ('high', 'medium' and 'low', Table 1 and Fig. 1a). Out of the 45 MARCATS, 25 are labeled with 'high' agreement, that is to say, they fulfil the following criteria regarding the annual mean and the seasonality (Table 1 and (Table S1). In addition, 7 'high' agreement MARCATS also show a data density > 50 % (13 MARCATS if we lower the data coverage to > 30 %, Fig. 1a (Table S1 and Table S2 except M22, M33 and M9 for the nutrient phasing) and are therefore excellent potential candidates for an analysis of the processes controlling the coastal pCO2 dynamics. 6 additional MARCATS regions fulfil the criteria related to the seasonal pCO2 evaluation while they fail to fulfil the annual 400 mean pCO2 bias threshold of 20 µatm. These 'medium' agreement regions (Table 1 and

Methodological limitations
While our results show a relatively good agreement between MOM6-COBALT and coastal-SOM-FFN regarding the spatial 410 and temporal pCO2 distribution over the global coastal ocean, the comparison remains challenging for several reasons.
First, while the climatology of Laruelle et al. (2017, coastal-SOM-FFN) is currently the best available product for a modeldata comparison, it has its own limitations. For instance, in some regions, particularly coastal upwellings such as the Moroccan (M22) and Peruvian (M4) upwellings, the pCO2 fields generated by the coastal-SOM-FFN do not reproduce well the high and 415 variable pCO2 values measured in-situ (see e.g., Friederich et al., 2008 andMcGregor et al., 2007). Such poor performance of the coastal-SOM-FFN algorithm in these types of systems were already identified by Laruelle et al. (2017). Indeed, upwelling regions are still relatively poorly monitored and expand partly beyond the coastal domain used by Laruelle et al. (2017), leading to locally skewed calibration of the SOM-FFN. Deficiencies in the observation-based product can thus partly explain the large model-data bias (106 µatm, largest of all MARCATS) calculated in the Peruvian upwelling region. Moreover, although the 420 Surface Ocean CO2 Atlas database (SOCAT) has expanded significantly over the past few years, some regions are still poorly monitored. In the coastal regions where no observational data exist (e.g., in the Black Sea, the Sea of Okhotsk, the Bay of Bengal, Fig. 4b) in the SOCAT database used here (SOCATv6, Bakker et al., 2016), it is difficult to evaluate the performance of the SOM-FFN and, thus, of an OGCM in reproducing the pCO2 field. In addition, for certain regions subjected to complex dynamic biogeochemical settings (e.g., upwelling, seasonal cover of sea-ice, influenced by rivers, marginal seas), the pCO2 425 field reconstructed by the SOM-FFN suffers from poor performance, which can partly be explained by the lack of observational data. This lack of observations could partly explain why MOM6-COBALT-coastal-SOM-FFN pCO2 biases exceed 20 µatm in these regions. The seasonal model evaluation against raw SOCATv6 is limited at the MARCATS scale and mainly performed against coastal-SOM-FFN due to the very few coastal regions that contain a continuous climatological seasonal pCO2 cycle (Fig. 6a) in the SOCATv6 database. This study highlights the regions (Fig.1a, e.g., Indian ocean margins, Peruvian 430 upwelling, marginal seas) where new observational data are most urgently needed, specifically collected during periods of the years that are currently not covered to improve our understanding of the CO2 exchange between coastal regions and the atmosphere at the regional and global scales. In addition, only one global continuous pCO2 climatology derived by the SOM-FFN method currently exists for the coastal ocean. It would therefore be beneficial for the community to develop other observation-based climatologies relying on other interpolation techniques, as currently the case for the open ocean. 435 Second, the model-data comparison should also be analyzed in the light of the current limitations in the model itself. OGCMs have been designed for global ocean applications and the coarse spatial resolution of these models, on the order of 0.5° in the present study, cannot resolve accurately mesoscale and sub-mesoscale processes as well as tidal mixing in shelf regions even with a model configuration including parameterizations for these processes. The coastal currents are also not always well 440 resolved because of the coarse resolution of the shelf bathymetry. These small-scale hydrodynamic features are known to affect the spatio-temporal variability of pCO2 and the air-sea CO2 exchange (Bourgeois et al., 2016;Kelley et al., 1971;Lachkar et al., 2007;Laruelle et al., 2010). Therefore, although MOM6-COBALT runs at 0.5°, discrepancies between coastal-SOM-FFN and MOM6-COBALT in narrow EBCs such as the Peruvian Upwelling Current (M4) and along South west Africa (M33) could also be explained by the limited spatial resolution of the model. Moreover, OGCMs such as MOM6-COBALT have a 445 relatively simple representation of the biogeochemistry which does not fully captures some of the important processes of the carbon dynamics in coastal waters such as sea-ice temporal dynamics (Adcroft et al., 2019), neritic calcification (O'Mara and Dunne, 2019), or terrestrial and marine organic matter decomposition and burial (Lacroix et al., 2021). Moreover, the largest biases observed in marginal seas can partly be explained by large fluvial inputs and oceanic water flows through fine scale topography (e.g. straits) that are poorly represented in global OGCMs. 450 Finally, the annual mean/seasonal pCO2 biases between the coastal-SOM-FFN and MOM6-COBALT can also be traced back to divergences in the environmental fields simulated by the model compared to observations (Table S1 and Table S2). For instance, in most marginal seas, the model poorly resolves the annual mean and seasonal cycle of SSS and nutrients compared to the observations. These discrepancies impact the simulated pCO2 via the controls of the SSS on the CO2 solubility and of 455 nutrients on the biological pump and CO2 uptake. In the tropical W Atlantic (M7) which is under the influence of the Amazon River, the model simulates lower annual mean SSS (and therefore lower pCO2) than the observations. In the tropical E Pacific (M3) and in South-East Asia (M38), the poor agreement between simulated and observed seasonal pCO2 cycle could be explained by significant biases in the nutrient seasonal cycles (low Pearson correlation coefficient). Interestingly however, some regions reveal significant biases in the major environmental fields but not in the pCO2 (e.g., Tropical W Atlantic, M7) 460 while in other regions, the reverse is observed (e.g., the Mediterranean (M20) and W Arabian (M27) Seas and in New Zealand (M36)). Also, for some regions biases in environmental fields do not affect the pCO2 as expected. For instance, along the East coast of the U.S (M10), MOM6-COBALT simulates larger SST compared to observations while the simulated pCO2 is lower compared to coastal-SOM-FFN on an annual mean. This clearly shows that biases in environmental fields are not sufficient to explain fully the biases in pCO2 diagnosed between MOM6-COBALT and coastal-SOM-FFN. 465

Processes governing the seasonal pCO2 variability
Our second objective is to examine the drivers of the pCO2 seasonality in three well sampled and contrasted coastal regions where the model to coastal-SOM-FFN agreement is satisfactory: The East coast of North America (M10), the Norwegian Basin (M16) and the Californian Current (M2). The East coast of North America is a sink of atmospheric CO2 that has been extensively studied over the past decade (e.g., Fennel et al., 2019;Laruelle et al., 2015;Shadwick et al., 2010Shadwick et al., , 2011Signorini 470 et al., 2013). The pCO2 spatio-temporal dynamics in this MARCATS is particularly well captured by MOM6-COBALT ('high' agreement, Fig. 1a), despite an annual mean SST bias of 1.3 °C on the data-model comparison in this region (Table S1).
Because the SST amplitude and seasonal phasing are in agreement between the model and data (Table S2), the bias on the mean SST does not impact the seasonal pCO2 cycle (Pearson correlation coefficient > 0.5 and bias < 20 µatm on the seasonal pCO2 amplitude, Table 1). We also selected the Californian Current because it is a source of CO2 to the atmosphere, and 475 similarly to the East coast of the U.S, it ranks among one of the best monitored coastal regions in the world (e.g., Evans et al., 2011;Fennel et al., 2019;Hales et al., 2012;Turi et al., 2014). In this region, the model is classified as 'medium' agreement (Table 1 and Fig. 1a). Indeed, the simulated seasonal cycle of pCO2 is in relatively good agreement with coastal-SOM-FFN and Table 1), despite biases in the annual mean pCO2 compared to observations (Fig. 3f) and a phase shift in the seasonality of SSS and nutrients (Pearson correlation coefficient < 0.5). However, the Californian Current is also one of the 480 few coastal regions where an analysis of the processes controlling the pCO2 seasonality has already been performed using a regional biogeochemical model and sequential simulation removing processes one after the other (Turi et al., 2014), which can hence be compared to our analysis. Finally, the choice of the Norwegian Basin is motivated by the good performance ('high' agreement) of the model and the intense atmospheric CO2 sink that occurs in this contrasted region.

Seasonality along the East coast of North America 485
The seasonal evolution of pCO2 averaged over the East coast of the U.S (M10) is represented in Fig. 7a. Ocean pCO2 is minimum in winter (February/March ~ 331 µatm), it increases through spring and peaks in summer (August, ~ 400 µatm) before decreasing again in the fall. Figure 7b reveals the complex interplay of the four ocean internal processes (thermal, biology, ocean circulation, and fresh water flux) on the seasonal pCO2 variability which can either act in synergy or oppose each other. 490 The thermal effect (thermal, red line on Fig. 7b) increases pCO2 from early spring to summer by decreasing the solubility of CO2. In contrast, the solubility of CO2 increases in autumn and winter, inducing a decline in pCO2. The largest changes in pCO2 associated with the change in SST occur during spring (29 µatm month -1 in June) and fall (-26 µatm month -1 in November). This thermal effect was already identified by Signorini et al. (2013) in their observational study and further 495 confirmed by Cai et al. (2020). These authors highlighted that lowest pCO2 was generally reported in winter or at the beginning of spring and highest pCO2 in summer or autumn, despite significant temporal and spatial heterogeneity between the different sub-regions of the East coast of the U.S (Scotian shelf, the Gulf of Maine, the Georges Bank/Nantucket shoals, the Middle Atlantic Bight, and the South Atlantic Bight). The effect of biology above the mixed layer depth (bio, green line) reduces pCO2 throughout the year revealing that primary production exceeds organic matter degradation in the surface layer all year long. 500 The largest pCO2 decrease associated with biology is observed in early spring (values of -68 µatm month -1 in April) which is well documented (e.g., Shadwick et al., 2010Shadwick et al., , 2011Signorini et al., 2013). The transport of chemical species by ocean circulation (circ, blue line) increases pCO2 and tends to oppose biology year-round except at the end of fall/beginning of winter. This pCO2 increase induced by the circulation term is maximum in April (26 µatm month -1 ). Throughout the year, the contribution of fresh water fluxes (fw, pink line) remains minor compared to the other terms (maximum absolute value of 9 505 µatm month -1 in January). For each month/season, the air-sea CO2 exchange term (CO2 flux, black line) counteracts change in pCO2 associated with ocean internal processes taking place in surface seawater (sum of bio, circ, thermal and fw). The CO2 flux term increases pCO2 at the sea-surface (acting as an atmospheric CO2 sink) throughout the year except during summer (between July and September) where it decreases sea surface pCO2 and releases CO2 towards the atmosphere (acting as an atmospheric CO2 source). This simulated atmospheric CO2 uptake all year long except for the summer season is also in 510 agreement with previous literature (Fennel et al., 2019;Laruelle et al., 2015;Signorini et al., 2013). The study of Laruelle et al. (2015) has nevertheless shown that in spring, the southern part of the Eastern North American coast is quasi neutral and that in fall, some regions such as the Gulf of Maine or the Georges Bank acts as a CO2 source. The temporal change of pCO2 (pCO2 change, cyan line) is the result of the non-perfect balance between the internal processes and the air-sea CO2 flux.

515
We evaluate the rate of change tied to each process during the marked peak-to-peak pCO2 increase observed between winter and summer (from 331 µatm in February to 400 µatm in August, Fig. 7a). A positive rate of change (in µatm month -1 ) indicates that the process contributes to an increase in pCO2 between winter and summer (February-August). This process-based analysis reveals that the winter-to-summer pCO2 increase in the East coast of the U.S (M10) mainly results from thermal (rate of change = +5 µatm month -1 ) and ocean circulation (rate of change = +4 µatm month -1 ) influences combined with a large reduction of 520 the biological CO2 uptake (rate change of +7 µatm month -1 , Fig. 7b). The importance of the thermal and circulation effects as well as the presence of a strong biological drawdown are in line with results from past studies (e.g., Laruelle et al. (2015), Shadwick et al. (2010Shadwick et al. ( , 2011, Signorini et al. (2013) and Cai et al. (2020). Our results which identifies the reduction of biological carbon uptake as a key control of pCO2 seasonality agree with the studies of Shadwick et al. (2010Shadwick et al. ( , 2011, but slightly diverge compared to those of Signorini et al. (2013) or Laruelle et al. (2015), which found that the thermal effect was 525 the dominant driver. This difference is largely explained by the different levels of details in the decomposition method. While most model studies, including ours, use seasonal change in SST, SSS, DIC and ALK, observational approaches cannot isolate the compounding changes tied to biological activity from those of ocean transport.

Seasonality in the Norwegian basin and in the Californian Current
The pCO2 seasonal cycle in the Norwegian Basin (M16) and the Californian Current (M2) simulated by MOM6-COBALT are 530 represented in Fig. 7c and Fig. 7e, respectively. The Norwegian Basin shows a near-constant pCO2 value (~ 330 µatm) throughout the year except in spring when it drops by 30 µatm (minimum pCO2 value of 300 µatm in June). The phasing of the seasonal pCO2 cycle in the Californian Current is similar to that along the East coast of U.S, with a minimum pCO2 value of 366 µatm in March followed by an increase that reaches a maximum pCO2 value of 433 µatm in August and then decreases again at the beginning of the fall. 535 The decomposition of the seasonal cycle into different processes for both the Norwegian Basin and the Californian Current ( Fig. 7d and Fig. 7f) reveal patterns that are qualitatively similar to those already diagnosed for the East coast of the U.S (Fig.   7b). For both shelf regions, the biological and circulation effects respectively remain negative and positive throughout the year, while the thermal effect increases pCO2 in spring and summer but decreases pCO2 in fall and winter. The fresh water term is 540 also minor compared to the other terms. Quantitatively, however, the amplitude of the different terms points to different first order control in the pCO2 seasonality for each region. The amplitudes are calculated here using the marked peak-to-peak change in pCO2 which occurs between February and June in the Norwegian basin and between March and August in the Californian Current.
In the Norwegian basin, the strong winter to summer pCO2 decreases (43 µatm, Fig. 7c) is mainly associated with the large and rapid CO2 uptake associated with the spring phytoplankton bloom (biological rate of change = -45 µatm month -1 in average between February and June and with a maximum pCO2 uptake of -175 µatm month -1 in June, Fig. 7d). This biological drawdown is only partly compensated by the supply of high pCO2 water masses by the ocean circulation (rate of change = +24 µatm month -1 ). This dynamics is consistent with the fact that the Norwegian Basin is one of the most productive region of the 550 world characterized by a well-documented, intense spring bloom (e.g., Findlay et al., 2008). In addition, the effect of thermal changes only plays a comparatively minor role here (rate of change = +7 µatm month -1 ).
In contrast to the East coast of the U.S and the Norwegian Basin, the analysis performed in the Californian Current reveals that circulation is the main driver of the winter-to-summer pCO2 increases (68 µatm, Fig 7e). The upwelling of high-pCO2 555 waters increases surface pCO2 year-round. Its influence is however weaker in winter than in summer, thereby explaining the pCO2 increase observed between February and August (rate of change = +12 µatm month -1 , Fig. 7f). This large contribution from circulation is consistent with the simulations of Turi et al. (2014), which identified the ocean transport associated with upwelling in the Californian Current as the dominant process, and the higher intensity of the summer upwelling and its impact on pCO2 was also reported in prior work (e.g., Evans et al., 2015;Fiechter et al., 2014;Turi et al., 2014). In this region, biology 560 also opposes the effect of ocean circulation, with upwelled deep water bringing nutrients to the surface and stimulating phytoplankton productivity (e.g., Evans et al., 2015;Fiechter et al., 2014;Turi et al., 2014). However, it plays a minor role in the pCO2 increase (rate of change ~ 0 µatm month -1 ) as well as for the thermal effect (rate of change = +4 µatm month -1 ).

Conclusions
In this study, an OGCM (MOM6-COBALT) which is primarily designed for the open ocean was used to examined sea surface 565 pCO2 seasonality in the coastal domain. We first evaluated the ability of the model to reproduce the spatial and temporal dynamics of key environmental variables, such as SST, SSS and sea surface nutrients against in-situ observations. The spatiotemporal variability of coastal pCO2 was also evaluated using direct coastal pCO2 observations from the SOCAT database (SOCATv6, Bakker et al., 2016), and a global observational continuous monthly pCO2 climatology available at high spatial resolution (coastal-SOM-FFN, Laruelle et al., 2017). 570 Our model-data comparison showed a relatively good agreement on the environmental variables spatio-temporal distribution except for some coastal regions mainly located in marginal seas. Our results also revealed a relatively good agreement between pCO2 from MOM6-COBALT, coastal-SOM-FFN and SOCATv6, both in time and space, and most of the discrepancies between the three products are found in regions with poor data coverage, such as in the Bay of Bengal, The Sea of Okhotsk or 575 in the Hudson Bay (Fig. 1a). This study highlights the regions (Fig.1a, e.g., Indian ocean margins, Peruvian upwelling, marginal seas) where new observational data are most urgently needed, specifically data collected during different periods of the year that are currently missing to improve our understanding of the CO2 exchange between coastal regions and the atmosphere at the regional and global scales. From the model-data evaluation, we identified regions where the MOM6-COBALT model shows highest agreement in reproducing the spatial and seasonal pCO2 variability, and where the different processes governing 580 the pCO2 dynamics can be examined with reasonable confidence ('high' and 'medium' agreement regions in Table 1 and Fig.   1a).
We also adapted a novel method to quantify the contributions of the different physical and biological processes governing the sea surface pCO2 seasonality in the coastal domain. This method goes one step further than past coastal studies (e.g., Signorini 585 et al., 2013;Turi et al., 2014) where the processes attribution was only based on the seasonal changes in DIC, ALK, SST and SSS or/and combined with a series of sequential simulations isolating one term after the other. In particular, our simulations are non-sequential and allow accounting for the co-variations between the different variables impacted by each process and how their simultaneous evolution modulates in quantitative terms the pCO2 dynamics. Our approach, which is illustrated in three coastal regions (the East coast of the U.S, the California Current and the Norwegian Basin), allows to decipher the 590 complex interplay between ocean transport of chemical species (DIC, ALK and SSS), biological drawdown, fresh water fluxes (dilution/concentration effects) and thermal changes (air-sea fluxes and transport of temperature) on the pCO2 dynamics.
Depending on the season and region, these terms can reinforce or oppose each other and act to strengthen or dampen the amplitude of pCO2 seasonal variations that control the air-sea CO2 exchange. Along the East Coast of the U.S and in the Californian Current, pCO2 increases from winter-to-summer. In the former region, this increase is controlled by a subtle 595 balance between biological drawdown, thermal changes and ocean circulation, while in the Californian Current, the circulation due to the upwelling (supplying pCO2-rich waters to the surface) drives the increase in pCO2. In contrast, in the Norwegian Basin, biological drawdown dominates the marked spring pCO2 decrease observed in the region. These differences in the quantitative controls of pCO2 dynamics from one region to another support our proposed analysis at the broad scale of the 45 MARCATS regions that together compose the global coastal ocean. 600 A handful of observational-based studies analyzed the seasonal variability of pCO2 in the global coastal ocean (Cao et al., 2020;Chen and Hu, 2019;Laruelle et al., 2017). The mechanistic understanding of seasonal pCO2 variations was, and remains limited by the amount of available observations. The modeling approach tailored for the coastal ocean presented in this manuscript complements observational studies and help improve our quantitative understanding of the underlying physical 605 and biological drivers of the coastal pCO2 dynamics. The comparison of the model performance to a state-of-the-art coastal pCO2 database and continuous pCO2 data product also lends confidence in our model results for a large fraction of the global coastal domain. The coastal ocean is under tremendous anthropogenic pressure (e.g. climate, land-use change and agriculture, pollution, urbanization; e.g., Mackenzie et al., 2005;Regnier et al., 2013;Seitzinger et al., 2005). Understanding the interplay between physics, biology and thermal processes and how they control coastal pCO2 worldwide will be key to assess how their 610 future changes impact air-sea CO2 exchange in coastal environments.

Acknowledgements
L. Resplandy and E. Liao acknowledge the Cooperative Institute for Modeling the Earth System between NOAA GFDL and Princeton University, the Sloan Research foundation and the Princeton Catalysis Initiative. G. G. Laruelle is research associate of the F.R.S-FNRS at the Université Libre de Bruxelles. P. Regnier received financial support from BELSPO through the 615 project ReCAP, which is part of the Belgian research program FedTwin and from the European Union's Horizon 2020 research and innovation program VERIFY (grant agreement no. 776810) and ESM2025 projects

Data availability
The Surface Ocean CO₂ Atlas (SOCAT) is an international effort, endorsed by the International Ocean Carbon Coordination 620 Project (IOCCP), the Surface Ocean Lower Atmosphere Study (SOLAS) and the Integrated Marine Biosphere Research (IMBeR) program, to deliver a uniformly quality-controlled surface ocean CO₂ database. The many researchers and funding agencies responsible for the collection of data and quality control are thanked for their contributions to SOCAT. Every previous version of the SOCAT database can also be accessed from the following page: https://www.socat.info/index.php/previousversions/. The coastal-SOM-FFN pCO2 datasets description and dataset can be downloaded from Laruelle et al. (2017) and the 625 atmospheric CO2 concentration data (xCO2) derived from the Earth System Research Laboratory (Conway et al., 1994;GLOBALVIEW-CO2, 2004). The SST and SSS used for the evaluation the model were extracted from the NOAA OI SST V2 (Reynolds et al., 2007)        1: Model vs coastal-SOM-FFN agreement level. For each MARCAT, the agreement ('high', 'medium' and 'low') is attributed from the pCO2 spatio-temporal analysis. Regions where the model fulfils criteria on the annual mean and seasonality are labelled as 'high' agreement regions (i.e., annual mean mismatch < 20 µatm between MOM6-COBALT and coastal-SOM-FFN, Pearson 930 correlation coefficient > 0.5 and seasonal amplitude mismatch < 20 µatm, dots in Fig. 1a). High* agreement region can present a bias > 20 µatm on the comparison with SOCATv6 (see Table S1). 'Medium' agreement regions represent MARCATS where the model only fulfils seasonal criteria (seasonal amplitude and phase, dashed in Fig. 1a). Other's regions ('low' agreement) do not fulfil the two criteria associated to the seasonality (no symbol in Fig. 1a). Regions with 'high' agreement are considered as the most robust for an in-depth analysis of the processes driving the coastal pCO2 dynamics and are highlighted in bold on the