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₂ (pCO₂), aragonite saturation state (Ω_Ar), and sea–air CO₂ flux. Model results indicate that seasonal changes in surface pCO₂ are strongly controlled by temperature across 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₂ fluxes are consistent with previous observation-based estimates that show CO₂ uptake during winter–spring, and CO₂ outgassing during summer–fall. Annually, our model indicates a basin-wide mean CO₂ uptake of 0.35 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, and a northern GoM shelf (< 200 m) uptake of 0.93 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. The observation and model-derived patterns of surface pCO₂ and CO₂ fluxes show good correspondence; thus this study contributes to improved constraints of the carbon budget in the region.

1 Introduction

The global ocean is absorbing approximately one-third of the anthropogenic CO₂ 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 state of calcium carbonate (e.g., Orr et al., 2005). This process, commonly known as ocean acidification, has deleterious impacts on calcifying organisms, such as corals, coralline algae, shellfish, and shell-forming plankton (Doney, 2012). Ocean acidification is disturbing marine ecosystems worldwide (e.g., Mostofa et al., 2016), demanding urgent societal responses to address coastal ecosystem impacts. Therefore, a better understanding of the past and current carbon system variability at global and regional scales is crucial to better monitor and predict ocean and ecosystem responses to enhanced CO₂ levels.

Significant progress has been made in the understanding of ocean carbon dynamics in coastal waters of the United States during the last 15 years or so. However, many aspects remain poorly understood and described (e.g., Chavez et al. 2007; Wanninkhof et al., 2015; Fennel et al., 2019). Uncertainties in carbon system patterns are particularly large in the Gulf of Mexico (GoM), a low-latitude semi-enclosed basin surrounded by the coasts of the southern United States and eastern Mexico (Fig. 1). The GoM encompasses diverse biogeochemical regimes, from the warm and oligotrophic open GoM, strongly influenced by the Loop Current and mesoscale eddies, to wide and productive continental shelves, influenced by river runoff- and wind-driven coastal currents (e.g., Dagg and Breed, 2003; Zavala-Hidalgo et al., 2006; Wang et al., 2013; Muller-Karger et al., 2015; Anglès et al., 2019). Therefore, multiple dynamics modulate the GoM carbon chemistry, which makes reducing uncertainties in these biogeochemical patterns a challenging task.

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

Most observational studies on carbon dynamics in the GoM have been conducted on the Louisiana–Texas shelf (e.g., Cai, 2003; Lohrenz et al., 2010, 2018; Guo et al., 2012; Cai et al., 2013; Huang et al., 2012; 2015; Hu et al., 2018). In this region, the Mississippi–Atchafalaya river system (MARS) has a strong influence, delivering a significant amount of freshwater, carbon, and nutrients, the latter fueling high biological production (Green et al., 2008; Lehrter et al., 2013). Enhanced primary production during spring and summer periods increases carbon uptake near the MARS delta, which results in decreased surface partial pressure of CO₂ (pCO₂) and increased ocean uptake of CO₂ (Lohrenz et al., 2010, 2018; Guo et al., 2012; Huang et al., 2015; Hu et al., 2018). Subsequent sinking and remineralization of large amounts of organic carbon over the Louisiana–Texas shelf, concurrent with strong water column stratification, results in bottom acidification during the summer (Cai et al, 2011). The variability in carbon chemistry for other GoM areas has been less examined, but an increasing number of observations from dedicated research programs (e.g., Gulf of Mexico Ecosystem and Carbon Cycle, or GOMECC) and ship of opportunity (SOOP) programs are contributing to a reduction in the spatial and temporal data gaps. Robbins et al. (2014) derived estimates of sea–air CO₂ fluxes over the entire GoM, concluding that the GoM basin is a CO₂ sink. Recently, Robbins et al. (2018) described pCO₂ patterns on the west Florida shelf, indicating that this region is mainly a CO₂ source with significant spatial and seasonal variability.

