Downscaling CMIP6 Global Solutions to Regional Ocean Carbon Model: Connecting the Mississippi, Gulf of Mexico, and Global Ocean

Coupled physical-biogeochemical models can significantly reduce uncertainties in estimating the spatial and temporal patterns of the ocean carbon system. Challenges of applying a coupled physical-biogeochemical model in the regional ocean include the reasonable prescription of carbon model boundary conditions, lack of in situ observations, and the 10 oversimplification of certain biogeochemical processes. In this study, we applied a coupled physical-biogeochemical model (Regional Ocean Modelling System, ROMS) to the Gulf of Mexico (GoM) and achieved an unprecedented 20-year highresolution (5 km, 1/22°) hindcast covering the period of 2000-2019. The model’s biogeochemical cycle is driven by the Coupled Model Intercomparison Project 6-Community Earth System Model 2 products (CMIP6-CESM2) and incorporates the dynamics of dissolved organic carbon (DOC) pools as well as the formation and dissolution of carbonate minerals. Model 15 outputs include generally interested carbon system variables, such as pCO2, pH, aragonite saturation state (ΩArag), calcite saturation state (ΩCalc), CO2 air-sea flux, carbon burial rate, etc. The model’s robustness is evaluated via extensive model-data comparison against buoy, remote sensing-based Machine Learning (ML) predictions, and ship-based measurements. Model results reveal that the GoM water has been experiencing an ~ 0.0016 yr-1 decrease in surface pH over the past two decades, accompanied by a ~ 1.66 μatm yr-1 increase in sea surface pCO2. The air-sea CO2 exchange estimation confirms that the river20 dominated northern GoM is a substantial carbon sink. The open water of GoM, affected mainly by the thermal effect, is a carbon source during summer and a carbon sink for the rest of the year. Sensitivity experiments are conducted to evaluate the impacts from river inputs and the global ocean via model boundaries. Our results show that the coastal ocean carbon cycle is dominated by enormous carbon inputs from the Mississippi River and nutrient-stimulated biological activities, and the carbon system condition of the open ocean is primarily driven by inputs from the Caribbean Sea via Yucatan Channel. 25

DIC and alkalinity data, its reliability was questionable as temperature and salinity alone cannot fully describe the spatial and 65 temporal pattern of these carbon variables. Laurent et al. (2017) presented a regional model study of the eutrophication-driven acidification and simulated the recurring development of an extended acidified bottom waters in summer on the NGoM shelf.
They proved that the acidified waters were confined to a thin bottom boundary layer where the production of CO2 was dominated by benthic metabolic processes. Despite reduced Ω values being produced at bottom water due to acidification, these regions remain supersaturated with aragonite. Chen et al. (2019) presented a unified model to estimate surface pCO2 by 70 applying Machine Learning (ML) methods on remote sensing data and cruise pCO2 measurements. Their ML model confirmed that the GoM was a carbon sink. Recently Gomez et al. (2020) performed another GoM carbon model study covering the period of 1981-2014. Their model initial and boundary conditions were derived from a downscaled Coupled Model Intercomparison Project 5 (CMIP5) Modular Ocean Model (25km resolution, Liu et al., 2015). Their model results showed that GoM was a sink for atmospheric CO2 during winter-spring, and a source during summer-fall, producing a basin-wide 75 mean CO2 uptake of 0.35 mol m −2 yr −1 . Nevertheless, their model does not include DOC pool or calcification process, which are imperative to describe the dynamic of DIC and alkalinity in the ocean.
Despite the above carbon system regional modeling efforts, we notice that several processes that could contribute significantly to the carbon cycle in the GoM have not been investigated yet. The carbon cycle in the ocean is linked with the nutrient cycle 80 through photosynthetic activities, calcification, and OM remineralization (Anav et al., 2013;Farmer et al., 2021;King et al., 2015). OM remineralization could be the most critical mechanism regulating the ocean carbon system, followed by the CaCO3 cycle (Lauvset et al., 2020), with the remineralization of small detritus accounting for over 40% of the DIC production on the shelf (Laurent et al., 2017). Autochthonous nutrients from direct remineralization of OM determine the gradient of DIC in the euphotic layer (Boscolo-Galazzo et al., 2021;Boyd et al., 2019). During this process, the fast-sinking of OM, higher OM 85 particulate to dissolved ratio foster larger sedimentation rate and greater DIC removal of the euphotic layers; on the contrary, slower sinking and faster decomposition rate of OM favors nutrient and DIC retention in the euphotic layers (Davis et al., 2019;Mari et al., 2017;Turner, 2015). The remineralization of land-derived OM and CaCO3 precipitation are significant factors controlling air-sea CO2 flux (Mackenzie et al. 2004). Studies have revealed the formation of marine CaCO3 (Burton and Walter,1987;Inskeep and Bloom, 1985;Zhong and Mucci, 1989;Zuddas and Mucci, 1998)

and the dissolution of marine 90
CaCO3 mineral is Ω-dependent as well (Adkins et al. 2021). The Ω will be depressed with more CO2 dissolves in seawater and can be used as an indicator for the buffering capacity of the ocean carbonate system. Given that Ω influences the calcification rate of marine organisms and regulates the acidity of bottom waters, it should be considered in the CaCO3 cycle for a comprehensive carbon cycle assessment.

