Seasonal patterns of surface inorganic carbon system variables in the Gulf of Mexico inferred from a regional high-resolution ocean-biogeochemical model

. Uncertainties in carbon chemistry variability still remain large in the Gulf of Mexico (GoM), as data gaps limit our ability to infer basin-wide patterns. Here we configure and validate a regional high-resolution ocean-biogeochemical model for the GoM to describe seasonal patterns in surface pressure of CO 2 (pCO 2 ), aragonite saturation state ( Ω Ar ), and sea-air CO 2 flux. Model results indicate that seasonal changes in surface pCO 2 are strongly controlled by temperature across 15 most of the GoM basin, except in the vicinity of the Mississippi-Atchafalaya River System delta, where runoff largely controls dissolved inorganic carbon (DIC) and total alkalinity (TA) changes. Our model results also show that seasonal patterns of surface Ω Ar are driven by seasonal changes in DIC and TA, and reinforced by the seasonal changes in temperature. Simulated sea-air CO 2 fluxes are consistent with previous observation-based estimates that show CO 2 uptake during winter-spring, and CO 2 outgassing during summer-fall. Annually, our model indicates a basin-wide mean CO 2 uptake 20 of 0.35 mol m –2 yr –1 , and a northern GoM shelf (<200 m) uptake of 0.93 mol m –2 yr –1 . The observation and model-derived patterns of surface pCO 2 and CO 2 fluxes show good correspondence; thus this study contributes to improved constraints of the carbon budget in the region.


Introduction
The global ocean is absorbing approximately one third of the anthropogenic CO 2 released into the atmosphere from fossil fuel burning (e.g., Sabine et al., 2004;Gruber et al 2019), resulting in a sustained decline in seawater pH and the saturation and in situ observations are noted.For example, their model did not reproduce the decrease in surface pCO 2 linked to high primary production over the MARS mixing zone (Huang et al., 2015), and spatially averaged values of model pCO 2 were largely overestimated in the northern GoM during summer (by more than 100 µatm in several cases).In addition, the modeled sea-air CO 2 flux in the northern GoM (-0.32 mol m -2 yr -1 ) was about one third of the flux derived by Huang et al. (2015) and Lohrenz et al. (2018), while the modeled flux for the deep Gulf (-1.04 mol m -2 yr -1 ) was more than twice the flux derived by Robbins et al. (2014).In another modeling study, Laurent et al. (2017) examined near-bottom acidification driven by coastal eutrophication.Their model reproduced observed patterns in surface pCO 2 , sea-air CO 2 fluxes, pH, alkalinity, and DIC, but the model domain was limited to the Louisiana-Texas shelf.
Discrepancies between modeling results and observations, as well as the scarcity of biogeochemical modeling studies examining GoM-wide patterns, make additional modeling efforts necessary in order to reduce uncertainty in carbon patterns.In the present study, we use the outputs from a 15-component ocean-biogeochemical model for the GoM to characterize the seasonal variability of the inorganic carbon system variables at the ocean surface, with a focus on aragonite saturation state (Ω Ar ), pCO 2 , as well as sea-air CO 2 fluxes.This paper is structured such that we: 1) describe the ocean biogeochemical model and dataset used for the study; 2) validate the model based on observations from a coastal buoy, the GOMECC-1 cruise, and SOOP; 3) describe surface inorganic carbon system variables; 4) describe sea-air CO 2 fluxes in coastal and ocean domains; and 5) discuss the main model results in the context of previous observational and modeling studies.