Nevertheless, data gaps and observational constraints still limit our ability to infer carbon patterns in the ocean. Thus, regional ocean biogeochemical models that simulate carbon dynamics at multiple timescales are valuable tools to better understand the carbon system variability and its underlying drivers. In the GoM, several three-dimensional modeling studies addressing carbon cycle aspects have been conducted. Xue et al. (2016) used the Fennel biogeochemical model (Fennel et al. 2008; Fennel and Wilkin, 2009) to examine pCO₂ and sea–air CO₂ fluxes during 2005–2010. They reproduced observed spatiotemporal patterns across the GoM to some degree; however, some discrepancies between their model results and in situ observations are noted. For example, their model did not reproduce the decrease in surface pCO₂ linked to high primary production over the MARS mixing zone (Huang et al., 2015), and spatially averaged values of model pCO₂ were largely overestimated in the northern GoM during summer (by more than 100 µatm in several cases). In addition, the modeled sea–air CO₂ flux in the northern GoM (−0.32 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{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 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{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₂, sea–air CO₂ fluxes, pH, alkalinity, and dissolved inorganic carbon (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₂, as well as sea–air CO₂ 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₂ fluxes in coastal and ocean domains; and (5) discuss the main model results in the context of previous observational and modeling studies.

2 Model and data

2.1 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, produced by zooplankton excretion and organic matter remineralization, and affected by sea–air CO₂ fluxes. Changes in model TA are estimated using an explicit conservative expression for alkalinity (Wolf-Gladrow et al., 2007). Model CO₂ fluxes are derived using the Wanninkhof (2014) bulk flux equation. Details of the carbon module can be found in Sect. S1 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 a biogeochemical model (Dunne et al., 2013). 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 access:23 September 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 7 373 420 and 7 381 600. 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 (Sect. 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₂ was prescribed as a continuous nonlinear function, derived from the Mauna Loa monthly CO₂ time series (https://www.esrl.noaa.gov/gmd/ccgg/trends/, last access: 16 August 2018) using similar curve-fitting method that Thoning et al. (1989; Sect. S3).

Figure 1 Model snapshot of surface dissolved inorganic carbon (mmol m⁻³) during 1 May 2009. Regions used to describe model results are the western GoM shelf, the northern GoM shelf, the west Florida shelf, the Yucatan shelf, and 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. The magenta polygon demarks the region near the Mississippi Delta used to derive patterns in Fig. 7.

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₂ and Ω_Ar. The calculations were performed using the MatLab version of the CO2SYS program for CO₂ 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₄ dissociation constant of Dickson (1990b).

Figure 2 Time series of mole fraction of CO₂ (xCO₂), SST, and surface salinity derived from a surface mooring (Coastal Mississippi Buoy) and model outputs at 30^∘ N and 88.6^∘ W (location depicted as red star in Fig. 1). Simulated and observed monthly averages are shown as blue and red lines, respectively. Buoy data (6 h interval) are depicted in magenta.

For the present study, we focused on describing seasonal patterns in surface Ω_Ar, surface pCO₂, and sea–air CO₂ flux during 2005–2014 (i.e., the last 10 years of the model run). Ω_Ar represents the degree of saturation of calcium carbonate (CaCO₃) 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

$$\begin{array}{}\text{(1)}& {\mathrm{\Omega }}_{\text{Ar}}=\left[{\mathrm{Ca}}^{\mathrm{2}+}\right]\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]{\left(K_{\text{Ar}}^{\prime }\right)}^{-\mathrm{1}},\end{array}$$

where [Ca²⁺] is total calcium concentration, which is a function of salinity, [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] is total carbonate ion concentration, which is derived from the simulated DIC and TA, and $K_{\text{Ar}}^{\prime }$ is the apparent solubility product of the CaCO₃ 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 in Ω_Ar mainly depend on [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$], and are positively related to changes in the TA : DIC ratio (Wang et al., 2013).

2.2 Data

Surface measurements of mole fraction of CO₂ (xCO₂), 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) were retrieved from the NOAA National Center for Environmental Information (https://www.nodc.noaa.gov, last access: 4 March 2019). Vertical profiles for DIC, TA, temperature, and salinity off Tampa (Florida) and Louisiana were derived from measurements collected during the GOMECC-1 cruise; Wang et al., 2013), retrieved from NOAA-AOML (http://www.aoml.noaa.gov/ocd/gcc/GOMECC1, last access: 4 March 2019). Surface pCO₂ data were obtained from underway measurements collected onboard research cruises and multiple ships of opportunity, and compiled by Barbero et al. (in preparation). The pCO₂_GoM_2018 dataset, which contains more than 457 000 measurements in the GoM during 2005–2014 (Fig. S5), is available as a data package from NCEL.

3 Model–data comparison

We used data from the Coastal Mississippi Buoy to evaluate the model's ability to reproduce coastal patterns in xCO₂, 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 cross-shore 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₂, making the exact reproduction of observed buoy patterns challenging. The best match between simulated and observed xCO₂ was during 2011–2012, where xCO₂ ranged from about 230 ppm in spring to more than 400 ppm in fall.

The pCO₂GoM_2018 dataset was used to compare climatological seasonal patterns in pCO₂ (Fig. 3). Overall, simulated and observed pCO₂ 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₂GoM_2018 dataset (Fig. S6 in the Supplement). Indeed, January observations came from only one cruise, which largely increases observational uncertainty. A spatial visualization of the pCO₂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₂ minimum near the MARS region, and the large seasonal amplitude in the western Florida shelf.

Figure 3 Mean monthly patterns for the observed (red lines) and simulated (blue lines) surface pCO₂ over the (a) open GoM and (b) northern GoM regions (shown in Fig. 1). Light pink and cyan shading depict the observed and modeled interquartile interval, respectively. Gray shading depict the model's 5–95 % percentile interval. Observations are from ships of opportunity and research cruises conducted during 2005–2014 (ship tracks are shown in Fig. S4.1).

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 (Mississippi line), with the bias ranging from around 5 to 90 µmol kg⁻¹ for DIC and 5 to 40 µmol kg⁻¹ 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⁻¹ 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₂ and Ω_Ar for the 2005–2014 study period.

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 Fig. 1 as black stars). The range of the model's variables for June–August during 2000–2014 is also shown as cyan shading.

4 Surface pCO₂ and Ω_Ar seasonality

Model-derived patterns for surface pCO₂ showed significant seasonal variability across the GoM (Fig. 5). Minimum and maximum pCO₂ 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₂ 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₂ values during spring and summer (Fig. 5b and c). The continental shelf with the highest seasonally averaged pCO₂ was the west Florida shelf, where pCO₂ reached values greater than 450 µatm during the summer. Seasonality in modeled pCO₂ was strongly modulated by sea surface temperature (SST), such that the annual amplitude for these two variables displayed very consistent spatial patterns (Figs. 6a, b and S7). The largest annual signal for pCO₂ 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₂ and SST were strongly correlated in all regions except near the MARS delta (Fig. 6c).

Figure 5 Mean model surface pCO₂ (uatm) in winter (DJF), spring (MAM), summer (JJA), and fall (SON) during 2005–2014. The black contour depicts the 200 m isobath.

Figure 6 (a, b) Seasonal amplitude patterns for model surface pCO₂ 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₂ and SST. Black contour depicts the 200 m isobath.

The low pCO₂-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₂ along the salinity gradient linked to the Mississippi River plume (Fig. 7). The simulated surface pCO₂ patterns during spring and summer displayed a marked increase from middle to low salinities (Fig. 7a and d), which was also associated with an increase in DIC (Fig. 7b and e). The minimum pCO₂ 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₂ 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 middle salinities mainly driven by NCP, indicative of a biologically induced drawdown of pCO₂. During fall (Fig. 7g–i), as well as winter (not shown), NCP was much smaller than during spring–summer, and DIC was mainly controlled by sea–air exchange and advection plus mixing processes. As a consequence, model surface pCO₂ did not show a middle salinity minimum linked to phytoplankton uptake.

Figure 7 Mean patterns of simulated surface variables as a function of salinity near the Mississippi River (magenta polygon in Fig. 1) during spring (a–c), summer (d–f), and fall (g–i): (a, d, g) pCO₂ and pCO₂ normalized to 25 ^∘C; (b, e, h) dissolved inorganic carbon (DIC) and total alkalinity (TA); (c, f, i) budget terms for DIC: advection plus mixing (Adv + Mix), sea–air CO₂ flux (Sea–Air), and net community production (NCP). Thin dashed lines demarcate the interquartile interval (between percentiles 25 % and 75 %). Only results for salinities greater than 17 are shown, since the spatiotemporal resolution from the monthly model outputs did not resolve features at lower salinities well.

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 and b). Although the seasonal patterns for Ω_Ar and pCO₂ 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. 5b and c and 8b and c), when the highest Ω_Ar and lowest pCO₂ values were co-located near the MARS delta, and the lowest Ω_Ar and highest pCO₂ values were in the west Florida shelf and the western part of the northern GoM shelf. 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 (Fig. 9c). This pattern was consistent with enhanced biological uptake of DIC promoted by MARS's nutrient inputs, and the positive salinity impact on aragonite solubility.

Figure 8 Mean model surface aragonite state in winter (DJF), spring (MAM), summer (JJA), and fall (SON) during 2005–2014. The black contour depicts the 200 m isobath.

Figure 9 Correlation between surface aragonite saturation state and surface (a) temperature, (b) TA : DIC ratio, and (c) salinity. The black contour depicts the 200 m isobath.

To better describe the impact of SST in the simulated pCO₂ and Ω_Ar variability, we calculated average monthly climatologies for temperature-normalized pCO₂ 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₂ was evident when we compared the monthly climatologies for pCO₂ and pCO_{2_at25} (Fig. 10a and b). Surface pCO_{2_at25} displayed much weaker annual variation than surface pCO₂, 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₂ 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 and 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 on DIC and TA (Ω_Ar increases with TA and decreases with DIC, while pCO_{2_at25} has the opposite pattern).

Figure 10 Figure 10. Monthly climatology for model (a) pCO₂, (b) pCO₂ at 25 ^∘C, (c) aragonite saturation state (Ω_Ar), (d) Ω_Ar at 25 ^∘C, (e) dissolved inorganic carbon (DIC), and (f) total alkalinity (TA) in northern GoM shelf (nGoM; blue), west Florida shelf (wFL; green), western GoM shelf (wGoM; cyan), Yucatan shelf (black), and open GoM (oGoM, red). Patterns were derived for 2005–2014.

Simulated climatological patterns for DIC and TA (Figs. 10e, f; 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₂ (see Sect. 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 lower long-term average DIC and TA 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 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).

5 Sea–air CO₂ fluxes

Seasonal changes in surface model pCO₂, mainly driven by SST changes (Fig 6c), determined strong seasonal variability in simulated sea–air CO₂ fluxes. As a consequence, the GoM becomes a CO₂ sink in winter–spring and a CO₂ source in summer–fall (Fig. 11a–d). An exception to this pattern occurred close to the MARS delta, which is predominantly a CO₂ sink year-round. In this region, the pCO₂ drop induced by phytoplankton uptake during spring–summer (Fig. 7a and d) determined maximum uptake of atmospheric CO₂ at middle salinities (seen in the sea–air exchange term in Fig. 7c and f). The greatest model CO₂ uptake, above 7 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, occurred over the northern GoM shelf during winter, as this region experiences the lowest surface pCO₂ values induced by the coldest winter conditions in the region (Fig. S7). The greatest model CO₂ 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 and western GoM), and western Yucatan Peninsula during the summer, ranging from ∼ 2 to 3 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (Fig. 11c). Maximum SST values characterized summer conditions in these regions (Fig. S7). The annual mean pattern showed modeled CO₂ uptake ranging from −4 to −1 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in the northern GoM, and from −2 to 0 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ elsewhere (Fig. 11e). In addition, the pattern revealed areas where CO₂ 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).

Figure 11 Model sea–air CO₂ flux ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) patterns during 2005–2014. (a–d) Spatial mean patterns for (a) winter (DJF), (b) spring (MAM), (c) summer (JJA), and (d) fall (SON). (e) Spatial annual mean. (f) Monthly climatology for the northern GoM shelf (nGoM; blue), west Florida shelf (wFL; green), western GoM shelf (wGoM; cyan), Yucatan shelf (black), and open GoM (oGoM, red). Negative (positive) flux implies ocean uptake (degassing). The magenta contours in panels (a–e) depict 0 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, and black contours the 200 m isobath.