95
In this study, we dynamically downscaled the solution of the latest Coupled Model Intercomparison Project 6 (CMIP6) product to a regional model to inherit the climate perturbation signals (Liao et al., 2020) and the accumulative effect of carbon variables from the global solutions. Our regional model includes critical carbon cycle processes lacking in previous efforts, including https://doi.org/10.5194/bg-2021-339 Preprint.  the most up-to-date carbonate chemistry thermodynamic parameterization, phosphate cycling, formation & dissolution of CaCO3, and the inclusion of the DOC as a semi-labile carbon pool. The objective of this study is 1) to assess the feasibility 100 and robustness of utilizing downscaled global model products to drive a regional coupled physical-biogeochemical model, and 2) to examine the temporal trend of key variables of the carbon system (pCO2, pH, air-sea CO2 exchange, and Ω) of the surface ocean in the GoM. In addition, to evaluate the impact of MARS and the global ocean on GoM's carbon cycling, two perturbed experiments are designed. The following sections are organized: model setup is given in Sect. 2; in Sect. 3, we validate the model's performance against buoy, remote sensing-based ML solution, and ship-based measurements; the trend of key carbon 105 systems over the past two decades are given in Sect. 4; an assessment of the contribution of riverine inputs and global ocean are discussed in Sect. 5, together with a model uncertainty analysis.