Model
The biogeochemical model is similar to the one described by Gomez et al. (2018), but with an additional carbon module that simulates dissolved inorganic carbon (DIC) and total alkalinity (TA).The carbon module is based on Laurent et al. (2017) formulations, and considers a carbon to nitrogen ratio of 6.625 to link the carbon and nitrogen cycles.DIC is consumed by phytoplankton uptake, and produced by zooplankton excretion and organic matter remineralization, and affected by sea-air CO 2 fluxes.Changes in model TA are estimated using an explicit conservative expression for alkalinity (Wolf-Gladrow et al., 2007).Model CO 2 fluxes are derived using the Wanninkhof (2014) bulk flux equation.Details of the carbon module can be found in Section 1 in Supplement.A description of the model's nitrogen and silica cycle components is found in Gomez et al. (2018).
The coupled ocean circulation-biogeochemical model was implemented on the Regional Ocean Model System (ROMS; Shchepetkin and McWilliams, 2005).The model domain extends over the entire Gulf of Mexico (Fig. 1), with a horizontal resolution of ~8 km, and 37 sigma-coordinate (bathymetry-following) vertical levels.A third-order upstream scheme and a fourth-order Akima scheme were used for horizontal and vertical momentum, respectively.A multidimensional positive definitive advection transport algorithm (MPDATA) was used for tracer advection.Vertical turbulence was resolved by the Mellor and Yamada 2.5-level closure scheme.Initial and open-boundary conditions were derived from a 25 km resolution Modular Ocean Model for the Atlantic Ocean (Liu et al., 2015), which includes TOPAZ (Tracers of Ocean Phytoplankton with Allometric Zooplankton) as biogeochemical model (Dunne et al., 2010).The model was forced with surface fluxes of momentum, heat, and freshwater from the European Center for Medium Range Weather Forecast reanalysis product (ERA-Interim; Dee et al., 2011), as well as 54 river sources of freshwater, nutrients, TA, and DIC (http://waterdata.usgs.gov/nwis/qw,last accessed September 23 rd , 2018) (Aulenbach et al., 2007;He et al., 2011;Martinez-Lopez and Zavala-Hidalgo, 2009;Munoz-Salinas and Castillo, 2015;Stets et al., 2014).Monthly TA series for the MARS were derived from observations collected at the USGS stations 7373420 and 7381600.Following Stet and Striegl (2012), riverine DIC concentrations were calculated from observations of pH, TA, and temperature.Observational gaps in the Atchafalaya series were filled out using linear equations linking chemical properties at the Atchafalaya station to those at the Mississippi station (Section S2).For rivers other than the MARS, we used mean climatological DIC and TA values, as the availability of data for these rivers was insufficient to generate monthly series over the entire study period.The partial pressure of atmospheric CO 2 was prescribed as a continuous nonlinear function, derived from the Mauna Loa monthly CO 2 time series (www.esrl.noaa.gov/gmd/ccgg/trends/,last accessed August 16 th , 2018) using similar curve-fitting method that Thoning et al. (1989) (Section S3).
The ocean-biogeochemical model in Gomez et al. (2018) was spun-up for 40 years.In the present study, an additional 9-year spin-up for the carbon system components was completed, using the basin-model boundary conditions, ERA surface forcing, and river runoff from 1981-1983.After completing the spin-up, the model was run continuously from January 1981 to November 2014, with averaged outputs saved at a monthly frequency.DIC and TA, in conjunction with temperature and salinity, were used to derive the full set of inorganic carbon system variables, including pCO 2 and Ω Ar .The calculations were performed using the MatLab version of the CO2SYS program for CO 2 System Calculations (van Heuven et al., 2011), considering the total pH scale, the carbonic acid dissociation constants of Mehrbach et al. (1973) as refitted by Dickson and Millero (1987), the boric acid dissociation constant of Dickson (1990a), and the KSO 4 dissociation constant of Dickson (1990b).
For the present study, we focused on describing seasonal patterns in surface Ω Ar , surface pCO 2 , and sea-air CO 2 flux during 2005-2014 (i.e., the last 10 years of the model run).Ω Ar represents the degree of saturation of calcium carbonate (CaCO 3 ) phase aragonite, with Ω Ar values less than 1 indicating undersaturation (aragonite is thermodynamically unstable, which favors dissolution), and Ω Ar values greater than 1 indicating oversaturation (seawater favors aragonite precipitation).
Ω Ar is defined as: which is derived from the simulated DIC and TA, and K ' Ar is the apparent solubility product of the CaCO 3 phase aragonite in seawater, which increases with pressure and salinity, and decreases with temperature (Mucci, 1983;Millero, 1995).At a given pressure, temperature and salinity, changes on Ω Ar mainly depend on [CO 3 2-], and are positively related to changes in the TA to DIC ratio (Wang et al., 2013).

Data
Surface measurements of mole fraction of CO 2 (xCO 2 ), temperature, and salinity from the Central Gulf of Mexico Ocean Observing System (Coastal Mississippi Buoy) at 30°N and 88.6°W (Sutton et al., 2014;Fig. 1)

Model-data comparison
We used data from the Coastal Mississippi Buoy to evaluate the model's ability to reproduce coastal patterns in xCO 2 , temperature, and salinity in the northern GoM shelf (Fig. 2).Overall, simulated temporal surface patterns agreed with observations, especially considering that the buoy is located within a region highly impacted by river runoff, strong crossshore gradients, and high variability in salinity, DIC and TA.We can expect therefore that relatively small changes in river plume location (such as those derived from Mobile Bay and the Mississippi River) can significantly impact salinity and xCO 2 , making the exact reproduction of observed buoy patterns challenging.The best match between simulated and observed xCO 2 was during 2011-2012, where xCO 2 ranged from about 230 ppm in spring to more than 400 ppm in fall.
The pCO 2 GoM_2018 dataset was used to compare climatological seasonal patterns in pCO 2 (Fig. 3).Overall, simulated and observed pCO 2 patterns were in good agreement.In the open GoM region, there was a close match between model and observed patterns in July-December, with a relatively small model underestimation (~10 to 20 µatm) during February-June (Fig. 3a).In the northern GoM, the largest disagreement was observed in January-February (Fig. 3b), but this difference is most likely due to the reduced number of observations during winter in the pCO 2 GoM_2018 dataset (Fig. S6).
Indeed, January observations came from only one cruise, which largely increases observational uncertainty.A spatial visualization of the pCO 2 GoM_2018 observations and model outputs is presented for each calendar month in Fig. S6.The main spatial features were well reproduced by the model, including the pCO 2 minimum near the MARS region, and the large seasonal amplitude in the western Florida shelf.
We also compared vertical patterns in DIC, TA, temperature, and salinity derived from the model, with vertical profiles from the GOMECC-1 cruise (Fig. 4).The model reproduced the main patterns in DIC, TA, salinity, and temperature well, especially off Tampa.Monthly averaged model DIC and TA were underestimated in the upper 200 m off Louisiana, with the bias ranging from around 5 to 90 µmol kg -1 for DIC and 5 to 40 µmol kg -1 for TA, but the observations were within or close to the simulated variable's ranges during June-August 2000-2014.These model-observation differences could be partly due to misrepresentation of cross-shore transport in a region strongly influenced by the Mississippi river runoff.Also, TA and salinity were overestimated below 400 m at both stations by around 25 µmol kg -1 and 0.3, respectively, but this bias had a limited impact on the surface properties and fluxes examined (see following sections).Overall, our comparisons between model outputs and observations indicated that the model faithfully reproduced relevant inorganic carbon system features and patterns, and therefore was suitable for characterizing seasonal and spatial patterns of pCO 2 and Ω Ar for the 2005-2014 study period.

Surface pCO 2 and Ω Ar seasonality
Model derived patterns for surface pCO 2 showed significant seasonal variability across the GoM (Fig. 5).Minimum and maximum pCO 2 values were generally observed during winter and summer seasons, respectively, although large spatial differences were observed among the shelf regions.A notable model feature was observed in the central part of the northern GoM near the MARS delta, where pCO 2 displayed low values year-round (<350 µatm), with a seasonal minimum in spring.
Other coastal regions less impacted by riverine discharge displayed much higher pCO 2 values during spring and summer (Fig. 5b,c).The continental shelf with the maximum seasonally averaged pCO 2 was the west Florida shelf, where pCO 2 reached values greater than 450 µatm during the summer.Seasonality in modeled pCO 2 was strongly modulated by SST, such that the annual amplitude for these two variables displayed very consistent spatial patterns (Fig. 6a,b; Fig. S7).The greatest annual signal for pCO 2 and SST was within the northern GoM shelf and west Florida shelf, and the smallest was in the Loop Current region.Monthly time series of modeled pCO 2 and SST were strongly correlated in all regions except near the MARS delta (Fig. 6c).
The low pCO 2 -SST correlation near the MARS delta can be explained by the role that river runoff and enhanced primary production play as drivers of carbon system variability.This was evident in the variability of modeled pCO 2 along the salinity gradient linked to the Mississippi river plume (Fig. 7).The simulated surface pCO 2 patterns during spring and summer displayed a marked increase from mid to low salinities (Fig. 7a,d), which was also associated with an increase in DIC (Fig. 7b,e).The minimum pCO 2 values were about 285 µatm in spring and 320 µatm in summer, at salinities close to 30 and 27, respectively.To identify the drivers of DIC variability along the salinity gradient, we displayed the simulated budget terms for surface DIC as a function of salinity.These budget terms correspond to the sea-air CO 2 flux (Sea-air), the combined effect of advection and mixing (Adv+Mix), and the net community production (NCP), the latter representing the difference between primary production and respiration (i.e.biologically driven changes in DIC).The derived patterns for spring-summer showed model DIC losses at mid salinities mainly driven by NCP, indicative of a biologically induced drawdown of pCO 2 .During fall (Fig. 7g-i), as well as winter (not shown), NCP was much smaller than during springsummer, and DIC was mainly controlled by sea-air exchange and advection plus mixing processes.As a consequence, model surface pCO 2 did not show a mid salinity minimum linked to phytoplankton uptake.
The simulated patterns for surface Ω Ar (Fig. 8) revealed a significant meridional gradient from fall to spring, with minimum values in the inner shelves from northern GoM and west Florida (2.5-3.6), and maximum values over the Loop Current and west of the Yucatan Peninsula (3.9-4.1).During summer, the simulated surface Ω Ar reached its maximum near the MARS delta (>4.5), while relatively weak Ω Ar gradients were observed across the open GoM region.Surface Ω Ar generally displayed maximum values in summer and minimum in winter, though always well above the saturation threshold of 1.This seasonal variation in surface Ω Ar was strongly correlated to changes in the TA:DIC ratio and SST (Fig. 9a,b).
Although the seasonal patterns for Ω Ar and pCO 2 displayed a similar phase (maximum in summer, minimum in winter), the spatial variability of these two variables was opposite.This was most evident during spring-summer (Figs. and TA (r = 0.65 and 0.60, respectively, with wind leading by 1 month; Fig. S12a), influenced this seasonal pattern.The similar annual amplitude and phase for DIC and TA, as well as high TA values year-round, caused a relatively weak seasonal variability for pCO 2_at25 and Ω Ar_at25 on the Yucatan shelf.Still, a significant correlation between easterly winds and surface pCO 2_at25 (r = 0.55) was found in the northern Yucatan coast, with pCO 2_at25 usually peaking during spring (Fig. S12b).

Sea-air CO 2 fluxes
Seasonal changes in surface model pCO 2 , mainly driven by SST changes (Fig 6c), determined strong seasonal variability in simulated sea-air CO 2 fluxes.As a consequence, the GoM becomes a CO 2 sink in winter-spring and a CO 2 source in summer-fall (Fig. 11a-d).An exception to this pattern occurred close to the MARS delta, which is predominantly a CO 2 sink year-round.In this region, the pCO 2 drop induced by phytoplankton uptake during spring-summer (Fig. 7a,d) determined maximum uptake of atmospheric CO 2 at mid salinities (seen in the sea-air exchange term in Fig. 7c,f).The greatest model CO 2 uptake, above 7 mmol m -2 d -1 , occurred over the northern GoM shelf during winter, as this region experiences the lowest surface pCO 2 values induced by the coldest winter conditions in the region (Fig. S7).The greatest model CO 2 outgassing, disregarding local peaks near major river mouths like the Mississippi river, was observed on the west Florida shelf (northern inner shelf in particular), southern Texas shelf (northern-west GoM), and western Yucatan Peninsula during the summer, ranging from ~2 to 3 mmol m -2 d -1 (Fig. 11c).Maximum SST values characterized summer conditions in these regions (Fig. S7).The annual mean pattern showed modeled CO 2 uptake ranging from -4 to -1 mmol m -2 d -1 in the northern GoM, and from -2 to 0 mmol m -2 d -1 elsewhere (Fig. 11e).In addition, the pattern revealed areas where CO 2 outgassing occurred near the Mississippi River, Atchafalaya River, and Mobile Bay mouths, on the western Yucatan Peninsula, and nearshore over the west Florida shelf (Fig. 11e).
The estimated monthly patterns for modeled sea-air CO 2 flux revealed prevailing CO 2 outgassing during May-October in west Florida, western GoM, and Yucatan Peninsula, and June-October in the northern and open GoM (Fig. 11f).
The timing for the maximum CO 2 outgassing was June-July in the western GoM, August in west Florida and Yucatan, and

Simulated carbon patterns
Characterization of historical carbon system patterns are needed to advance our understanding of carbon dynamics, as well as to identify coastal ecosystem susceptibility to ocean acidification (Wanninkhof et al., 2015).Previous studies have described to some degree surface pCO 2 seasonality within the GoM (e.g.Lohrenz et al., 2010;2018;Robbins et al., 2018), but less has been done to describe seasonal patterns for other inorganic carbon system variables.In the present study, we focused our analysis on the seasonal cycles of surface pCO 2 and Ω Ar , but seasonal patterns of surface DIC and TA were also reported.We used a similar model to the one configured by Gomez et al. (2018) for the GoM, with an extra carbon module to simulate carbon dynamics, following model formulations described by Laurent et al. (2017).As shown in Section 3, the model simulated the main surface spatio-temporal patterns for the inorganic carbon system well.Compared to a previous basin-wide modeling effort (Xue et al., 2016), our model shows significantly less seasonal biases in surface pCO 2 , with relatively minor pCO 2 underestimation during spring (<20 µatm).Further model refinements could be required for improving the representation of carbon system dynamics.These include incorporating additional model components and processes, like dissolution and precipitation of calcium carbonate that will affect  (2015), which showed the most buffered surface waters off the MARS delta during summer.We found a strong positive correlation between the TA:DIC ratio and Ω Ar , which reflects the Ω Ar dependency to changes in [CO 3 2-].This is consistent with Wang et al. (2013), who reported spatial covariation of these two variables over the GoM and the eastern coast of USA.
We also found a strong positive correlation between SST and Ω Ar , which can be linked to the impact of temperature on aragonite solubility (aragonite solubility decreases with temperature) and sea-air CO 2 fluxes (warm conditions favor surface DIC decrease due to CO 2 outgassing, which increases the TA:DIC ratio).Comparison between monthly climatologies for surface Ω Ar and Ω Ar_at25 reveals that Ω Ar seasonality induced by changes in the TA:DIC ratio tends to be reinforced by temperature-induced changes.
Surface Ω Ar patterns can be useful to identify regions more vulnerable to ecosystem disturbances induced by surface ocean acidification.Our model indicates minimum surface Ω Ar ranging from 2.5 to 3.4 on the northern GoM and west Florida inner shelves during winter, and greater than 3.4 on the western GoM and Yucatan shelves.This suggests higher ecosystem resilience to surface ocean acidification in the latter regions.Surface Ω Ar patterns do not necessarily reflect vulnerability of coastal benthic organisms to ocean acidification, since Ω Ar values for surface and bottom layers can largely differ in regions where the water column is strongly stratified.This is the case for the Louisiana inner shelf during summer, which displayed maximum surface Ω Ar values (>4.2) linked to high biological uptake, but low bottom Ω Ar values (<2.6; not shown) due to bottom acidification induced by organic carbon remineralization and weak bottom ventilation (see Cai et al. (2011) and Laurent et al. (2017) for further discussion).However, our model outputs did not reveal such signature of bottom acidification on the west Florida, western GoM and Yucatan shelves, as these regions display relatively weak vertical stratification and lower eutrophication levels compared to the northern GoM shelf.
Sea-air CO 2 flux derived from the model output shows that the GoM is a CO 2 sink during winter-spring, and a CO 2 source during summer-fall.However, significant differences in the annual flux magnitude were observed among regions, which can be associated with distinct ocean-biogeochemical regimes.The northern GoM shelf, a river-dominated ocean margin strongly influenced by seasonal patterns in MARS runoff (McKee et al., 2004;Cai et al., 2013), is the coastal region with the lowest surface pCO 2 and the largest CO 2 uptake from the model.This pattern is due to the substantial cooling experienced by the northern GoM shelf during winter (linked to its northernmost location), and the enhanced biological uptake promoted by river runoff near the MARS delta during spring-summer.Our results support the framework proposed by Huang et al. (2015) for the Mississippi river plume during spring-summer, which indicates i) high pCO 2 levels and CO 2 outgassing at low salinities (<20), linked to the low productivity, high turbidity, and CO 2 oversaturated waters delivered by the Mississippi river; ii) minimum pCO 2 values and maximum atmospheric CO 2 uptake at mid salinities (20-33), as high phytoplankton production, induced by decreased water's turbidity and nutrient runoff, produces a drop in surface DIC, and iii) increased pCO 2 levels and sea-air CO 2 flux at high salinities (>33), as phytoplankton production declines offshore in the oligotrophic open GoM waters.In the west Florida and western GoM shelves, two coastal margins that are not strongly influenced by river runoff, temperature plays a dominant role as driver of pCO 2 and sea-air CO 2 flux seasonality.As a result, the annually integrated sea-air CO 2 flux (per m 2 ) in these two shelves represents only 31% and 23% of the simulated carbon uptake in the northern GoM, respectively.In the Yucatan Peninsula, temperature is likewise the main driver of model surface pCO 2 and CO 2 flux seasonality.The zero flux in this region results from a less pronounced winter cooling, which determines a relatively weak carbon uptake during winter-spring.However, wind-driven upwelling also plays a role by increasing model surface pCO 2 during spring, especially nearshore.Although previous studies have documented the impact of coastal upwelling on SST and surface chlorophyll in the Yucatan shelf (e.g.Zavala-Hidalgo et al., 2006), no study has addressed the associated impact on carbon chemistry, as insufficient inorganic carbon observations exist for this region.Further observational studies are required therefore to corroborate this dynamic.Finally, the simulated annual carbon uptake was weak for most of the GoM basin.Therefore, it is likely that relatively small disturbances in the pCO 2 drivers could turn the carbon sink regions into carbon sources.A potential mechanism for this change is ocean warming, since future ocean projections in the GoM suggest a significant SST increase (>2°C) due to anthropogenic climate change to the end of the twenty-first century (Liu et al., 2012;2015;Alexander et al., 2020;Shin and Alexander, 2020).This is a topic deserving examination for future modeling efforts.