The estimated monthly patterns for modeled sea–air CO₂ flux revealed prevailing CO₂ outgassing during May–October in west Florida, western GoM, and Yucatan Peninsula, and June–October in the northern and open GoM (Fig. 11f; Table 1). The timing for the maximum CO₂ outgassing was June–July in the western GoM, August in west Florida and Yucatan, and September in the northern and open GoM. The timing for the maximum CO₂ 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 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively. For the entire GoM basin, the simulated average annual flux and standard deviation was −0.97 and 2.78 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (−0.35 and 1.01 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), respectively. Integrated across the entire model domain, the resulting flux was −7.0 Tg C yr⁻¹.

6 Discussion

6.1 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₂ 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₂ 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 Sect. 3, the model simulated the main surface spatiotemporal 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₂, with relatively minor pCO₂ 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 TA, improving the representation of land–ocean biogeochemical fluxes (e.g., prescribing time evolving TA and DIC for rivers other than the MARS), and increasing the model's horizontal resolution to resolve sub-mesoscale dynamics. Our current model configuration represents an important advance in the model capabilities for the GoM, capturing realistically dominant seasonal patterns.

Table 2 Comparison between annual sea–air CO₂ fluxes ($\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) derived from our model results and previous studies in the Gulf of Mexico. Standard deviation is shown in parentheses. Negative flux implies ocean CO₂ uptake, and positive flux CO₂ outgassing (shown in bold). Shelf regions are depicted in Fig. 1.

