Modeling p CO 2 Variability in the Gulf of Mexico

22 A three-dimensional coupled physical-biogeochemical model was used to 23 simulate and quantify temporal and spatial variability of sea surface p CO 2 in the Gulf of 24 Mexico (GoM). The model was driven by realistic atmospheric forcing, open boundary 25 conditions from a data-assimilative global ocean circulation model, and observed 26 freshwater and terrestrial nutrient and carbon input from major rivers. A seven-year 27 model hindcast (2004–2010) was performed and validated against ship measurements. 28 Model results revealed clear seasonality in surface p CO 2 and were used to compute 29 carbon budgets in the Gulf. On average, the GoM was found to be a CO 2 sink with a flux 30 of 1.11±0.84 × 10 12 mol C yr -1 , which, together with the enormous fluvial carbon input, 31 was balanced by the carbon export through the Loop Current. Two model sensitivity 32 experiments were performed: one without biological sources and sinks and the other 33 using river input from the 1904-1910 period as simulated by the Dynamic Terrestrial 34 Ecosystem Model (DLEM). It was found that biological uptake was the primary driver 35 making GoM an overall CO 2 sink and that the sub-regional carbon budget was 36 susceptible to changes in river forcing. When the 1904-1910 river conditions were 37 applied, the northern GoM became a CO 2 source instead.


Introduction
Human consumption of fossil fuels has resulted in continuously increasing levels of atmospheric CO 2 since the Industrial Revolution began around 1750 (Solomon et al., 2007).If the increasing trend continues, the projected pCO 2 by the end of the 21st century (970 ppm in A1F1 scenario; Stocker et al., 2014) could be nearly triple the present level.In the face of different climate scenarios, a better understanding of the oceans' role in regulating the global carbon cycle is crucial, because oceans act not only as receivers of the enormous carbon loading from coastal rivers (Cai et al., 2011a;Bauer et al., 2013) but also as vast carbon reservoirs via the "carbon pump" mechanism (Sabine et al., 2004;Sabine and Tanhua, 2010).On regional scales, the marine carbon cycle tends to be more complicated and shows contrasting behaviors in different areas (coastal vs. open ocean, low latitude vs. high latitude, etc.) and during different seasons (e.g., Lohrenz et al., 2010, for the northern Gulf of Mexico (GoM); Jiang et al., 2008, for the South Atlantic Bight;Signorini et al., 2013, for the North American east coast; Tsunogai et al., 1999, for the East China Sea).Quantifying the ocean carbon budget is therefore a difficult task.Coupled physical and biological models are useful tools for understanding complex biogeochemical processes and estimating carbon and nutrient fluxes in coastal oceans where spatial and temporal heterogeneities are high and data are sparse (e.g.Fennel and Wilkin, 2009;Fennel, 2010;Fennel et al., 2011;and He et al., 2011).
Our study focuses on the carbon cycle in the GoM.One unique feature of the Gulf is that it receives enormous riverine nutrient and carbon inputs, both organic and inorganic, the majority of which are from the Mississippi-Atchafalaya River system.Excessive nutrient loading causes coastal eutrophication, which triggers not only the well-known hypoxia phenomenon (a.k.a. the "dead zone", Rabalais et al., 2002) but also a newly revealed coastal ocean acidification problem (Cai et al., 2011b).However, the carbon cycling associated with such enormous terrestrial carbon and nutrient inputs remains unclear: on the one hand extensive riverine carbon input results in CO 2 over-saturation in coastal waters, which serve as a CO 2 source to the atmosphere (e.g.Lohrenz et al., 2010;Guo et al., 2012); on the other hand, enhanced primary production in the river plume due to significant inputs of inorganic nutrients induces a net influx of CO 2 although the Mississippi River plume region is an overall heterotrophic system that breaks down organic carbon (Murrell et al., 2013;Huang et al., 2013Huang et al., , 2015)).Further offshore, the circulation in the GoM is largely influenced by the energetic Loop Current.Large anticyclonic eddies aperiodically pinch off from the Loop Current (Sturges and Leben, 2000), which, along with the wind-driven cross-shelf circulation and other mesoscale and submesoscale processes, enhance material exchanges between the eutrophic coastal waters and oligotrophic deepocean waters (e.g., Toner et al., 2003).Indeed, a recent observational study suggested a significant dissolved inorganic carbon export (DIC, ∼ 3.30 × 10 12 mol C yr −1 ) from the GoM shelves to the Loop Current waters (Wang et al., 2013).
While global inorganic carbon budgets have been made available through joint seawater CO 2 observations (e.g.World Ocean Circulation Experiment and Joint Global Ocean Flux study; Sabine et al., 2004;Feely et al., 2004;Orr et al., 2005), they are too coarse to represent CO 2 variability in the GoM (Gledhill et al., 2008).Other recent efforts were able to provide GoM subregional carbon assessments based on limited in situ observations (e.g.Cai et al., 2003, Lohrenz et al., 2010, and Huang et al., 2013, 2015, focused on the Mississippi River plume and the Louisiana Shelf; Wang et al., 2013, covered three cross-shelf transects in the northeastern GoM but only for one summer).Significant uncertainties exist in such budget estimations due to large temporal and spatial gaps presented in the observations (e.g.Coble et al., 2010;Hofmann et al., 2011;Robbins et al., 2014).In this regard, coupled physical-biogeochemical models are capable of representing the biogeochemical cycle with realistic physical settings (e.g., ocean mixing and advection) and providing an alternative means for a Gulf-wide carbon budget estimation.
Here we present a GoM pCO 2 analysis based on the results of a coupled physical-biogeochemical model simula-tion.Our objective was to simulate the CO 2 flux at the air-sea interface (which at present is based on observational analyses alone and subject to large uncertainty), as well as its variability in relationship with river plume dynamics and dominant oceanic processes in different regions of the GoM.