CO 2 flux comparison
Table 2 shows mean CO 2 fluxes derived from our model, previous regional studies for the GoM, and global datasets.The (2014), and Bourgeois et al. (2016).Annual CO 2 fluxes for the GoM basin displayed a significant dispersion, ranging from -0.72 to +0.20 mol m -2 yr -1 .However, the three regional studies providing basin-wide estimates (including ours) agree that the GoM is a carbon sink.We obtained an average value of -0.35 mol m -2 yr -1 , which is comparable with Robbins et al.
(2014) and Xue et al. (2016) estimates.In contrast, two out of three basin fluxes derived from global gridded datasets, Takahashi et al. (2009) andLandshützer et al. (2016), suggest that the GoM is a weak CO 2 source.This discrepancy between regional and global studies most likely reflects inaccuracy in global datasets, due to the low density of pCO 2 observations in the GoM basin and coarse grid resolutions (5° latitude x 4° longitude in Takahashi et al. 2009 and1° latitude x 1° longitude in Landshützer et al. 2016).
We obtained fluxes that are in reasonable agreement with observation-based fluxes for most of the sub-regions depicted in Figure 1.In the open GoM region, our mean flux (-0.33 mol m -2 yr -1 ) is about 70% of the flux derived by Robbins et al. (2014).For all four GoM shelf regions combined (west Florida, northern GoM, western GoM, and Yucatan), our estimated flux (-0.39 mol m -2 yr -1 ) is 20% above the value reported by Laruelle et al. (2014).In the northern GoM, our simulated flux (-0.93 mol m -2 yr -1 ) is remarkably similar to the reported fluxes of Huang et al. (2015) and Lohrenz et al.
(2018) (-0.95 and -1.1 mol m -2 yr -1 , respectively).In the Yucatan Peninsula, our zero flux condition is close to the weak uptake condition derived by Robbins et al. ( 2014) (-0.09 mol m -2 yr -1 ).The major disagreement between our estimates and previous studies is on the west Florida and western GoM shelves.We determined that these two regions are carbon sinks (-0.30 and -0.22 mol m -2 yr -1 , respectively), whereas observational studies by Robbins et al. (2014;2018), as well as the modeling study by Xue et al. (2016), estimated a mean carbon outgassing condition.Some overestimation in our modeled CO 2 uptake is possible, as the model surface pCO 2 in the open GoM tended to be underestimated during late winter and spring.However, the observational uncertainty in Robins et al. (2014;2018) also needs to be considered.The dataset of underway pCO 2 measurements, used to generate the observed bulk CO 2 fluxes, has very limited spatial coverage over the western GoM.Also, this dataset has a reduced number of winter observations in west Florida and other GoM regions (only 8% of the GoM data were collected in December-February, less than 2% during January).A correct estimation of the winter flux is important, as this season largely determines the sign of the annual flux.Indeed, excluding winter, our simulated spring to fall flux for west Florida is positive (+0.12 mol m -2 yr -1 ).
The simulated fluxes largely differ from the fluxes reported by Xue et al. (2016), which was the only previous regional modeling study describing basin wide patterns in the GoM.They obtained a three times stronger uptake in the open GoM, and much weaker uptake on the shelf regions (e.g.their simulated annual flux for the northern GoM shelf was one third of our estimation).We believe these differences in CO 2 fluxes can be mainly explained by pCO 2 biases in the model used in Xue et al. (2016).Indeed, their model underestimated surface pCO 2 in the open GoM, and thus obtained a marked pCO 2 minimum over the Loop Current region (see their Fig.13a), a pattern not supported by SOOP observations (Fig. S6).
In addition, their model largely overestimated surface pCO 2 on the northern GoM and west Florida inner shelves, especially during summer-fall, not reproducing well the marked pCO 2 drop that is observed close to the MARS delta.