1: Regional study; 2: global study; 3: model-based; 4: observational-based; 5: gridded dataset; 6: Margins and Catchments Segmentation (MARCATS) dataset.

Simulated patterns in surface pCO₂ 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₂ 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₂-SST correlation pattern derived from the model is consistent with previous observational studies, which suggested an increased correlation between pCO₂ 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 by Gledhill et al. (2008). Our model results also agree with observations by Guo et al. (2012), Wang et al. (2013), and Wanninkhof et al. (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 on changes in [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{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 the 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₂ fluxes (warm conditions favor surface DIC decrease due to CO₂ 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 signatures 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₂ flux derived from the model output shows that the GoM is a CO₂ sink during winter–spring, and a CO₂ source during summer–fall. However, significant differences in the annual flux magnitude were observed among regions, which could 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₂ and the largest CO₂ 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₂ levels and CO₂ outgassing at low salinities (< 20), linked to the low productivity, high turbidity, and CO₂ oversaturated waters delivered by the Mississippi River; (ii) minimum pCO₂ values and maximum atmospheric CO₂ uptake at middle salinities (20–33), as high phytoplankton production, induced by the water's lower turbidity and decreased nutrient runoff, produces a drop in surface DIC, and (iii) increased pCO₂ levels and sea–air CO₂ 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₂ and sea–air CO₂ flux seasonality. As a result, the annually integrated sea–air CO₂ flux (per m²) 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₂ and CO₂ 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₂ 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₂ 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 through the end of the twenty-first century (Liu et al., 2012, 2015; Alexander et al., 2020; Shin and Alexander, 2020). This is topic deserves examination in future modeling efforts.

6.2 CO₂ flux comparison

Table 2 shows mean CO₂ fluxes derived from our model, previous regional studies for the GoM, and global datasets. The regional-scale studies are Robbins et al. (2014; 2018), Huang et al. (2015), Xue et al. (2016), and Lohrenz et al. (2018). The global-scale studies include Takahashi et al. (2009), Rödenbeck et al. (2013), Landshützer et al. (2016), Laruelle et al. (2014), and Bourgeois et al. (2016). Annual CO₂ fluxes for the GoM basin displayed a significant dispersion, ranging from −0.72 to +0.20 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{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 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, which is comparable with Robbins et al. (2014) and Xue et al. (2016) estimates. In contrast, two of three basin fluxes derived from global gridded datasets, Takahashi et al. (2009) and Landshützer et al. (2016), suggest that the GoM is a weak CO₂ source. This discrepancy between regional and global studies most likely reflects inaccuracy in global datasets, due to the low density of pCO₂ observations in the GoM basin and coarse grid resolutions (5^∘ latitude × 4^∘ longitude in Takahashi et al. 2009 and 1^∘ latitude × 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 Fig. 1. In the open GoM region, our mean flux (−0.33 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{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 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) is 20 % above the value reported by Laruelle et al. (2014). In the northern GoM, our simulated flux (−0.93 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) is remarkably similar to the reported fluxes of Huang et al. (2015) and Lohrenz et al. (2018; −0.95 and −1.1 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{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 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{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 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{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₂ uptake is possible, as the model surface pCO₂ 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₂ measurements, used to generate the observed bulk CO₂ 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 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{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 in 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₂ fluxes can be explained mainly by pCO₂ biases in the model used in Xue et al. (2016). Indeed, their model underestimated surface pCO₂ in the open GoM, and thus obtained a marked pCO₂ 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₂ on the northern GoM and west Florida inner shelves, especially during summer–fall, not reproducing well the marked pCO₂ drop that is observed close to the MARS delta.

7 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₂ than previous modeling efforts in the region. We described seasonal patterns in surface pCO₂ and Ω_Ar. Both variables show maximum values during late summer and minimum during winter and early spring. The seasonal cycle for pCO₂ 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₂, and sea–air CO₂ flux seasonality. Simulated fluxes show CO₂ uptake prevailing during winter–spring, and CO₂ outgassing during summer–fall. The integrated annual flux for the GoM basin is −0.35 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (−4.2 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). The largest model CO₂ uptake is in the northern GoM shelf, linked to the most intense winter cooling, and significant biological uptake during spring–summer. The weakest CO₂ 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 spatiotemporal coverage in the underway pCO₂ measurement dataset over those two regions, especially during wintertime.