Model setup
Our model is built on the platform of Coupled Ocean-Atmosphere-Wave-Sediment Transport modeling system (COAWST; 110 Warner et al., 2010). COAWST is an open-source community model which includes the Model Coupling Toolkit to allow data exchange among three state-of-the-art numerical models: Regional Ocean Modelling System [ROMS, svn 820, Haidvogel et al., 2008;Shchepetkin and McWilliams, 2005], the Weather Research and Forecasting model [WRF, Skamarock, et al., 2005], and the Simulating Waves Nearshore model [SWAN, Booij et al., 1999]. The carbon model presented in this study is based on a well-validated coupled physical-biogeochemical model by Zang et al. (2019 and, which covers the entire GoM waters 115 (Gulf-COAWST, Fig. 1). The Gulf-COAWST has a horizontal grid resolution of ∼5 km and 36 sigma-coordinate (terrainfollowing) vertical levels. A third-order upstream horizontal advection and fourth-order centred vertical advection are used for momentum and tracer advection. The biogeochemical model is developed based mainly on the pelagic N-based biogeochemical model Pacific Ecosystem Model for Understanding Regional Oceanography (NEMURO, Kishi, et al., 2007;Kishi et al., 2011). In this study, we extend the number of the NEMURO's state variables from seven to seventeen, including 120 nitrate (NO3), ammonium (NH4), two phytoplankton groups (small phytoplankton, diatom), one explicit phosphorus pool PO4, two pools of silicon (particulate silica, and silicic acid), three zooplankton groups (micro-zooplankton, meso-zooplankton, and predator-zooplankton) and two pools of detritus containing carbon (particulate organic nitrogen [PON], particulate inorganic carbon [CalC], a dissolved organic carbon (DOC) pool, a dissolved organic nitrogen pool (DON), oxygen (O2) to regulate the availability of aerobic OM remineralization, dissolved inorganic carbon (DIC) for bookkeeping the sum of dissolved inorganic 125 carbon species in water, and total alkalinity (TA) for bookkeeping the acid neutralization capacity of water. The stoichiometry between carbon and nitrogen in the OM production and remineralization is set to 6.625 following Fennel (2008).  The revised biogeochemical model incorporates key processes regulating the carbon model, including primary production, river DIC, PON and DOC delivery, sediment carbon burial, CO2 air-sea flux, CaCO3 cycling, and OM remineralization ( Fig.  135 2). Widely used carbon system variables, such as pCO2, pH, Ω, etc., are used as carbon system state indicators. The carbon module that takes in DIC, TA, PO4, Si (dissolved inorganic silicon), salt, temp, for calculating pCO2, pH, ΩArag, and ΩCalc, largely followed the recommended best practices (Dickson, Sabine, & Christian, 2007;Eyring et al., 2016;Orr et al., 2017;Zeebe and wolf-Gladrow, 2001), with an updated parameter prescription for dissociation constants for carbonic acid (K1) and bicarbonate ion (K2) (Millero, 2010), and solubility products for aragonite KA and calcite KC (Mucci, 1983). 140 Inorganic carbonate mineral (mainly CaCO3) forms during the photosynthetic activities of some phytoplankton species and foster aggregation of detritus and their sinking. The rate of CaCO3 production followed a dynamic ratio regarding primary production of small phytoplankton with low-temperature inhibition and enhancement during bloom conditions (Moore et al., 2004). The production and dissolution of CaCO3 are important processes for ocean acidity regulation, as its production (by 1 145 unit) nominally takes away a unit of [CO3 2-] from water, which reduces the alkalinity and DIC by 2 and 1 unit, respectively. This process routinely happens during photosynthetic activities of some phytoplankton species (such as coccolithophores, parameterized implicitly as a portion of small phytoplankton in this model) and other marine calcifiers. Carbonate minerals produced in the euphotic zone could be treated as equivalent storage of alkalinity and are usually transported towards the ocean sediment through sinking. Aragonite and calcite are two common mineral phases of CaCO3 secreted by marine organisms and 150 are included in the model. ΩCalc and ΩArag are calculated as the equilibrium product of Ca 2+ and CO3 2-. When Ω > 1, calcification is thermodynamically favored, and when Ω < 1, dissolution is thermodynamically favored. As defined in Eq. (1) (Millero, 1982;Millero, 1995), and [CO3 2-] is calculated through the model carbon module. Ksp is the stoichiometric solubility product mainly dependent on pressure, temperature, and salinity. Ksp is defined for aragonite and calcite as KA and calcite KC, 155 respectively.
In our model, the sediment pool of sinking particles is a simplified representation of burial and benthic remineralization processes, where the flux of sinking materials out of the bottommost grid point is added to the sediment pool and entered the 160 burial pool (remain inactive) with a dynamic ratio, the active sediment pools undergo enzyme-aided decomposition at rates regulated by temperature and oxygen, and release corresponding influx of ammonium, DIC, and alkalinity at the sediment/water interface. Our model uses a CO2 production ratio of 0.138 between sediment aerobic respiration and denitrification (Fennel et al., 2006) and an alkalinity production ratio of 1.93 between pyrite burial and denitrification . Upon being sunk to acidified regions, the dissolution of CaCO3, regulated by Ω, can consume dissolved CO2 and 165 neutralize the acid.
The bulk formula for air-sea gas exchange for CO2 is used following Wanninkhof (1992). Air-sea CO2 flux is calculated as Eq. (2).
Where FCO2 is the air-sea CO2 flux in the unit of mmol CO2 m -2 d -1 . Sc is the Schmidt number (nondimensional) (calculated following Wanninkhof, 2014), s is the solubility of CO2 in seawater in mol CO2 m -3 µatm -1 (calculated following Weiss, 1974), and ∆ pCO2 is the air and sea pCO2 difference in µatm. The term k660 is the quadratic gas transfer coefficient in cm h -1 (converted to m d -1 ). We calculated the air-sea CO2 flux using the relationships of gas exchange with wind speed at 10 m over 175 the sea surface (U10) following Wanninkhof (1992). We used the ocean convention for the CO2 flux, i.e., a positive flux is defined as the ocean being a sink of atmospheric CO2. Air pCO2 level is prescribed using fitted curve from column-averaged dry-air mole fraction of atmospheric carbon dioxide from 2002 to present derived from satellite product (merged dataset from SCIAMACHY/ENVISAT, TANSO-FTS/GOSAT, and OCO-2 [https://cds.climate.copernicus.eu/; Dils et al., 2014]). We performed a 20-yr model hindcast covering the period of 1 January 2000 -31 December 2019. The model setup was similar to that of Zang et al. (2020), with ocean physical initial and boundary conditions downscaled from the 1/12° data assimilated Hybrid Coordinate Ocean Model (HYCOM/NCODA, GLBu0.08/expt_19.1, expt_90.9, expt_91.0, expt_91.1, expt_91.2, and expt_93.0 [https://www.hycom.org; Chassignet et al. 2003]). Atmospheric forcings include ground level or sea 185 surface downwelling shortwave/longwave radiation, ground-level or sea surface upwelling shortwave/longwave radiation, surface air pressure, surface air temperature, relative humidity, precipitation, wind at 10 m were extracted from the NCEP Climate Forecast System Reanalysis (CFSR) (Saha et al. 2010) and Climate Forecast System Version 2 (CFSv2) (Saha et al. 2011 historical,r1i1p1f1,190 nominal resolution 100 km, [Danabasoglu, 2019]). Due to the global model ending in December 2014, the biogeochemical boundary condition of 2014 was used repeatedly for the period of 2015-2019. O2 is absent from the CESM2-WACCM-FV2 source, and is interpolated from the World Ocean Atlas 2018 (WOA18) product García et al., 2019).
Freshwater and terrestrial nutrient inputs from 47 major rivers discharged to the GoM are applied as point sources. River discharge and water quality data of rivers in the U.S. are collected from US Geological Survey (USGS) stations 195 (https://maps.waterdata.usgs.gov), river DOC is prescribed following the values reported by Shen et al. (2012), with additional reference from several other works (Reiman & Xu, 2019;Stackpoole et al., 2017;Xu & DelDuco, 2017).
Mexican river discharge data are collected from BANDAS (https://www.gob.mx/conagua), water quality data of Mexican rivers is prescribed as the average of that of the Mississippi and Atchafalaya rivers. River nutrient and carbon load are reconstructed from available USGS observations. Missing river alkalinity values are interpolated from climatological values, 200 and missing river DIC values are calculated from pH and alkalinity using the MATLAB program CO2SYS (Lewis & Wallace, 1998). Validations of the model's performance in physics, nutrient cycle, and primary production can be found in Zang et al., 2019 and 2020. In this study, we focus on the model's performance in the carbon cycle, which is presented in the next section.
The simulation is set to start on 2000-01-01 and end in 2019-12-31. Since model simulated DIC concentration in the water 205 column is sensitive to initial conditions (Hofmann et al., 2011;Xue et al., 2016), using initial condition from a snapshot

Validation
This section focus on the validation of the model results via comparison against autonomous mooring systems with surface pCO2 measurements, ship-based measurements from the Gulf of Mexico Ecosystems and Carbon Cruise transects (GOMECC, Barbero et al., 2019;Wanninkhof et al., 2013;Wanninkhof et al., 2016), and pCO2 underway measurements (data downloaded 220 from https://www.ncei.noaa.gov/access/oads/). Direct observations of the GoM carbon system have been recognized as unbalanced among seasons due to fewer data points available in winter compared to other seasons. To overcome the sporadic direct measurement dataset, we also performed a model-data comparison against the remote sensing-based ML product of sea surface pCO2 by Chen et al. (2019).

Model-buoy comparisons 225
Temporal variability of sea surface pCO2 was recorded by the autonomous mooring system at two sites (CoastalMS and coastal to 2020-08-12; 2020-08-12 to 2021-08-25; 2021-08-25 to 2021-11-29) is mutually influenced by the Mississippi River and the coastal ocean. The high-frequency measurements provide a time-resolved picture of year-round changing pCO2 values.
Temperature and salinity can influence the chemical equilibrium in the carbonate system, therefore shifting the pCO2 values.
Validating the temperature and salinity at these two mooring sites is a prerequisite before looking into the surface pCO2 levels. 235 In Fig. 3, the top four panels compare the sea surface temperature (SST) and salinity (SSS) between model and buoy measurements and show matching seasonal variability. At CoastalMS, the range for sea surface pCO2 is 150~600 µatm. Sea surface pCO2 records are more volatile at CoastalLA with a maximum value > 800 µatm and a minimum value around 150 µatm. Following the salinity drop, pCO2 at the CoastalMS site is simultaneously reduced, demonstrating the river's influence on both salinity and pCO2. At CoastalLA, however, the pCO2 level does not necessarily follow the trend of salinity, implying 240 complex controlling factors in addition to the river inputs. The bottom two panels of Fig. 3 show a good agreement between measured and simulated sea surface pCO2. We notice model-data discrepancies in April 2018-04 at CoastalLA and July 2011 at CoastalMS and ascribe such bias to the uncertainty in the riverine DIC inputs prescription and the limited model horizontal resolution (~5 km).

Model-cruise comparisons
Cruise carbon measurements include underway water pCO2 data and conductivity-temperature-depth (CTD) bottle results.   Figure 4 shows the vertical profiles of observed DIC, TA, and their ratio collected at the LA transects (-90°W, 27.5°N -29.1°N, shown in Fig.1) overlaid with the model solution. All three transects were taken during August when nutrient supply from the MARS was high. Model simulated profiles at the transects match well with the in situ CTD data. Relative low surface DIC concentration (< 2150 mmol m -3 ) above 200 m isobath demonstrates 260 the river's influence in NGoM. The general increasing trend of DIC with depth confirms the presence of a biological pump, where inorganic carbon is utilized during photosynthesis in the euphotic layer. Subsequently, the generated OM sinks into deeper waters while being remineralized along the way. The TA profiles show more variation compared with DIC, where generally a lower TA concentration (< 2380 millequivalent m -3 ) could be found at the surface as the direct dilution from river discharge, follow by a quick increase to ~2380 millequivalent m -3 in the euphotic layer due to the active photosynthetic 265 activities, which generate alkalinity. Further deep, the TA profiles show a decreasing trend between 200 and 700 m, which could be explained by the water column respiration and nitrification. The TA profiles show a slow increase from 800 m and deeper, which coincides with the alkalinity generation processes in sediment and possibly dissolution of carbonate minerals, both adding to the bottom water TA. The TA/DIC ratio has a maximum at the surface due to low DIC concentration and decreases with depth as DIC concentration increased. 270

Model-ML pCO2 product comparisons 275
Direct comparison between cruise measurement of ocean surface pCO2 and daily averaged model result might suffer from systematic bias due to the sparsity of curies data, both temporally and spatially. The ML model generates surface pCO2 from demonstrating substantial influence from the Loop Current. Compared with the ML model, our model produces lower pCO2 estimates over NGoM during winter and fall, higher pCO2 estimates over WF during summer, and stronger influence from the Caribbean Sea. Chen et al. (2019) reported that no satellite data was found for pCO2 <145 μatm or >550 μatm during their model development. This can also be a factor when considering the differences between the two products. 285 Besides buoy records, transects, and the ML products, we also perform an extensive model-data comparison using available ship-based underway pCO2 measurements from the Ocean Carbon and Acidification Dataset (https://www.ncei.noaa.gov/access/oads/). Our model can capture both spatial and temporal variability in the underway pCO2 290 dataset as well ( Fig. S1 and Fig. S2). These extensive model-data comparisons give us the confidence that our model, driven by carbon boundary conditions from the global model, can reproduce temporal, spatial, and vertical variability of the CO2 dynamics in the GoM.

Result
In this section, we present the spatial and temporal pattern of key carbon system variables, namely pCO2, pH, Ω, and air-sea 295 CO2 flux simulated over the past 20 years in the GoM. In this study, we emphasize the surface carbon condition in two regions: NGoM and Open GoM, where most existing in situ data are distributed. We perform a linear fit of the time series of the key carbon system variables in each region and show the fitted relationships in Fig. 6. The slopes of the fitted linear plot give estimations on the change rate of each carbon variable over the past two decades.

pCO2
We simulate a generally increasing trend in surface pCO2 level for both NGoM and Open GoM, with an increasing rate of 1.61 µatm yr -1 and 1.66 µatm yr -1 , respectively (Fig. 6). Seasonal ocean surface pCO2 variation is primarily affected by temperature 305 variations. To evaluate the pCO2 trend without temperature effects, we decouple the thermal and nonthermal components of pCO2 at the ocean surface using Eq. (4) and (5) and further extract the pCO2 variation due to gross primary production and airsea CO2 flux. The temperature sensitivity of CO2, γT = 4.23% per degree Celsius, proposed by Takahashi et al. (1993), is used to perform the thermal decoupling by multiplying the mean pCO2 value with SST induced difference described in Eq. (4). The nonthermal counterpart is obtained by removing the thermal effect from the pCO2 time series using Eq. (5). pCO2 variations 310 due to gross primary production are estimated from the carbon module based on the DIC consumed by gross primary production and denoted as pCO2 GPP . pCO2 variations due to air-sea CO2 flux are calculated from the carbon module based on the DIC change from the air-sea exchange and denoted as pCO2 flux . Figure 7 shows the seasonal and spatial patterns of four decoupled pCO2 components, namely the pCO2 th , pCO2 nt , pCO2 GPP , and pCO2 flux . The pCO2 th patterns in the top row (a)(b)(c)(d) of Fig. 7 reflect the fluctuation of pCO2 due to thermal effects.
Over the four seasons, a general pattern of rising pCO2 th from spring to summer and a gradual reduction from summer onwards can be observed. The NGoM shelf exhibits the lowest pCO2 th values during winter, while WF shows elevated pCO2 th values during summer. The higher pCO2 th values dwelling in the southern part of the Yucatan shelf reveal the warm water flowing 320 into the GoM from the Caribbean Sea. The second row (e)(f)(g)(h) of Fig. 7 shows the nonthermal component of pCO2. The relatively high pCO2 nt during winter on the NGoM shelf, compared to that of the Open GoM, shows the strong solvation effect of CO2 with low SST, contributing to a high DIC/TA ratio and strong carbon uptake. The lower two rows of Fig. 7 show pCO2 changes due to the gross primary production and CO2 air-sea exchange, respectively. The pCO2 GPP reflects the intensity of primary production in terms of pCO2 reduction. pCO2 GPP has larger magnitudes in NGoM during spring and summer and is 325 gradually attenuated during and after fall. The large magnitudes of pCO2 GPP during summer in NGoM waters and Open GoM region following the extension of MARS plume, suggesting strong biological CO2 removal in those regions. These results show that gross primary production has a stronger regulation on surface pCO2 during spring and summer in river-dominated waters and upwelling regions. While a minor contribution from gross primary production can be seen during winter, on the flat and shallow WF, and in Open GoM regions southern of Loop Current. The pCO2 flux reflects the intensity of air-sea CO2 330 exchange attempting to mitigate the air-sea disequilibrium caused by local physical and biological processes. The relatively high value of pCO2 nt and low magnitude of pCO2 GPP , as well as low river discharge (minimal river water mixing) in the WF during winter, indicate a strong CO2 uptake from the atmosphere due to low SST. This analysis agrees with the low pCO2 th and the high pCO2 flux values in the WF during winter shown in Fig. 7(d)(p). Situations during seasons other than winter are more complicated due to active biological activities and mixing. The Mississippi delta region has a high pCO2 th value during 335 summer, however, combined with effects of mixing and strong primary production (large magnitude of pCO2 GPP ), this region acts as a strong carbon sink that exhibits a high value of pCO2 flux compensated from the atmosphere. Figure 7 demonstrates that most of the time around the year, the surface pCO2 pattern is not dominated by a single factor but a synthesis of multiple controlling processes. The result of pCO2 decomposition agrees with the current view of the pCO2 dynamic and carbon uptake patterns in the GoM, which is strong carbon uptake during winter due to thermal effect and high biological CO2 drawdown 340 during spring and summer under the riverine influence.

pH 345
Ocean surface pH in the GoM shows a clear decreasing trend, with a 0.0020 yr -1 decrease over the NGoM region and a 0.0015 yr -1 decrease over the Open GoM region. Figure 8 shows the seasonal pattern of ocean surface pH over the GoM. Spatial and seasonal pH patterns show larger variation over the NGoM, especially on the inner shelf (depth<50m). The pH level in the surface water is closely associated with temperature, photosynthetic activities, and water mixing. The high pH value on the NGoM shelf reveals the strong influence of riverine alkalinity export and nutrient-stimulated primary production.

Aragonite and Calcite Saturation State
Aragonite undersaturation occurs ([CO3 2-] < 66 μmol kg⁻¹) before calcite undersaturation ([CO3 2-] < 42 μmol kg⁻¹) (Feely, Doney, & Cooley, 2009;Feely et al, 2002), as a result, ΩCalc is approximately 50% higher than ΩArag and their spatial and 360 seasonal variations are very similar as shown in Fig. 9. Variations in temperature, alkalinity, and pCO2 impose important controls on ΩArag. The multiyear variability of ΩArag at the ocean surface is shown in Fig. 6(c). The NGoM region shows a smaller decreasing trend in ΩArag (0.0045 yr -1 ) compared to that of Open GoM (0.0068 yr -1 ). Noted the data contained in Fig.   6 does not include water from the shallow shelf waters (water depth < 10 m), therefore the trend in NGoM does not incorporate the condition in coastal estuaries. The spatial distribution of ΩArag across the GoM depicts a healthy level of Ω and a low risk 365 of ocean acidification (Fig. 9). While the coastal ocean generally has a relatively high Ω level, some coastal locations warrant special attention when evaluating their tendency towards calcium mineral dissolution. These locations include coastal regions that experience a large load of riverine OM inputs (e.g., the Mississippi River delta in summer) and the upwelling regions that receive relatively higher acidity water from the bottom ocean (e.g., west of Yucatan). These regions show significant Ω reductions compared to the surrounding waters and are the potential primary victim of ocean acidification. The influence from 370 the river on Ω is complex, on one hand, a high nutrient level of river discharge could stimulate a high photosynthetic rate which consumes DIC and increases Ω, on the other hand, photosynthesis favors calcification which consumes carbonate ion and reduces Ω. Therefore, Ω is subject to increase with stronger photosynthesis and decrease with stronger calcification. Therefore the magnitude of the overall effect will depend upon photosynthetic rates and the calcification rate. In this work, two phytoplankton groups are modeled, diatom (silicifying), and small phytoplankton (implicitly include the calcifying 375 coccolithophores, foraminifera, and dinoflagellates), of which only the small phytoplankton group has the potential to conduct calcification (Raven & Giordano, 2009). Besides being regulated by temperature and small phytoplankton concentration, the calcification rate also depends on the composition of the phytoplankton population. The small phytoplankton group has a survival advantage at relatively low nitrogen concentrations and could be grazed by two zooplankton groups (mesozooplankton and micro-zooplankton), whereas diatom is more nutrient-demanding and can be grazed by three zooplankton 380 groups (predator zooplankton, meso-zooplankton, and micro-zooplankton). The competitive phytoplankton evolution shaped the relative rates between photosynthesis and calcification on the NGoM shelf during summer. The reduced Ω to the east of the Mississippi River delta is a combined result from high DIC water entrained by Loop Current eddies west of the delta, and an increased ratio of small phytoplankton in offshore waters.

Air-sea CO2 Flux
Air-sea CO2 flux is calculated from daily averaged model data over 2001 to 2019 ( Table 2). The GoM overall is a CO2 sink with a mean flux rate of 0.62 mol C m -2 yr -1 (11.77 Tg C yr -1 ), which is commensurate with the reported value of 11.8 Tg C 390 yr -1 by Coble et al. (2010). The greatest carbon uptake rate occurs in winter (1.97 mol C m -2 yr -1 ), while the weakest carbon uptake is present during fall (0.16 mol C m -2 yr -1 ). The strongest carbon efflux is simulated in summer (-0.57 C m -2 yr -1 ). On average, water in the NGoM acts as a sink throughout the years, and the water in the Open GoM acts as a weak source during summer (and fall for 2002, 2004, 2006, and 2009) and a sink during the rest of the year. The direction and magnitude of the air-sea exchange can be seen in Fig. 10, where a positive number indicates the ocean is a carbon sink. The NGoM is a very 395 strong CO2 sink year-round and Open GoM is a source of CO2 during summer but a sink in the rest of the year (except during fall in a few years), as shown in Fig. 10(b). There are clear trends and patterns on multiyear CO2 air-sea flux as shown in Fig.   10 In the last three rows of Table 2, the annual air-sea CO2 flux is listed together with that reported by Xue et al. (2016)  with the flux estimates of Gomez, et al. (2020), we rescale our estimates to the gas transfer velocity parameterization used in their work (Wanninkhof, 2014) and produced mean estimations of 1.59, 0.52, 0.50 mol m -2 yr -1 for NGoM, Open GoM and Gulf-wide, respectively. Gomez et al. (2020) simulate smaller carbon sink values possibly due to the overestimation in NGoM shelf surface pCO2 (especially near the MARS plume region and WF, the model-data discrepancy can be viewed in the supporting information Fig. S6 in their paper). The observation-based studies by Huang et al. (2015) yielded an annual sink 415 of 0.96 ± 3.7 mol C m -2 yr -1 , based on the monthly satellite products QuikSCAT wind data (12.5 km resolution). Lohrenz et al.
(2018) estimated an annual sink of 1.1 ± 0.3 mol C m -2 yr -1 , using gas transfer velocities estimated for each 8-day period. These two observation-based studies suggest a larger carbon sink than previous model work and are more consistent with the estimation of this work. Another reason that this study estimated a stronger air-sea flux is because previously reported values are based on an earlier timespan when the ocean was a weaker carbon sink due to a relatively low atmospheric pCO2 compared 420 with more current conditions. After comparing with multiple sources, we believe the air-sea CO2 flux estimated in this work is more reliable than previous GoM model studies.

Discussion
In this study, we demonstrate that the regional high-resolution carbon model can reproduce the spatial and seasonal patterns of ocean surface pCO2 and generate reliable TA/DIC profiles on the NGoM shelf. We have detected a consistent acidification trend in the GoM over the past two decades. In this session, we further diagnose river discharge and global ocean's impacts 430 on the GoM carbon system via a comparison between the control experiment (His) and the two perturbed experiments (Bry and NoR), followed by a discussion of model uncertainty.

Contribution from River and Global Ocean
Two sets of perturbed experiments (Bry, NoR) are performed to complement the control experiment (His). One with the 435 clamped boundary conditions that repeat the DIC and TA level as that of the year 2000, the other with the river forcing removed to examine the impact of fluvial input on the coastal carbon system. Figure 11 shows the multiyear averaged rate of change and mean levels of the four carbon variables (pCO2, pH, ΩArag, and CO2 flux) simulated by the three experiments. Table 3 summarizes the mean levels of pCO2 over the NGoM and Open GoM. The definition of the pCO2 th , the pCO2 nt , the pCO2 GPP , and the pCO2 flux can be found in Sect. 4.1. The pCO2 mixing is defined in Eq. (7), which reflects the pCO2 level due to the 440 heterogeneous water mixing. It can also be considered as the pCO2 level determined by the water with a multiyear mean temperature and without the influence from gross primary production or air-sea CO2 flux.
The most salient difference among the three experiments is the significant elevation of annual mean pCO2 level (in Fig. 11) in 445 the NGoM during the NoR experiment, combined with a significantly reduced carbon sink during summer (in Fig. 12, from 0.287 to -0.093 Tg season -1 using His as a benchmark). The difference can be better resolved by the pCO2 decomposition results shown in Table 3, where a drastic change in water carbon system emerges in NGoM during the NoR experiment (as compared to the other experiments with river input), evidenced by the large pCO2 mixing value deviated from that of His and Bry experiment in the NGoM region during spring, summer, and winter. The low values of pCO2 nt in NGoM during summer can 450 be explained by a strong biological drawdown of CO2 associated with the high productivity fuelled by riverine nutrient supply. dominating impact of river input to the NGoM carbon system in terms of gross primary production, surface pCO2 level, and air-sea CO2 exchange.

460
The Bry experiment is intended to resemble a GoM carbon system evolving with local carbon forcing without the influence of secondary carbon accumulation transported from the global ocean. From Fig. 11, we can see the directions (+/-signs) of the trend of all four carbon variables are preserved for all three experiments. This means under both perturbations, the regional ocean continued to have an acidifying tendency with increasing pCO2, decreasing pH, decreasing ΩArag, and increasing CO2 air-sea flux. This result is intuitively expected as the local atmospheric pCO2 keeps increasing at a faster rate than that of the 465 ocean, which should determine the direction of ocean carbon system change. Results from Table 3 show a close resemblance in the magnitude and seasonal pattern between Bry and His experiment in the Open GoM region, with a small yet steady reduction in pCO2 mixing for the Bry experiment among all seasons. The small reductions in pCO2 mixing of the Bry experiment compared to that of His reflect the contribution from extraneous carbon accumulation from the global ocean that is included in the His experiment. As expected, slightly greater CO2 sink values are reported in Fig. 12 for Bry than His. Since the Bry 470 experiment has a smaller carbon accumulation in the Open GoM region compared to that of the His experiment, the ocean surface requires a slightly greater carbon uptake to reach equilibrium with the atmosphere. The acidification trend results can be seen from Fig. 11(b), where intensified acidification trends are simulated for His experiment (-0.0015789 yr -1 ) compared with that of Bry (-0.0015396 yr -1 ) and NoR (-0.0015243 yr -1 ). Combining the results from the three experiments, we conclude that, in addition to elevated atmospheric CO2 level, both inputs from MARS and global oceans contribute to the overall 475 acidification trend in GoM, with the impacts from MARS mainly limited to the NGoM shelf region and global ocean's impacts spanning in the Open GoM.  A likely warmer climate combined with heavier precipitation and greater river discharge is predicted in the following years for MARS (Dai et al., 2020;Fischer & Knutti, 2015;Frei et al., 1998;Tao et al., 2014), although climate change might reduce precipitation for some low and middle latitudes regions (Arora & Boer, 2001;Na, Fu, & Kodama, 2020). A warmer climate will reduce the momentum of the Loop Current, and less tropic water (reduce by about 20-25%) will be introduced into the 490 GoM from the Caribbean (Liu et al., 2012). As a consequence, Loop Current might penetrate less into the NGoM and reduce the upwelling along the NGoM and WF slope. Stronger river discharge with nutrient loads will exacerbate the NGoM acidification in bottom water (Laurent et al., 2018) while increasing the surface water biological CO2 utilization and removal, creating larger river plume regions that exhibit distinct carbon footprint compared to its surrounding waters. Such predictions resemble the perturbation prescribed in the Bry and NoR experiments, where reduced global ocean impact can be assessed by 495 the difference between Bry and His experiments, and impacts from increased river discharge can be assessed by the difference between His and NoR experiments. We anticipate the Open GoM to be a stronger carbon sink in the future under the projection of Loop Current weakening. And the NGoM will continue to be a strong carbon sink with the sink region expanded in response to predicted greater river discharge and smaller momentum in Loop Current.

Model-Data Discrepancy and Model Uncertainty
Field samples of the carbon system give us synoptic knowledge of the carbon cycle in the ocean. However, carbon system attributes are subject to large fluctuation due to temperature, salinity, mixing, and biological activities, current observations of 505 the carbon system at the sea surface or vertically along transects are far from enough to reveal the carbon system evolution in the GoM. As ship-based observations are limited by the spatial coverage and temporal coverage, mooring observations have a high-frequency (~3 h) in time but only cover limited geological locations. ML model derived from remote sensing and cruise data inherit the bias from satellite ocean color products and ship-based measurements, and more importantly, it assumes the training data contains all information that defines the system it is trying to predict, which is not necessarily true. One benefit 510 of the numerical models is to offer information to bridge fragmentary knowledge, fill in the gaps between observation and reality. However, the marine carbon cycle is admittedly a complex process. Several simplifications and parameterizations are needed to perform a numerical simulation. Nevertheless, specifications for some key processes may warrant further investigation and better parameterization: 1) The multiple alkalinity generation processes in the sediment pool in current experiments are linked linearly with the aerobic decomposition of PON with a fixed ratio, which can potentially induce large 515 bias during high PON concentration. The anoxic zone chemistry component can be added to properly simulate the carbonate system in oxygen-deficient conditions (Raven, Keil, & Webb, 2021), which can prevail in bottom boundary layer waters in coastal regions in NGoM, especially during summer. Adding in anoxic zone chemistry will also allow a more diversified prescription for TA generation, which plays a key role in the understanding of sediment pH dynamics (Gustafsson et al. 2019;Middelburg, Soetaert & Hagens, 2020). 2) In our model, the density-related fragmentation/flocculation of detritus OM is 520 simplified with one particulate and one dissolved pool, each with a fixed sinking rate. Coagulation and flocculation can transform DOC into particulate OM, or subsequently form large aggregates, whose remineralization rate can be much faster (Ploug et al., 1999). The remineralization-sinking dynamic determined the fate of OM decomposition (and water column DIC profile) and should be allowed to have more degree of freedom in future model development. 3) Calcification in this work can reflect the primary factors regulating marine calcification. However, important feedback from water acidity on the calcification 525 is omitted due to the overall supersaturation with aragonite in GoM shelf waters. Therefore, the modeled CaCO3/PON ratio could not reflect the decreasing trends of the CaCO3/PON ratio under acidification (Zondervan et al., 2001). 4) Phytoplankton groups can play different roles in carbon cycling given their different sizes, sinking rates, calcification rates, etc., and their relative ratio would be critical to the carbon dynamic (Le Moigne et al. 2015;Poulton et al. 2007). The interplay between zooplankton grazing and phytoplankton bloom in this work captured the seasonal dynamic but only had fixed modes toward 530 nutrient levels. High nutrient concentration favors the success of diatom, and lower nutrient level gives small phytoplankton a competitive advantage in NGoM (Aké-Castillo & Vázquez, 2008;Chakraborty & Lohrenz, 2015;Qian et al., 2003;Strom & Strom, 1996). More phytoplankton groups and possible predation avoidance mechanisms could be added to the model to give the bloom pattern (and subsequently the carbon export) more variance (Liszka, 2018;Rost & Riebesell, 2004). 5) Adding the higher trophic level biology could be the next step to improve the model. Marine fishes are reported to produce precipitated 535 carbonates within their intestines at high rates and contribute to TA increase in the top 1000 meters of ocean waters (Wilson et al., 2009). 6) This model did not include sediment silicate weathering and carbon flux through atmospheric deposition, which can potentially be important sources/sinks of carbon to the ocean waters as well (Jurado et al., 2008;Wallmann et al. 2008

Conclusions 540
Building upon a Global Climate & Earth system model, the high-resolution regional carbon model can reliably simulate the spatial and temporal pattern of the surface ocean carbon system in the GoM. In this study, we show for the first time a solid validation of a regional carbon model via direct comparison against high-frequency CO2 buoys, TA/DIC vertical profiles along the coastal transects, remote sensing-based ML model product, and underway pCO2 measurements (surface). We calculated the decadal trends of important carbon system variables such as pCO2, pH, air-sea CO2 exchange, and Ω over the NGoM and 545 Open GoM regions.
The GoM surface pCO2 values experience a steady increase from 2001 to 2019, with an increasing rate of 1.61 µatm yr -1 in NGoM and 1.66 µatm yr -1 in Open GoM, respectively. Correspondingly, the ocean surface pH is declining at a rate of 0.0020 yr -1 and 0.0015 yr -1 for NGoM and Open GoM, respectively. The surface Ω over the NGoM and Open GoM region remain 550 supersaturated with aragonite during the time span of the model but with a slightly decreasing trend. The carbon sink of both NGoM and Open GoM region have increasing trends and will continue to increase at a faster pace in the coming years under the prospect of climate change and the rising atmospheric pCO2.
We decouple the influence on surface pCO2 into thermal, nonthermal components and further analyze the surface pCO2 changes 555 due to gross primary production and air-sea CO2 flux. We find that the low temperature during winter and the biological uptake during spring and summer are the primary drivers making GoM an overall CO2 sink. During the modeled period of 2001-2009, the GoM overall is a CO2 sink with a mean flux rate of 0.62 mol C m -2 yr -1 (11.77 Tg C yr -1 ). The NGoM region is a CO2 sink year-round and is very susceptible to changes in river forcing. The Open GoM region is dominated by thermal effects and converts from carbon sink to a source during summer. 560 Historical simulation (His) and perturbed tests (Bry, NoR) are performed to determine whether observed changes in the GoM carbon system are driven by secondary effects of carbon accumulation from the global ocean or local forcing such as river inputs. The results show that, in addition to the increasing atmospheric pCO2 over the GoM, the spatial distribution and trend in carbon system variables could only be explained when the effects of carbon accumulation via boundary conditions and the 565 impact from river discharge are included. With a projected warming climate, we anticipate the GoM to be a stronger carbon sink due to an elevated river discharge and reduced impact from the global ocean.

Competing interests:
The authors declare that they have no conflict of interest. 575