Summary and Conclusions
We configured a coupled ocean-biogeochemical model to examine inorganic carbon chemistry patterns in the GoM.The model was validated against observations from a coastal buoy, research cruises, and ships of opportunity, showing smaller seasonal and regional bias for surface pCO 2 than previous modeling efforts in the region.We described seasonal patterns in surface pCO 2 and Ω Ar .Both variables show maximum values during late summer and minimum during winter and early spring.The seasonal cycle for pCO 2 is strongly controlled by temperature, while Ω Ar follows changes in the TA:DIC ratio and temperature.Model results also indicated that river runoff and wind-driven circulation significantly influence coastal DIC and TA patterns in coastal regions, impacting Ω Ar , pCO 2 , and sea-air CO 2 flux seasonality.Simulated fluxes show CO 2 uptake prevailing during winter-spring, and CO 2 outgassing during summer-fall.The integrated annual flux for the GoM basin is -0.35 mol m -2 yr -1 (-4.2 g C m -2 yr -1 ).The largest model CO 2 uptake is in the northern GoM shelf, linked to the most intense winter cooling, and significant biological uptake during spring-summer.The weakest CO 2 uptake is in the Yucatan Peninsula, mainly a consequence of the relatively warm conditions experienced by this region during winter-spring, and to a less degree wind-driven upwelling of DIC-rich waters.Sub-regional estimates are in general consistent or close to previous observational studies, with the exception of the west Florida and western GoM shelves.We suggest that part of these discrepancies could be related to the still reduced spatio-temporal coverage in the underway pCO 2 measurement dataset over those two regions, especially during wintertime.
5b,c; 8b,c), when the highest Ω Ar and lowest pCO 2 values were co-located near the MARS delta, and the lowest Ω Ar and highest pCO 2 values were in the west Florida and northern-west GoM shelves.The annual amplitude of Ω Ar displayed a similar pattern to the annual amplitude of surface salinity, especially over the northern GoM, indicating a strong influence of river discharge on Ω Ar seasonality(Figs.S8 and S9).The correlation between Ω Ar and salinity showed negative values over the northern GoM and eastern part of the open GoM.This pattern was consistent with enhanced biological uptake of DIC promoted by MARS's nutrient inputs (Fig.9c).To better describe the impact of SST in the simulated pCO 2 and Ω Ar variability, we calculated average monthly climatologies for temperature-normalized pCO 2 and Ω Ar at 25°C (pCO 2_at25 and Ω Ar_at25 , respectively), and compared them with non-normalized patterns in five regions designated as the northern GoM shelf, west Florida shelf, western GoM shelf, Yucatan shelf, and open GoM (Fig.10a-d; regions depicted in Fig.1).Surface pCO 2_at25 and Ω Ar_at25 were calculated with the CO2SYS program, using the simulated DIC, TA, and salinity patterns, and 25°C (which is close to the average SST over the GoM basin).The strong influence of SST on model pCO 2 was evident when we compared the monthly climatologies for pCO 2 and pCO 2_at25 (Fig.10a,b).Surface pCO 2_at25 displayed much weaker annual variation than surface pCO 2 , and the timing for the seasonal maxima and minima largely differed.Indeed, surface pCO 2_at25 peaked during January-February in the northern GoM, during March in the west Florida and western GoM regions, and during February in the open GoM regions, i.e. when pCO 2 was at or near its lowest levels.The comparison between Ω Ar and Ω Ar_at25 also revealed significant temperature influence on model Ω Ar seasonality (Fig.10c,d).Specifically, SST amplified the annual variation in Ω Ar , while having a relatively weak impact on the Ω Ar seasonal phase.Both Ω Ar and Ω Ar_at25 were inversely related to pCO 2_at25 , reflecting the variables dependency to DIC and TA (Ω Ar increases with TA and decreases with DIC, while pCO 2_at25 has the opposite pattern).Simulated climatological patterns for DIC and TA (Fig.10e,f; Figs.S10 and S11) allowed us to investigate the importance of DIC and TA as drivers of pCO 2_at25 and Ω Ar_at25 seasonality.In the open GoM, west Florida, and western GoM regions, changes in TA were small, so the seasonal pattern in Ω Ar was mainly due to DIC changes.Maximum surface DIC values during late winter and early spring can be linked to increased uptake of atmospheric CO 2 (see Section 5) and enhanced vertical mixing, promoted by surface cooling and winds.Alternatively, both DIC and TA played an important role modulating Ω Ar seasonality in northern GoM and Yucatan Peninsula shelves.In the former, the annual variation of DIC and TA was strongly modulated by river runoff, which is mostly associated with the MARS.Whether the MARS dilutes ocean DIC and TA depends on the season.Alkalinity in the Atchafalaya river was lower than the open GoM alkalinity year-round, whereas Mississippi alkalinity was lower than open GoM alkalinity during December-June and greater the rest of the year (Fig.S3a).The DIC of the Atchafalaya was smaller than open GoM DIC during December-May and greater from June to November, while Mississippi DIC was greater or equal to the open GoM DIC year-round (Fig.S3b).We did not prescribe time-evolving DIC and TA for rivers other than the Mississippi River, but according to USGS records most of these other rivers have long-term average DIC and TA smaller than the oceanic values.Consequently, low TA values in the northern GoM during spring can be explained by a dilution effect, linked to maximum river discharge in the northern GoM during winter-spring.Low DIC values during spring-summer can be associated with high biological uptake, promoted by riverine nutrients and enhanced solar radiation, along with dilution (especially in spring) linked to high discharge of low DIC waters delivered by major river inputs, like the Atchafalaya River and Mobile Bay.This is not the case for the Mississippi river, which had DIC values greater than the open GoM.Along the Yucatan Peninsula, simulated surface DIC and TA patterns showed maximum values in summer and minimum values in winter.Coastal upwelling of DIC and TA-rich waters along the northern Yucatan Peninsula coast, reflected in a significant correlation between easterly (alongshore) winds and both DIC September in the northern and open GoM.The timing for the maximum CO 2 uptake was January in the northern GoM, west Florida, and Yucatan Peninsula, and February in the western and open GoM.The model annual flux for the northern GoM, west Florida, western GoM, Yucatan, and open GoM are -2.56,-0.81, -0.60, 0.0, and -0.90 mmol m -2 d -1 , respectively.For the entire GoM basin, the simulated average annual flux and standard deviation was -0.97 and 2.78 mmol m -2 d -1 (-0.35 and 1.01 mol m -2 yr -1 ), respectively.Integrated across the entire model domain, the resulting flux was -7.0 Tg C yr -1 .