Method
Our analysis uses solutions from a coupled physicalbiogeochemical model covering the GoM and South Atlantic Bight waters (Xue et al., 2013, model domain see Fig. 1).The circulation component of the coupled model is the Regional Ocean Modeling System (ROMS; Haidvogel et al., 2008;Shchepetkin and McWilliams, 2005;Hyun and He, 2010) and is coupled with the biogeochemical module described in Fennel et al. (2006Fennel et al. ( , 2008Fennel et al. ( , 2011)).The nitrogen cycling parameterization has seven state variables: two species of dissolved inorganic nitrogen (DIN hereafter: nitrate (NO 3 ) and ammonium (NH 4 )), one functional phytoplankton group, chlorophyll as a separate state variable to allow for photoacclimation, one functional zooplankton group, and two pools of detritus representing large, fast-sinking particles and suspended, small particles.The carbon cycle is connected to the nitrogen cycle via a C-to-N ratio of 6.625 for the organic components (phytoplankton, zooplankton, large and small detritus).The sediment component of the biogeochemical model is a simplified representation of benthic remineralization processes, where the flux of sinking organic matter out of the bottommost grid box results immediately in a corresponding influx of ammonium and DIC at the sediment / water interface.The parameterization accounts for the loss of fixed nitrogen through sediment denitrification based on the linear relationship between sediment oxygen consumption and denitrification reported by Seitzinger and Giblin (1996) and only accounts for the portion of denitrification that is supported by nitrification of ammonium in the sediment (referred to as coupled nitrification / denitrification).
A 7-year (1 January 2004-31 December 2010) model hindcast was performed, driven by NCEP's high-resolution combined model and assimilated atmospheric dataset (North American Regional Reanalysis, www.cdc.noaa.gov),open boundary conditions for ocean model (temperature, salinity, water level, and velocity) from a data-assimilative global ocean circulation model (HYCOM/NCODA, Chassignet et al., 2007), and observed freshwater and terrestrial nutrient input from 63 major rivers (Aulenbach et al., 2007;Milliman and Farnsworth, 2011;Fuentes-Yaco et al., 2001;and Nixon, 1996).Model validations (physics, nutrients, and chlorophyll) and a nitrogen budget were reported in Xue et al. (2013Xue et al. ( , 2015)).In this study, we have focused on the carbon cycle in the GoM.As in Xue et al. (2013), we considered the first year of the simulation (2004) as model spin-up; all results presented here use model output from 2005 to 2010.The carbonate chemistry of the coupled model is based on the standard defined by the Ocean Carbon Cycle Model Intercomparison Project Phase 2 (Orr et al., 1999).There are two active tracers, DIC and alkalinity, to determine the other four variables of the carbonate system (i.e.pCO 2 , carbonate ion concentration, bicarbonate ion concentration, and pH; Zeebe and Wolf-Gladrow, 2001).Details of the formulas used in the simulation are provided in the Appendix A.
Similar to the results reported by Hofmann et al. (2011), we found that the model-simulated DIC concentration in the water column was very sensitive to the initial conditions.Although there were many historical measurements in the GoM, these data were limited to the northern GoM shelf regions and thus were insufficient to initialize the model.Instead, we tested model sensitivity using three sets of initial and open boundary conditions, which were derived using the empirical salinity-temperature-DIC-alkalinity relationships described in Lee et al. (2000Lee et al. ( , 2006)); Cai et al. (2011a), and Wang et al. (2013), respectively.Among them, the initial condition prescribed following Lee et al. (2000, 2006, Fig. 2; for details see Appendix B) provided the best model-data comparison.For the open boundary condition, we found simulated surface pCO 2 exhibited very limited variance (< 5 %) regardless of which conditions were applied.To be consistent with the setup of the initial condition, the results presented here were driven by boundary conditions derived from Lee et al. (2000Lee et al. ( , 2006)).For particular organic carbon, we set a small, positive value for both phytoplankton and zooplankton along the open boundaries.
The carbon cycle parameterizations used in this study followed the same approach and values as in Fennel et al. (2008), Fennel andWilkin (2009), andFennel (2010).For gas exchange calculations we followed the formulas in Wanninkhof (1992; for details see Appendix C).For air pCO 2 , we utilized the Atmospheric Infrared Sounder (AIRS Science Team, 2008) monthly gridded observation dataset and averaged them over the study area.We applied the curve-fitting method using a C language program named CCGCRV (http://www.esrl.noaa.gov/gmd/ccgg/mbl/crvfit/crvfit.html,Fig. 3), and the air pCO 2 in the gas exchange calculation was prescribed as  Lee et al. (2000Lee et al. ( , 2006)).
To account for riverine inputs, we constructed climatological monthly alkalinity time series by averaging all available US Geological Survey (USGS) observations for each major river, including the Mississippi, Atchafalaya, Mobile, and Brazos in the GoM.Because direct riverine DIC measurements were not available, we approximated riverine DIC inputs using the corresponding alkalinity value plus 50, following the observational study by Guo et al. (2012).The fluvial DIC input to the GoM was estimated as ∼ 2.18 × 10 12 mol C yr −1 , the majority of which was delivered by the Mississippi-Atchafalaya River (∼ 1.80 × 10 12 mol C yr −1 , Fig. 4, comparable with the estimation in Cai et al., 2003).
The results of three model experiments covering the period of 2004-2010 are presented in this study, in which, Experiment 1 (Exp1) was a "control run", with observed riverine inputs from USGS and biological sources and sinks of DIC and alkalinity in the water column; Experiment 2 (Exp2) was a "no-biology run", where all biological sources and sinks of DIC and alkalinity were disabled, similar to the experiment described in Fennel and Wilkin (2009); and Experiment 3 (Exp3) had the same set up as Exp1, but the riverine inputs (water, nutrients, and carbon of the Mississippi-Atchafalaya River) were taken from the Dynamic Land Ecosystem Model (DLEM; Tian et al., 2015) simulation for the period of 1904-1910 (Fig. 4).Specifically, we used the monthly model outputs of water, NO 3 , NH 4 , and alkalinity from DLEM as riverine inputs to drive the ocean model in Exp1.Also in Exp3 the air pCO 2 was set to the 1904-1910 condition derived by formula (1).The purpose of Exp2 was to qualitatively examine the role of biological processes in regulating regional pCO 2 variability, whereas Exp3 examined the coastal carbon cycle's response to alternations in river inputs as a result of land-use change within the Mississippi watershed (the first 10 years of the 20th century vs. those of the 21st century).Although we applied riverine and air pCO 2 estimated for the period of 1904-1910, the purpose of Exp3 was not to reproduce the pCO 2 for that period as changes of other variables over the past 100 years were not considered (e.g.air temperature, ocean and food web conditions).

Validation of the control run
We utilized the ship-based sea surface pCO 2 database compiled by the Lamont-Doherty Earth Observatory (LDEO Version 2014, > 180 000 data points in the Gulf over 2005-2010;Takahashi et al., 2015) and Huang et al. (2015) for model validation (see locations of ship measurements in Fig. 5).The ship measurements by Huang et al. (2015) were taken in October 2005; April, June, and August 2006; May and August 2007; January, April, July, and November 2009; and March 2010 and contain > 78 000 data points.To alleviate the spatial and temporal mismatches associated with these in situ measurements, we computed their temporal and spatial mean using a 10-day temporal binning for temporal processing and then compared them with model-simulated pCO 2 time series (Fig. 6).To facilitate our analysis, the GoM was divided into five subregions: (1) Mexico Shelf (MX Shelf), (2) Western Gulf of Mexico Shelf (WGoM Shelf), (3) Northern Gulf of Mexico Shelf (NGoM Shelf), (4) West Florida Shelf (WF Shelf), and ( 5) the open ocean, which is > 200 m water depth (regional definitions followed Benway and Coble, 2014, maps of subregions see Fig. 1).The data points falling in each of the subregions was first grouped by a 10-day temporal binning and then spatially averaged to get a mean value for each subregion.On the NGoM Shelf, the control simulation was able to capture the measured pCO 2 in 21 out of the 26 data groups (the mean value of in situ measurements fell within

Results
In this section, we present model-simulated sea surface pCO 2 and air-sea CO 2 flux in the five subregions.Because few data existed and large pCO 2 gradients were found in both in situ measurements and model simulation in shallow waters, areas that are shallower than 10 m were excluded from our analysis.

Temporal variability of sea surface pCO 2
Spatially averaged model-simulated pCO 2 on the NGoM Shelf exhibited clear seasonality, with highest values (∼ 500 ppm) around August and lowest values (∼ 300 ppm) around February (Fig. 6a).Notably, spatially averaged pCO 2 on the NGoM Shelf was not coincident with high river carbon and nutrient inputs (Fig. 3).Peaks in pCO 2 generally occurred 2 to 3 months later than the annual maximum in river input.The maximum riverine input during 2005-2010 was observed in June 2008 when a major flood occurred (Fig. 4a), yet no significant elevation of pCO 2 was seen in the model simulation.Gulf-wide spatially averaged pCO 2 (Fig. 4b) had a temporal pattern similar to that on the NGoM Shelf, with high pCO 2 values (∼ 425 ppm) in August and low values (∼ 350 ppm) in February.Averaged pCO 2 on the NGoM Shelf was generally 50 ppm higher than that in the entire Gulf.

Model simulations of air-sea CO 2 flux
The simulated carbon flux was calculated from a multiyear model mean (2005)(2006)(2007)(2008)(2009)(2010).We found that the GoM overall was a CO 2 sink with a mean flux rate of 0.71 ± 0.54 mol C m −2 yr −1 (∼ 1.11 ± 0.84 × 10 12 mol C yr −1 , Table 1 and Fig. 7).Examining region by region, we found that the open ocean, occupying ∼ 65 % of the GoM by area, acted as a CO 2 sink (1.04 ± 0.46 mol m −2 yr −1 of C) during most of the year except in summer.The greatest carbon uptake occurred in winter (2.44 ± 0.49 mol C m −2 yr −1 ).It is evident that waters around the Loop Current act as a sink throughout the year, whereas the western part of the open ocean waters shifted from acting as a CO 2 source in summer and fall to a sink in winter and spring.
Compared with the open ocean, air-sea flux on the continental shelf was more location-dependent and varied from season to season.Among the four shelf subregions, the MX Shelf has the largest area.It acted as a strong CO 2 sink in winter and spring (0.49 ± 0.28 and 0.97 ± 0.28 mol C m −2 yr −1 ) and then a carbon source in summer and fall (−0.96± 0.38 and −0.76 ± 0.45 mol C m −2 yr −1 ).Waters along the eastern side of the MX Shelf were a sink during most of the year, while to the west the shelf was a source in summer and fall.On an annual scale, this region was a sink with an air-sea flux of 0.19 ± 0.35 mol C m −2 yr −1 .To the north, the WGoM Shelf has the smallest area among the four shelf subregions.It acted as a CO 2 source during spring, summer, and fall (−0.24± 0.59, −1.69 ± 0.43, and −1.06 ± 0.34 mol C m −2 yr −1 ) and a strong CO 2 sink during winter (1.62 ± 0.32 mol C m −2 yr −1 ).On an annual scale the WGoM region was a CO 2 source with a degassing rate of 0.34 ± 0.42 mol C m −2 yr −1 .
The NGoM Shelf shifted from acting as a CO 2 source in summer and fall (−1.42 ± 0.74 and −0.79 ± 0.63 mol C m −2 yr −1 ) to a sink in winter and spring (1.01 ± 0.89 and 2.49 ± 0.70 mol C m −2 yr −1 ).The most prominent feature here was the continuous, strong degassing in the coastal waters around the Mississippi-Atchafalaya River mouths.However, as the water becomes deeper, the NGoM Shelf water shifted from acting as a sink during winter and spring to a source during summer and fall.Despite the extensive degassing in the coastal water, the NGoM Shelf overall was a CO 2 sink on a yearly basis (0.32 ± 0.74 mol C m −2 yr −1 ).Similarly, the WF Shelf also shifted from acting as a CO 2 source in summer and fall (−1.26± 0.53 and −1.73 ± 0.67 mol C m −2 yr −1 ) to a sink in winter and spring (1.19 ± 0.38 and 0.28 ± 0.33 mol C m −2 yr −1 ).The degassing in the inner shelf was strong enough to make the WF Shelf a CO 2 source on a yearly basis (−0.38 ± 0.48 mol C m −2 yr −1 ).
Despite the salient spatial and temporal variability, the GoM was an overall CO 2 sink, mainly because of the strong uptake in the open ocean.For validation purposes, we compared (in Table 1) model-simulated air-sea flux against an estimation based on observations, which utilized all available measurements collected within the GoM from 2005 to 2010 (Robbins et al., 2014).Our control-run estimations generally agree with in situ measurements in all five subregions in terms of the ocean's role as a CO 2 source or sink.There is some discrepancy in the magnitude of the estimated flux, specifically in the open ocean subregion.We note that Robbins et al. (2014) used monthly mean pCO 2 and wind fields in their calculation as opposed to the 10-day interval we used here.Therefore, to facili-

Net community production (NCP)
As NCP plays an important role in regulating water CO 2 con-  mer, with the highest value (2.62 mmol N m −3 ) simulated in summer 2008 when there was a major flooding event.Compared with the NGoM condition (0.53 mmol N m −3 ), mean surface NCP in the open ocean was relatively small, with a multiyear mean value of 0.11 mmol N m −3 .In addition, the Gulf-wide mean surface NCP exhibited peaks in late winter and early spring, mainly incurred by the strong upwelling along the western side of the Yucatán Strait (Fig. 8a and d).
Compared with the surface NCP, the magnitude of bottom NCP was found to be small and is thus not shown.

Model sensitivity experiments: no-biology simulation (Exp2)
To qualitatively test the role of biological processes in regional CO 2 variability, a no-biology run was conducted, in which all biology sources and sinks of DIC and alkalinity were disabled similar to the experiment described in Fennel and Wilkin (2009).The experiment produced higher surface pCO 2 than the control run.pCO 2 is strongly elevated around the Mississippi River delta on the NGoM Shelf dur-ing spring and summer.For the open ocean, the pCO 2 increase was mainly confined within the Loop Current and was strongly impacted by Caribbean waters flowing in through the Yucatán Channel (Fig. 10).To assess the influence of NCP on CO 2 variation, we plotted the pCO 2 difference between the control run (Exp1) and no-biology run (Exp2) against the surface NCP from the control run in Fig. 11.In the NGoM, the pCO 2 difference between the control run and no-biology run was strongly correlated with NCP (r = 0.80), indicating a regional biological carbon removal.For the open ocean, the pCO 2 difference shows no correlation with NCP, and we speculate that the biological carbon removal in this region was incurred by not only local NCP but also remote processes.As shown in Fig. 9, the poor correlation between pCO 2 and local NCP could be the result of the high pCO 2 water from the Caribbean, which will be discussed in Sect.5.2.The multiyear mean sea surface pCO 2 from the nobiology run was elevated by 88.0 ppm (from 393.1 to 466.5 ppm) for the NGoM Shelf and 56.0 ppm (from 375.1 to 463.1 ppm) for the entire Gulf (Fig. 6, spatially averaged over the subregions).This pCO 2 increase was not temporally uniform.On the NGoM Shelf, pCO 2 increases in the nobiology run were clearly higher during spring-summer (with increases of 84.1 and 95.6 ppm) than during fall-winter (with increases of 57.3 and 56.0 ppm).On the Gulf-wide scale, the pCO 2 increase was stronger during summer (97.1 ppm) than the other seasons (86.5, 87.6, and 80.9 ppm for spring, fall, and winter).For air-sea flux, the elevated surface pCO 2 turns all five subregions into a carbon source throughout the year, resulting in a net outflux rate of 2.10 mol C m −2 yr −1 (Table 1).

Model sensitivity experiments: historical river forcing (Exp3)
The purpose of Exp3 was to examine coastal carbon dynamics' response to different river conditions.Figure 4 shows that river discharge and DIC inputs during years 1904-1910 as simulated by the DLEM are comparable with those at present (2004)(2005)(2006)(2007)(2008)(2009)(2010).The multiyear mean value of freshwater discharge is 25 700 m 3 s −1 for 1904-1910 and 23 900 m 3 s −1 for 2004-2010.The Mississippi-Atchafalaya delivered 1.51 × 10 12 mol C yr −1 during 1904-1910 and 1.70 × 10 12 mol C yr −1 during 2004-2010, which is comparable to the increase over the preceding century reported by Raymond et al. (2008), i.e., a 0.24 × 10 12 mol C yr −1 increase in an average discharge year.However, NO 3 inputs during 1904-1910 were < 30 % of current inputs (18.12 vs. 63.18× 10 9 mol N yr −1 ).Limited N input led to a smaller primary production not only on the NGoM Shelf but also in the adjacent waters on the WGoM and WF shelves.Due to the smaller primary production the coastal ocean was a weaker CO 2 sink during spring and summer (Fig. 12) and the NGoM Shelf a year-long carbon source with a net outflux rate of 0.61 mol C m −2 yr −1 (Table 1).A close examination of the spring and summer conditions on the NGoM Shelf shows that differences in primary production between Exp1 and Exp3 occur mainly along the Texas and Louisiana coasts.Primary production was significantly elevated in the control run because of enhanced NO 3 inputs (Fig. 12a and  c).Elevated primary production brought down the sea surface pCO 2 .During spring, enhanced primary production and www.biogeosciences.net/13/4359/2016/Biogeosciences, 13, 4359-4377, 2016 decreased CO 2 was simulated along the Louisiana and Texas coast (Fig. 12b), while during summer, when coastal circulation was influenced by westerly winds, the decreased pCO 2 was more confined within waters along the Louisiana coast.

Discussion
Prior to this investigation, the carbon dynamics in the GoM have been poorly characterized and had a high degree of uncertainty.This study provides one of the first attempts to simulate GoM-wide carbon fluxes and exchanges using a coupled physical-biogeochemical model.We next discuss the factors controlling sea surface pCO 2 variability on the riverinfluenced NGoM Shelf and the Loop Current-influenced open ocean.The relationship between pCO 2 and other hydrographic variables as well as model uncertainty are also considered.

NGoM Shelf
The Mississippi-Atchafalaya River and associated plume play the most important role in determining the pCO 2 distribution on the NGoM Shelf.The large input of fluvial DIC and alkalinity introduces carbonate saturation in the coastal waters and, conversely, nutrients from the river enhance local primary production, which results in DIC removal and thus reduces sea surface pCO 2 (e.g.Lohrenz et al., 2010;Guo et al., 2012;Huang et al., 2013Huang et al., , 2015)).Such biological removal of CO 2 was also confirmed by the elevated pCO 2 values in the no-biology run in this study.Although the river plume's influence on CO 2 flux has been addressed by prior observational studies, large uncertainties were also found regarding whether the NGoM Shelf is a CO 2 sink or source over a longer time period.For instance, Huang et al. (2013) found a large difference between the pCO 2 distributions in April 2009 and in March 2010.Such a difference was attributed to the variations in river plume extension influenced by local wind conditions and river discharge.In a later communication, based on ship measurements from 11 cruises, Huang et al. (2015) concluded that the NGoM Shelf acted as a net CO 2 sink, but with a large uncertainty (influx rate: 0.96 ± 3.7 mol m −2 yr −1 ).
Model results in this study revealed significant spatial and temporal gradients in sea surface pCO 2 as well.The multiyear mean (2005-2010) pCO 2 distribution was characterized by high values in the coastal waters (Fig. 13a), accompanied by low-salinity (Fig. 13c), high DIN, and high DIC (Fig. 13d  and e).The pCO 2 value was significantly lower as water became deeper, where the ocean acted as a CO 2 sink during most times of the year (Figs.7a through d).The surface pCO 2 distribution on the NGoM Shelf was highly correlated with surface salinity (r value: −0.81) and DIN concentra-  tion (r value: 0.80) throughout the year, while its correlations with surface temperature and DIC concentration were significant only for part of the year (for detailed season-byseason correlations see Table 2).Although our model suggests that the shelf-wide pCO 2 distribution was positively correlated with DIN concentration, this is not contrary to findings of the above-mentioned observational studies; that is, the high DIN stimulates primary production should be negatively correlated with sea surface pCO 2 .Instead, the high DIN concentration, together with the low salinity, was a signal of rich DIC from the riverine inputs and, potentially, the light-limited conditions due to the high suspended sediment and dissolved organic matter concentrations within the river plume.In other words, CO 2 outgassing from oversaturated plume water overwhelmed the CO 2 influx induced by "biological pump" in the areas near the river mouths.
To further link pCO 2 dynamics with biological processes stimulated by river inputs, we plotted the pCO 2 and DIC averaged over spring and summer seasons (high flow from the Mississippi) against surface salinity of the control run and no-biology run in Fig. 14.Seawater pCO 2 decreased almost linearly as salinity increased in the no-biology run (Fig. 14b).During spring and summer when river discharge and DIC inputs were high, the high-pCO 2 and low-salinity waters around the Mississippi River delta (86-88 • W, reddish points) can be easily differentiated from the highsalinity and low-pCO 2 waters on the Texas Shelf (92-95 • W, bluish points).The DIC-salinity relationship for waters around the Mississippi River delta (reddish points in Fig. 14d) fell below the conservative mixing relationship for the river end member calculated using in situ data collected in the spring and summer of 2008 by Cai et al. (2011a).For locations to the west, the DIC-salinity relationship reflected a mixture of waters from the Texas Shelf (bluish points) and those from the Atchafalaya River (yellowish-greenish points) likely with differing end members.
When biological processes were included, the shelf water exhibited large spatial and seasonal variability (left panels).A pCO 2 minimum was simulated in mid-salinity waters (30-33 psu) during spring and summer, which is consistent with the curve derived by Huang et al. (2015) using ship measurements.Compared with the no-biology run, pCO 2 was reduced significantly and exhibited a wider range in the control run.The biological removal of sea surface CO 2 was most salient in waters around the Mississippi River delta.The difference in pCO 2 between waters around the delta and the Texas Shelf became more salient.The DIC-salinity relationship for locations around the Mississippi River delta (reddish points in Fig. 14c) indicated a significant carbon removal along the salinity gradient.For waters on the Texas Shelf, the DIC-salinity relationship was confined to higher salinities and slightly increased compared with the no-biology run (bluish points in Fig. 14c).The DIC increase on the Texas Shelf in the control run could be linked with the benthic respiration in this region proposed by Hetland and Di-Marco (2007).

Open ocean
In the open ocean, the distribution of surface pCO 2 was strongly related to the surface DIC (r value: 0.93) and alkalinity throughout the year (r value: −0.85, for detailed season-by-season correlations see Table 2).An influence of DIN and primary production was evident in fall and winter months when wind-induced upwelling was strong (Xue et al., 2013).The dependence of pCO 2 on DIC and alka-linity makes the Loop Current an important factor controlling the regional air-sea CO 2 flux.In addition to a relatively high temperature, the Loop Current water is also characterized by low DIC and high alkalinity (Wang et al., 2013, and references therein).The multiyear mean sea surface temperature in Fig. 13b shows persistent warm water mass in the form of the Loop Current, which carries the carbonate characteristics of the Caribbean water (i.e.low DIC and high alkalinity, Fig. 13e and f).Surface pCO 2 in this warm water mass was significantly lower than surrounding shelf waters (Fig. 13a), making the Loop Current a strong CO 2 sink throughout the year (Fig. 7a-d).Any changes in the Caribbean water's carbonate characteristics will affect the carbon budget in the GoM as well as waters further downstream in the Gulf Stream.This is also illustrated by the high pCO 2 difference between the control run and no-biology run in Fig. 10 as well as the poor correlation between the pCO 2 drop (difference between control and no-biology runs) and NCP in the open ocean (Fig. 11b).

Carbon budget estimation and model uncertainty
Based on our model simulations, we conclude that the GoM is an overall CO 2 sink, taking up 1.11 ± 0.84 × 10 12 mol C yr −1 from the air.This estimation is comparable to those based on in situ observations, e.g.1.48 × 10 12 mol C yr −1 (Coble et al., 2010) and 0.30 × 10 12 mol C yr −1 (Robbins et al., 2014).These recent estimates are in stark contrast to the earlier SOCCR report (Takahashi et al., 2007), which found the GoM to be a CO 2 source (1.58 × 10 12 mol C yr −1 , the GoM and Caribbean Sea combined).In addition, we estimated that the GoM received ∼ 2.18 × 10 12 mol C yr −1 from rivers, the majority of which was from the Mississippi-Atchafalaya River (∼ 1.80 × 10 12 mol C yr −1 ).These two DIC sources (air, ∼ 1.11 × 10 12 mol C yr −1 , plus river, ∼ 2.18 × 10 12 mol C yr −1 ) is comparable to the DIC transported out of the GoM by the Loop Current (∼ 3.30 × 10 12 mol C yr −1 ; Wang et al., 2013).Such a balance cannot be achieved using the CO 2 flux estimated by Robbins et al. (2014).Nevertheless, here our intent is not to close the carbon budget, considering the large uncertainties involved and discussed below.Indeed, the ultimate CO 2 source and/or sink term would be dependent on the relative contribution of both DIC and nutrients to the upper layer of the ocean as well as the biogeochemical alteration therein (Dai et al., 2013).We notice that, during summer months, our model simulated a higher surface pCO 2 than ship measurements on the NGoM Shelf (Fig. 6a).As discussed in Sect.5.1, a large part of the strong CO 2 degassing was simulated on the Texas Shelf.However, a close examination of the distribution of available ship measurements indicates that data points on the Texas Shelf are fairly sparse and sporadic (Fig. 5), which may partially explain the mismatch between model and ship measurements in Fig. 6a.For instance, in the summer of 2010 when more ship measurements were available on the NGoM Shelf, both model and observation indicated a high pCO 2 in the summer.In addition, pCO 2 in the Mississippi plume was very sensitive to river DIC inputs.Our specification of riverine DIC (e.g.alkalinity plus 50) was based on limited measurements and may not reflect the true seasonal and interannual variability of alkalinity-DIC relationship.The current model resolution (∼ 5 km) may not be high enough to reproduce small-scale circulation patterns associated with the Mississippi River plume.The complexity of the food web and uncertainty in model parameterization (e.g.rudimentarily represented denitrification, remineralization, particular organic matters, the lack of phosphate and silicate components) warrant further investigation.

Summary
A coupled physical-biogeochemical model was used to hindcast surface pCO 2 in the GoM from January 2004 to December 2010.Favorable comparisons were found when validating model solutions against ship measurements on the Gulfwide scale, indicating that this coupled model can reproduce observed pCO 2 variability in the GoM.Time series of spatially averaged pCO 2 for both shelf and open ocean waters exhibit significant seasonal variability, with high values in August and low values in February.Model-simulated pCO 2 values were elevated by 56 and 88 ppm for the entire Gulf and the NGoM Shelf, respectively, when the biological sources and sinks of carbon were disabled (i.e., the no-biology run).Without biological processes, the GoM shifts to a strong carbon source with a outflux rate of 2.10 mol C m −2 yr −1 .Another sensitivity test examining river conditions from the 1904-1910 period (reduced NO 3 and comparable DIC) supported the view that the impact of river inputs were mainly limited to the NGoM Shelf, which under the conditions of the simulation acted as a CO 2 source with an outflux rate of 0.61 mol C m −2 yr −1 .
The Mississippi-Atchafalaya River plume is the dominant factor controlling the pCO 2 distribution on the NGoM Shelf.Although the NGoM Shelf is overall a CO 2 sink, highsurface pCO 2 was simulated in relatively shallow waters, induced by both oversaturated plume water.pCO 2 in the open ocean is controlled largely by the low-DIC, high-alkalinity Loop Current water from the Caribbean Sea.
Our model simulations characterize the GoM as an overall CO 2 sink, taking up ∼ 1.11 ± 0.84 × 10 12 mol C yr −1 from the air.Together with the enormous riverine input (∼ 2.18 × 10 12 mol C yr −1 ), this inorganic carbon influx was comparable with the DIC export through the Loop Current estimated by an earlier study.More accurate model predictions of water column DIC concentration will require more in situ data for improved specification of riverine DIC inputs, model DIC initial conditions, and further process studies to refine model parameterizations so as to better account for complex carbon dynamics in the coastal ocean.

Data availability
The operational mode of the SABGOM model is located at http://omgsrv1.meas.ncsu.edu:

Figure 1 .
Figure 1.Domain of the South Atlantic Bight and Gulf of Mexico (SABGOM) ROMS model with water depth in color (unit: m).Also shown are the five subregions used in this study, which are Mexico Shelf (MX), Western Gulf of Mexico Shelf (WGoM), Northern Gulf of Mexico Shelf (NGoM), West Florida Shelf (WF), and open ocean.Also shown is a schematic for the Loop Current.

Figure 3 .
Figure 3. Satellite-observed monthly pCO 2 (AIRS) averaged over the Gulf of Mexico (red stars) and the pCO 2 air used in model airsea CO 2 flux calculation (blue line), which is generated using the curve-fitting software CCGCRV.

Figure 4 .
Figure 4. Comparisons between the 2005-2010 riverine DIC and NO 3 conditions observed by USGS (red line) and the 1904-1910 river condition simulated by the Dynamic Land Ecosystem Model (black line; Tian et al., 2015).

Figure 5 .
Figure 5. Locations of in situ measurements from the LDEO database (blue) and Huang et al. (2015, grey) in the period of 2005-2010.

Figure 6 .
Figure 6.Time series of spatially averaged pCO 2 (control run in blue and no-biology run in red) (a) on the Northern Gulf of Mexico Shelf and (b) in the entire Gulf of Mexico, overlaid with in situ observations (in black) from Huang et al. (2015) and Takahashi et al. (2015).

Table 1 .
Comparison between observed and modeled air-sea CO 2 flux.Observations are taken from Robbins et al. (2014), whereas the model results are 7-year (2005-2010) model means * .Subregions Mexico Shelf Western Gulf Northern Gulf West Florida Shelf Open Ocean Gulf-wide * * tate the comparison of results, we recalculated the flux using a monthly mean pCO 2 and wind fields and obtained a flux estimate of 0.31 ± 0.35 mol C m −2 yr −1 for the open ocean subregion and 0.12 ± 0.23 mol C m −2 yr −1 for the entire GoM.These values are comparable to those in Robbins et al. (2014, 0.48 ± 0.07 mol C m −2 yr −1 for the open ocean and 0.19 ± 0.08 mol C m −2 yr −1 for the entire GoM).

Figure 7 .
Figure 7. Six-year (2005-2010) model (control run) mean air-sea CO 2 flux in the Gulf of Mexico during (a) spring, (b) summer, (c) fall, and (d) winter.Blue color indicates where the ocean is a sink for CO 2 ; red color indicates where the ocean is a source.

Figure 9 .
Figure 9.Time series of spatially averaged net community production (a) on the Northern Gulf of Mexico Shelf and (b) in the entire Gulf of Mexico (unit: mmol N m −3 ).

Figure 11 .
Figure 11.Scatter plots of the multiyear mean pCO 2 drop (no-biology run minus control run) and surface NCP in NGoM (left) and open ocean (right).

Figure 14 .
Figure 14.Six-year (2005-2010) spring-summer mean condition of model-simulated sea surface pCO 2 and DIC against salinity for the control run (a, c) and no-biology run (b, d) on the NGoM Shelf; also shown are longitude with colors (note that the Mississippi River delta is located around 89 • W and Atchafalaya River delta is located around 91 • W).Also shown in (c) and (d) are conservative mixing relationships for river end members from Cai et al. (2011a).
of the model mean).Specifically, agreement between model and observations was better during spring, fall, and winter, than during summer.The model overestimated pCO 2 in June 2006, August 2007, and July 2009.These discrepancies will be discussed in later sections.On the Gulf-wide scale, the control run reproduced the observed seasonality.Decent model-data agreements were found in 24 out of the 26 data groups.These subregional and Gulf-wide comparisons indicate that the coupled physicalbiogeochemical model is generally capable of resolving temporal and spatial variations in observed pCO 2 , allowing us to use this 7-year hindcast to further characterize the air-sea CO 2 flux.