Figure 1 .Figure 2 .Figure 3 .
Figure 1.Model snapshot of surface dissolved inorganic carbon (mmol m -3 ) during May 1 st of 2009.Regions used to describe model results are the western GoM shelf, the northern GoM shelf, the west Florida shelf, the Yucatan shelf, and 575 open GoM.Shelf regions are delimited offshore by the 200 m isobath.Black stars depict the location of two GOMECC stations at the Mississippi (M) and Tampa (T) lines used to validate the model.Red star depicts the location of the Coastal Mississippi Buoy (CMB).Blue circles indicate USGS stations 7373420 and 7381600 at the Mississippi (MS) and Atchafalaya (AT) rivers, respectively.Magenta polygon demarks the region near the Mississippi delta used to derive patterns in Fig. 7. 580

Figure 4 .Figure 5 .
Figure 4. Comparison between profiles of dissolved inorganic carbon (DIC), total alkalinity (TA), salinity, and temperature from monthly model outputs (blue lines) and GOMECC-1 data (red dots) for the most oceanic station on the (a) Tampa and (b) Mississippi lines (station locations are shown in Figure 1 as black stars).The model's variables range for June-August during 2000-2014 is also shown as cyan shade.

Figure 6 .Figure 7 .
Figure 6.(a,b) Seasonal amplitude patterns for model surface pCO2 and SST.The seasonal amplitude is the difference between the maximum and minimum values from monthly climatologies at each grid point (c) Correlation between surface model pCO 2 and SST.Black contour depicts the 200 m isobath.610
Wang et al. (2013), and Wanninkhof et al.soscale dynamics.Our current model configuration represents an important advance in the model capabilities for the GoM, capturing realistically dominant seasonal patterns.Simulated patterns in surface pCO 2 across the GoM show maximum values in spring-summer and minimum in winter, with seasonally averaged values ranging from around 250 to 500 µatm.Seasonal variability in SST was the main driver of surface pCO 2 seasonality across the GoM, except for the region around the MARS delta, where river runoff and biological uptake of DIC played a significant role during spring-summer.The pCO 2 -SST correlation pattern derived from the model is consistent with previous observational studies, which suggested an increased correlation between pCO 2 and SST away from the Mississippi-Atchafalaya mixing zone, in open GoM waters (e.g.Lohrenz et al., 2018).Simulated patterns in surface Ω Ar showed maximum values in late summer and minimum in late winter, with most values ranging from 3 to 4.4 units.The meridional and cross-shore gradients for model surface Ω Ar are consistent with patterns observed byGledhill et al. (2008).Our model results also agree with observations byGuo et al. (2012),Wang et al. (2013), and Wanninkhof et al.
TA, improving the representation of landocean biogeochemical fluxes (e.g.prescribing time evolving TA and DIC for rivers other than the MARS), and increasing the model'

Table 1 .
Mean CO 2 flux derived from monthly model outputs during 2005-2014.Standard deviation is in parenthesis.Negative flux implies ocean CO 2 uptake, and positive flux CO 2 outgassing (shown in red).Shelf regions are depicted in Fig.

Table 2 .
Comparison between annual sea-air CO 2 fluxes (mol m -2 yr -1 ) derived from our model results and previous studies in the Gulf of Mexico.Standard deviation is in parenthesis.Negative flux implies ocean CO 2 uptake, and positive flux CO 2 outgassing (shown in red).Shelf regions are depicted in Fig.1.