Hydrologic controls on CO2 chemistry and flux in subtropical lagoonal estuaries of the northwestern Gulf of Mexico

Estuaries are generally considered a source of CO2 to the atmosphere, although with significant uncertainties in magnitude and controlling factors between and within estuaries. We studied four northwestern Gulf of Mexico estuaries that experience extreme hydrologic conditions between April 2014 and February 2017 to determine the role of dry/wet cycle on estuarine CO2 system. Annual air–water CO2 flux ranged from 2.7 to 35.9 mol·C·m−2·yr−1; CO2 flux declined by approximately an order of magnitude along with declining river discharge. Episodic flooding made CO2 flux differ between dry (−0.7 to 20.9 mmol·C·m−2·d−1) and wet (11.6–170.0 mmol·C·m–2·d–1) conditions. During wet condition, increases in dissolved inorganic carbon (DIC) and total alkalinity (TA) significantly elevated CO2 degassing. Furthermore, ventilation of river‐borne CO2 strengthened degassing when estuaries became overwhelmingly river‐dominated. During flood relaxation, all estuaries experienced heightened productivity, evidenced by DIC and TA consumption in the mid‐salinity range (10–30). When prolonged drought led to hypersalinity (>36.5), biogeochemical and evaporative effects enhanced DIC and TA consumption and CO2 degassing. Due to flooding and high wind speeds, these estuaries were a strong CO2 source during spring and summer. Then they transitioned to a weak CO2 source or sink during the fall. Low temperatures further depressed CO2 efflux during winter except when a pulse of freshwater input occurred. This study demonstrates that changes in the hydrologic condition of estuaries, such as dry/wet cycle and river discharge gradient, will greatly alter air–water CO2 flux and estuarine contribution to the global carbon budget.

In the highly dynamic transition zone where the land meets the ocean, estuaries are generally net heterotrophic and therefore act as a CO 2 source to the atmosphere with substantial spatiotemporal variations in CO 2 flux (Laruelle et al. 2010;Chen et al. 2013). Coastal lagoons, a major type of estuary, are typically shallow (< 5 m) and have limited exchange with the adjacent ocean (Boynton et al. 1996). The estimated CO 2 flux (17.3 AE 16.6 molÁCÁyr −1 Ám −2 ) from coastal lagoons is almost the same as that from fjords (17.5 AE 14.0 molÁCÁyr −1 Ám −2 ), the world's largest estuarine type; together lagoons and fjords account for two-thirds of the total global estuarine area (the former accounts for 23.6% and the later accounts for 42.7%; Laruelle et al. 2010;Bauer et al. 2013). North America has 34% of the world's lagoons (Cromwell 1971), most of which are located in the Gulf of Mexico (GOM; Dürr et al. 2011). Yao and Hu (2017) suggested that Mission-Aransas Estuary in the northwestern GOM (nwGOM) is a CO 2 source with an estimated CO 2 flux at 12.3 AE 3.3 molÁCÁm −2 Áyr −1 . The lack of data for the rest of the nwGOM region makes it difficult to accurately quantify regional estuarine CO 2 flux, especially since the region is subject to a strong latitudinal climatic gradient.
There are seven major estuaries along the 600 km of Texas coastline in the nwGOM (Longley 1994). In each of these estuaries, there is a secondary bay (or upper estuary in this study) that directly receives freshwater inflow from at least one major river and a larger primary bay (or lower estuary in this study) that has restricted connection to the nwGOM due to the barrier island chain. Despite the similar geomorphic structure and physiography, these estuaries are remarkably hydrologically diverse due to the presence of a climatic gradient (Montagna et al. 2012). This gradient (wetter in northeast and drier in southwest) causes a difference in freshwater balance of almost two orders of magnitude, decreasing from northeast to southwest (Longley 1994;Montagna et al. 2012). Both recorded data (Texas Water Development Board, http://www.twdb.texas.gov) and a climate projection (Milliman et al. 2008) indicate that these estuaries are subject to prolonged periods of drought punctuated by periods of intense flooding, resulting in extreme changes in hydrologic conditions in relatively short periods of *Correspondence: hyao1@islander.tamucc.edu This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. time (weeks to months). The fast and intense hydrologic transitions make the nwGOM coast an ideal place to study hydrologic influence on estuarine CO 2 chemistry and flux.
To better understand the role of subtropical lagoons in regional and global air-water CO 2 flux budgets and to investigate the influence of the dry-wet cycle on the estuarine carbonate system under the context of global climate change, we characterized CO 2 chemistry and flux based on a 3-yr data set (2014)(2015)(2016)(2017) in four estuaries. We examined spatial and temporal variations in CO 2 flux and identified the primary controls on both carbonate chemistry and CO 2 flux.

Study sites
Four nwGOM estuaries ( Fig. 1)-Lavaca-Colorado Estuary (LCE), Guadalupe Estuary (GE), Mission-Aransas Estuary (MAE), and Nueces Estuary (NE)-were sampled between April 2014 and February 2017. These shallow estuaries (1-5 m in depth) are located in the middle of the nwGOM coast along a freshwater inflow gradient (Russell and Montagna 2007). We sampled LCE, GE, and NE seasonally in April, July, October, and January of each year from 2014 to 2017 (with the exception of 2017, when we sampled in February instead of January). These sampling months represent spring, summer, fall, and winter, respectively. Sampling occurred more frequently (biweekly to monthly) in MAE (Yao and Hu 2017). River samples were collected every other month between October 2015 and June 2017, with additional field campaigns immediately following flooding events. Carbonate chemistry of river end-members for these estuaries were then derived from the averages of corresponding riverine data. Ocean end-member was studied for all four estuaries, and its dissolved inorganic carbon (DIC) and pCO 2 were calculated using pH and total alkalinity (TA) from Texas Commission on Environmental Quality (https://www.tceq.texas.gov/) quarterly field campaigns.
Our sampling stations spanned from the river mouth(s) to tidal inlet in each estuary. Samples were taken from the surface and bottom of the water column in the estuaries (and from the surface only in the rivers) using a Van Dorn water sampler. Samples were preserved following the standard protocol for ocean CO 2 studies (Dickson et al. 2007). All field sampling occurred during the daytime. In-situ data (temperature, pressure, and dissolved oxygen [DO]) were obtained using a calibrated YSI 6600 V2 data sonde.
Monthly river discharge data were obtained from the U.S. Geological Survey (USGS) real-time streamflow record (http://waterdata.usgs.gov/tx/nwis/rt). Hourly wind speed and barometric pressure data were obtained from National Oceanic and Atmospheric Administration's (NOAA) weather stations along the coast (https://tidesandcurrents.noaa.gov/stations. html). Daily mean wind speed and barometric pressure were calculated for the sampling days. Wind speeds from anemometers that were typically installed at $ 3 m height above water were converted to wind speeds at 10 m height using the wind profile power law (Hsu et al. 1994).

Chemical analyses
Water samples were analyzed for carbonate system parameters-TA, DIC, and pH. Briefly, DIC and TA were both analyzed at 22 AE 0.1 C. DIC was determined by acidifying 0.5 mL of sample with 10% phosphoric acid and quantifying the extracted CO 2 on an AS-C3 DIC analyzer (Apollo SciTech). TA was determined using the Gran Titration (Gran 1952) on an AS-ALK2 alkalinity titrator (Apollo SciTech). Both DIC and TA analyses had a precision of AE 0.1%. Certified reference material (CRM Batch#142, 156, 159;Dickson et al. 2003) was used throughout our analyses for data quality control and assurance. pH was measured using either a spectrophotometric method (when salinity > 20, measured at 25 AE 0.1 C, Carter et al. 2013) or an Orion™ Ross pH electrode (when salinity < 20, measured at 25 AE 0.1 C, calibrated with National Bureau of Standards (NBS) buffers at 4.01, 7.00, and 10.01). pH measurements had a precision of AE 0.004 for the spectrophotometric method and AE 0.01 for the electrode. All potentiometric pH values from the electrode were converted from the NBS scale to the total scale. Salinity was measured using a benchtop salinometer calibrated with MilliQ water and a known-salinity CRM.
The partial pressure of CO 2 in the water (pCO 2,water ) and pH at in-situ temperature were calculated in the CO2SYS program (Lewis and Wallace 1998) using DIC and lab-measured pH at 25 C. We chose DIC and pH as the input variables to avoid possible errors associated with noncarbonate alkalinity when TA is used in speciation calculations (Abril 2015). Carbonic acid dissociation constants (K 1 , K 2 ) from Millero (2010) and the bisulfate dissociation constant from Dickson (1990) were used. In our previous study (Yao and Hu 2017), a salinity-dependent ΔTA (i.e., TA measured -TA calculated ) was observed. Orr et al. (2018) suggested 2.6-3.2% uncertainty in pCO 2 when pairing DIC and pH for calculations in CO2SYS, which equated to approximately an error of 8-16 μatm in calculated pCO 2 in this study (based on the annual average pCO 2 ). In addition, our calculated pCO 2,water matched well (AE 20 μatm) with real-time monitoring values using a calibrated SAMICO 2 sensor (McCutcheon et al. unpubl.).

Air-water CO 2 flux calculation
We used Eq. 1 to calculate the air-water CO 2 flux: where k (mÁd −1 ) is the gas transfer velocity as a function of wind speed, K 0 (molÁm −3 Áatm −1 ) is the gas solubility at measured in-situ temperature and salinity (Weiss 1974), and pCO 2,air can be calculated from: Yao et al. Estuarine CO 2 chemistry and flux where P b (atm) is the barometric pressure from NOAA weather stations ( Fig. 1), P w (atm) is the water vapor pressure calculated using salinity and temperature (Weiss and Price 1980), and xCO 2,air (ppm) is the mole fraction of atmospheric CO 2 in dry air. xCO 2,air data were obtained from a Mississippi-coast CO 2 buoy (https://www.pmel.noaa.gov/co2/story/Coastal+MS). Yao et al. Estuarine CO 2 chemistry and flux A positive flux (F value) means CO 2 degassing from the estuary to the atmosphere. Gas transfer velocity k was parameterized using wind speed and the equation from Jiang et al. (2008), which was derived from Raymond and Cole (2001) using more high wind speed (> 6 mÁs −1 ) data: where U (m s −1 ) is the wind speed at 10 m height, and Sc SST is the Schmidt number of CO 2 at in-situ temperature (Wanninkhof 1992). Area-weighted annual average CO 2 flux was calculated in each estuary following the equation: where F avg is annual average CO 2 flux (in mmolÁm −2 Ád −1 or molÁm −2 Áyr −1 ), F i is area-weighted CO 2 flux of sampling trip i, and d i indicates interval days in between two consecutive trips from i to i + 1 (Yao and Hu 2017).

Statistical analyses
According to river discharge (Fig. 2a) and corresponding estuarine salinity behavior (Table 1), the time periods were categorized into two hydrologic conditions (dry and wet). Two-way ANOVA was conducted to examine the response of the estuarine carbonate system to the dry/wet cycle and the difference between estuaries. A significant ANOVA model with significant interaction between the factors (dry/wet and estuary) indicates that the change in the dependent variable (the tested carbonate system parameter) in response to one factor depends on the level of the other factor. For example, if air-water CO 2 flux is the dependent variable and a significant interaction is found, then the effect of the dry/wet cycle on CO 2 flux varies significantly between estuaries. If a significant interaction was identified in a significant two-way ANOVA (p < 0.05), a follow-up Fig. 2. River discharge (a) and air-water CO 2 flux (b) in studied estuaries, shaded areas represent stand deviation within 95% confidence (unit: m 3 Ás −1 for river discharge, mmolÁCÁm −2 Ád −1 for air-water CO 2 flux). comparison of means was conducted for each level of a factor to remove the interaction. One-way ANOVAs were conducted to test for the dry condition only or the wet condition only, respectively, on whether the mean of the tested parameter varied between estuaries, and t tests were conducted to look for differences between wet and dry conditions in the means of the tested parameters within each estuary. The ANOVA assumptions of normality and homogeneity of variances were met, and a p value of < 0.05 was considered significant.

River discharge
During our study period, significant increases in river discharge were recorded in response to storms in spring-summer 2015, spring-summer 2016, and winter-spring 2017. Wet periods in LCE and GE corresponded to the average river discharge rates of 117.4 AE 111.2 and 69.6 AE 13.4 m 3 Ás −1 , respectively, which were about three times of those during dry periods (i.e., 41.9 AE 37.7 and 19.9 AE 9.3 m 3 Ás −1 , respectively). Although total discharge was much less in MAE and NE, these estuaries also experienced spikes in river discharge that could be categorized into the same dry and wet seasons, with increases from 0.5 AE 0.6 to 10.7 AE 11.6 and 4.9 AE 3.5 to 25.0 AE 34.5 m 3 Ás −1 , respectively. A freshwater inflow gradient is present across these estuaries as a result of changing watershed area and precipitation/evaporation balance (Table 3); during the course of our study, we observed an order of magnitude decline in freshwater inflow from northeast (annual discharge 116.4 AE 143.7 m 3 Ás −1 for LCE) to southwest (annual discharge 9.9 AE 19.2 m 3 Ás −1 for NE) estuary (Fig. 2a).

Salinity
Because of the variability in freshwater inflow balance, these estuaries exhibited significant temporal and spatial differences in salinity during our study period (Tables 1 and 2). Consistent with the freshwater inflow, mean salinity was lower in the north (19.1 AE 8.4 in LCE and 15.3 AE 9.6 in GE) and higher in the south (24.6 AE 9.6 in MAE and 31.2 AE 4.3 in NE). In addition, the fall season had the highest average salinity across all estuaries, possibly because of relatively low freshwater inflow and high evaporation. Hypersalinity (salinity greater than the ocean end-member salinity, 36.5) occurred in upper MAE and NE between summer and fall 2014, which marked the end of a multiyear drought.

DO concentration
No hypoxia (defined as DO < 62.5 μmol kg −1 ) was observed during our study period. DO was highest in the winter (316.2 AE 110.1 μmolÁkg −1 ) and lowest in the summer (207.1 AE 37.7 μmolÁkg −1 ) across all estuaries. In addition, DO decreased from northern (254.4 AE 88.5 μmol kg −1 annual mean in LCE) to southern estuaries (223.1 AE 33.9 μmolÁkg −1 annual mean in NE).

Estuarine carbonate system
TA and DIC varied greatly between seasons and between estuaries (Table 1, Fig. 3a,b). For example, GE had the highest average TA and DIC (3053.6 AE 354.3 and 2733.1 AE 391.7 μmolÁkg −1 , respectively), whereas LCE had the lowest average TA and DIC (2418.3 AE 444.4 and 2178.6 AE 388.5 μmolÁkg −1 , respectively). Both TA and DIC responded differently to the dry/wet cycle between estuaries, evidenced by the strong interaction of dry/wet and estuary effects in the two-way ANOVAs (p < 0.001, Table 2). Table 2. Two-way ANOVA tests examine mean differences in carbonate system parameters between estuaries and between dry/wet conditions. One-way ANOVAs test for differences between estuaries in parameter means in only dry conditions or only wet conditions, and t tests are for differences between dry/wet conditions within each individual estuary. F statistic means variance of the group means, p value means the probability that there is no difference between tested means in the model. ***p < 0.001; **p < 0.01; *p < 0.05.

Parameter
In LCE and MAE, TA and DIC decreased $ 200 μmolÁkg −1 after the flood. In GE, however, both TA and DIC increased $ 400 μmolÁkg −1 after the flood. In NE, neither of these parameters showed substantial changes (t test, p = 0.123 for DIC, p = 0.194 for TA; Table 2). The seasonal mean pH across our entire study area (including all stations across all four estuaries) was highest during the fall (8.133 AE 0.118) and lowest during the spring (8.007 AE 0.176). GE had a higher pH than other estuaries, with the highest seasonal pH (8.357 AE 0.234) observed in winter. In addition, the pH variation was significant under dry/wet cycle and responded differently between estuaries (both p < 0.001 from one-way and two-way ANOVAs; Table 2). For example, LCE and NE experienced large variations in pH between dry and wet conditions, with pH values decreasing after a flood by an average of 0.07 (LCE, p = 0.009 from t test) and 0.05 (NE, p = 0.030 from t test), respectively. In contrast, pH remained stable in GE (p = 0.310 from t test) and MAE (p = 0.397 from t test) despite changes in flooding condition.
Estuarine pCO 2 generally decreased from northern to southern estuaries, with the highest annual mean in LCE (595 AE 560 μatm) and lowest annual in NE (420 AE 107 μatm). Despite the large temporal variability in TA and DIC in each estuary, pCO 2 displayed a consistent seasonality throughout the studied estuaries (Fig. 3c); pCO 2 was relatively high in spring and summer (567 AE 428 and 585 AE 271 μatm, respectively) and lower in fall and winter (425 AE 96 and 390 AE 236 μatm, respectively). Again, variations of pCO 2 regarded to flooding condition displayed significant spatiality between estuaries (p < 0.001, one-way ANOVA; Table 2). From dry to wet condition, average pCO 2 almost doubled in LCE and GE (Table 1), even though the increase was statistically significant Fig. 3. Average seasonal variations of carbonate system parameters (DIC, TA, pCO 2 ) in studied estuaries, shaded area represents stand deviation within 95% confidence (unit: μmolÁkg −1 for DIC and TA, μatm for pCO 2 ).

Yao et al.
Estuarine CO 2 chemistry and flux (t test p < 0.001, Table 2) pCO 2 only increased $ 100 μatm in MAE while it showed no change in NE (t test, p = 0.930, Table 2). Thus, pCO 2 in LCE and GE had larger temporal variations than in MAE and NE (Fig. 3c).

Carbonate chemistry of river and ocean end-members
The river end-members displayed spatiotemporal variability in TA and DIC (Table 3). For example, TA and DIC concentrations were exceptionally high (3000-4500 μmolÁkg −1 ) during dry conditions at all river mouths and decreased substantially in wet conditions (both decreased 1500-3500 μmolÁkg −1 ). The only exception was NE, where minor changes in riverine TA and DIC (each $ 200 μmolÁkg −1 lower) were observed under wet conditions. Compared to the river end-members, the ocean endmember had little fluctuation in TA and DIC, with $ 200 μmolÁkg −1 decrease in wet conditions. River end-member pH was influenced by dry/wet cycle as well (p = 0.037 from t test, not shown), with the range significantly higher in dry (7.850-8.055) than wet condition (7.371-7.892). However, pH in the ocean end-member was only slightly greater (p = 0.500 from t test, not shown) in wet conditions (8.133 AE 0.058) than dry conditions (8.050 AE 0.071).

DIC and TA variations under dry and wet conditions
Our sampling period began in a drought period in 2014 with limited river discharge; in spring 2015 the drought was broken and the estuaries experienced a large increase in river discharge (Fig. 2a). For the remainder of the study duration, the entire area switched between dry and wet conditions intermittently.
Rivers in this semiarid area are known to be potentially important DIC and CO 2 sources (Butman and Raymond 2011;Zeng et al. 2011), and our study further illustrated temporal and spatial variations of riverine inputs (Fig. 4). Riverine DIC of LCE was 3174.5 AE 340.6 μmolÁkg −1 in dry conditions Table 3. Estuarine parameters and carbonate system of riverine and tidal inlet end-members for each estuary at dry (D) and wet (W) conditions. End-member carbonate system in each row indicates river end-member to corresponding estuary, tidal inlet is regarded as ocean end-member for all four estuaries. ( Table 3), close to the value reported in another study in this area (Zeng et al. 2011). All rivers were rich in DIC and TA compared to downstream estuaries and their ocean end-member, presumably due to high bedrock weathering and evaporation rates in this region (Zeng et al. 2011;Stets et al. 2014). However, during the wet period, intense flooding strongly diluted riverine DIC and TA (Table 3; Fig. 4b,f), this type of strong dilution effect was also observed in other studies in this region (Mooney and McClelland 2012;Montagna et al. 2018). This dilution effect varied significantly across the entire region, from the strongest in MAE (both TA and DIC were diluted more than half, $ 2800 μmolÁkg −1 ), to the weakest in NE (diluted by 200 μmolÁkg −1 ). Therefore, the contribution of riverine biogeochemistry to estuarine carbonate chemistry may vary between estuaries and between hydrologic regimes; this highlights the need for further research on river chemistry under meteorological and hydrologic controls. In estuaries, carbonate system variability is affected by the mixing of riverine and oceanic end-members as well as estuarine biogeochemical processes. Two end-member mixing models can be used to examine allochthonous and autochthonous dissolved constituents in order to analyze conservative/ nonconservative behaviors (Bianchi 2012). Estuarine DIC and TA clearly followed salinity change (in the range of 1000-2000 μmolÁkg −1 ; Fig. 4a-h), which reveals a clear mixing scenario from river mouths to tidal inlets, when considering (c) calculated DIC production/consumption in dry condition; (d) calculated DIC production/consumption in wet condition; (e) observed TA in dry condition; (f) observed TA in wet condition; (g) calculated TA production/consumption in dry condition; (h) calculated TA production/consumption in wet condition; (i) TA/DIC ratios in dry condition; (j) TA/DIC ratios in wet condition; (k) air-water CO 2 flux in dry condition; (l) air-water CO 2 flux in wet condition (unit: μmolÁkg −1 for DIC, TA, DIC estuarine , TA estuarine , mmolÁCÁm −2 Ád −1 for air-water CO 2 flux). the uncertainties from riverine inputs (standard deviations of riverine DIC and TA values were within the range of 150$550 μmolÁkg −1 for all estuaries; Table 3). Therefore, we followed the method in Jiang et al. (2008) to calculate DIC change due to river-ocean mixing: where (DIC r + o mix=i ) is DIC due to two end-members (river and ocean) mixing at station i; DIC river and Sal river are river endmember DIC and salinity; DIC ocean and Sal ocean are ocean endmember DIC and Salinity; Sal i is salinity at station i. Equation 5 is not applicable in estuaries under hypersaline conditions when the evaporation exceeds riverine input. Therefore, an evaporation-based equation was derived for the hypersaline water (S > 36.5): DIC i−1 also indicates surveyed data at station i from the previous field campaign, Sal i−1 is salinity at station i from the previous field campaign. Then produced/consumed DIC due to estuarine biogeochemical processes (DIC estuarine ) can be calculated as: Similarly, TA r + o mix , TA o mix , and TA estuarine can be estimated by simply replacing DIC with TA.
As a result of flooding, river-delivered organic matter was elevated, which subsequently enhanced estuarine respiration. Enhanced gross primary production due to increased nutrient input can also occur following flooding (Bruesewitz et al. 2013). Stimulated net ecosystem metabolism during wet periods contributed to large variation in DIC estuarine and TA estuarine (500-1000 μmolÁkg −1 , where positive values are net production and negative values are net consumption, Fig. 4d, h), especially in the upper estuaries. The positive DIC estuarine and TA estuarine values found for the majority of measurements during wet periods represented a heterotrophic state coupled with other biogeochemical processes. For example, ΔDIC vs. ΔTA (Δ represents the difference between values from two consecutive trips; Fig. 5) provided useful insights on drivers of DIC and TA dynamics (Sippo et al. 2016;. The majority of ΔDIC/ΔTA ratios in wet periods ( Fig. 5) were ratios that would be expected for a combination of aerobic respiration (−0.2), denitrification (0.9), and carbonate dissolution (2). All three processes have been well studied in this area during wet periods (Russell et al. 2006;Bruesewitz et al. 2013;Murgulet et al. 2018). This overall net heterotrophy also favored CO 2 emission (seen as positive air-water CO 2 flux in Fig. 4l). In addition, ΔDIC/ΔTA change in LCE and MAE indicated a strong dilution effect during a large flood (slope close to 1.2 when S < 5; Fig. 5), which was consistent with their riverine DIC and TA declines in wet conditions.
In contrast, a small number of negative DIC estuarine and TA estuarine near river mouths (S < 10, Fig. 4d,h) suggested a net consumption resulting from the transition to net autotrophy. Such consumption extended through the following period of flood relaxation. During dry conditions DIC and TA consumption could amount to a maximum of −2000 μmolÁkg −1 , particularly near the river mouth (Fig. 4c,g), and in these conditions, ΔDIC/ΔTA (Fig. 5) was mainly influenced by photosynthesis (−0.2), carbonate precipitation (2), and nitrification (−∞). Past studies in this region have shown high sustained levels of phytoplankton (Reyna et al. 2017) and particulate organic carbon (Mooney and McClelland 2012), indicating that autotrophy is favored. An and Joye (2001) also reported that nitrification is stimulated (10-fold higher, up to 12.85 mmolÁm −2 Ád −1 ) by benthic photosynthesis in a similar shallow estuary-Galveston Bay-adjacent to this region. It is reasonable to expect such enhanced nitrification occurred in the water column when phytoplankton accumulated after flood events. On the other hand, as all estuaries support commercial oyster production, carbonate precipitation may be an important process that could decrease both DIC estuarine and TA estuarine under dry conditions (Murgulet et al. 2018). Detailed explorations of C and N co-variation in response to hydrologic change are necessary to better understand the biogeochemistry at play rather than solely relying on stoichiometry.
Upper MAE and NE experienced hypersalinity (salinity exceeded the average ocean end-member salinity, S > 36.5; Fig. 4a,c,e,g,i,k) during May 2014 and October 2014. At hypersalinity CO 2 effluxes were variables but elevated (S > 35, Fig. 4k). This is consistent with our previous study, which revealed that evaporating estuarine water holds less dissolved CO 2 ; and CO 2 efflux is significant under high wind speeds (Yao and Hu 2017). This CO 2 efflux should contribute to DIC loss (DIC estuarine < 0; Fig. 4c). Consistent with a previous study that demonstrated net alkalinity loss at hypersalinity (Hu et al. 2015;Murgulet et al. 2018), we observed TA consumption with elevated salinity under dry conditions (Fig. 4e,  g). TA estuarine was expected to be negative in hypersaline conditions, as elevated ammonium release (up to 100 μmolÁm −2 Áh −1 ) and sulfide accumulation ($ 40 μmolÁkg −1 ) has been observed at the sediment-water interface in NE as a result of enhanced dissimilatory nitrate reduction to ammonium (Gardner et al. 2006). No hypoxia was observed in this study, indicating that nitrification and sulfide oxidation may have occurred in the water column. However, further studies that associate nitrogen and sulfur cycles with carbonate chemistry would be necessary fully reveal their roles.
River-borne, ocean-borne, and estuarine-generated CO 2 We    (Jiang et al. 2008): CO 2 ½ r + o mix=i is the aqueous CO 2 concentration if conservative mixing occurred between river and ocean end-members, calculated using DIC r + o mix and TA r + o mix :

Yao et al. Estuarine CO 2 chemistry and flux
[CO 2 ] estuarine is aqueous CO 2 due to estuarine biogeochemical reactions. [CO 2 ] i is actual aqueous CO 2 concentration based on field data.
As aqueous CO 2 concentrations were subject to water temperature changes, all [CO 2 ] categories were normalized to the average temperature (23.9 C) to eliminate the thermal effect (Jiang et al. 2008). TA in estuarine waters contains noncarbonate species, which would lead to an underestimation of CO 2 concentration if using TA and DIC as the input variables in CO 2 system speciation calculations; for example, including 50 μmolÁkg −1 (an average ΔTA from our surveys) of noncarbonate alkalinity in the TA input decreases calculated [CO 2 ] by 3.0 μmolÁkg −1 at salinity 23.3 and temperature 23.9 C (all average values). Nevertheless, this exercise would still be useful for qualitatively tracking the dynamic of CO 2 flux in the aquatic system, especially for [CO 2 ] mixing from multiple sources (Jiang et al. 2008).
Area-weighted average [CO 2 ] for each sampling trip was calculated (arithmetic mean of stations in each bay, multiplied by bay area), all three [CO 2 ] categories displayed significant temporal and spatial trends (Fig. 6, Table 2). When [CO 2 ] river accounted for more than 50% of [CO 2 ] the scenario was considered river-dominated, and when [CO 2 ] ocean accounted for more than 50% of [CO 2 ] the scenario was considered oceandominated. These estuaries displayed scenarios ranging between river-dominated and ocean-dominated, mostly depending on hydrologic conditions. For instance, river discharge was two orders of magnitude lower in NE and MAE than LCE and GE (Fig. 2a). This directly led to a [CO 2 ] river that was 9.2 times higher annually in GE (30.5 AE 9.1 μmolÁkg −1 ) than NE (3.3 AE 2.7 μmolÁkg −1 ) and a [CO 2 ] ocean 53.8% lower annually in GE (7.3 AE 4.0 μmolÁkg −1 ) than NE (15.8 AE 1.1 μmolÁkg −1 ). As a result, GE was river-dominated throughout the study duration (Fig. 6b), whereas NE was clearly Fig. 7. Correlation analyses between categorized aqueous [CO 2 ] x and air-water CO 2 flux in studied estuaries in wet and dry conditions. D, dry condition; W, wet condition (unit: μmolÁkg −1 for [CO 2 ] x ; mmolÁCÁm −2 Ád −1 for air-water CO 2 flux).
ocean-dominated throughout the study duration (Fig. 6d). NE likely maintained the ocean-dominated state due to anthropogenic intervention (i.e., dam construction and river fragmentation; Murgulet et al. 2016), which limit freshwater inflow and result in a more negative freshwater balance (Montagna et al. 2012). LCE and MAE exhibited both river-dominated and ocean-dominated scenarios depending on changing hydrologic conditions (Fig. 6a,c). The relative contribution of [CO 2 ] estuarine differed significantly between wet and dry conditions, while [CO 2 ] river only differed between wet and dry conditions in MAE and NE and [CO 2 ] ocean only differed between wet and dry conditions in NE (t test p < 0.05, Table 2).
Interestingly, all four estuaries responded differently to changes in flooding (Fig. 7). As expected, CO 2 flux was significantly related to [CO 2 ] river during flooding throughout the region (Fig. 7a-d). During the river-dominated state, CO 2 effluxes increased exponentially with elevated [CO 2 ] river in LCE, GE, and MAE (wet condition of Fig. 7a-c). A previous study suggested that river-borne CO 2 ventilation could contribute to such exponential increase of CO 2 emission in riverdominated estuaries (Borges et al. 2006). Borges et al. (2006) also suggested that river-borne CO 2 ventilation decreases exponentially with longer water residence, but even NE, where ocean-dominated condition prevailed most of the time, [CO 2 ] river still contributed to CO 2 flux in wet conditions (R 2 = 0.403, p = 0.038, Fig. 7d). However, relatively low correlation between [CO 2 ] river and CO 2 flux (for most estuaries R 2 < 0.2 both in dry and wet conditions; Fig. 7a-d) indicated that ventilation of river-borne CO 2 was not the main control on CO 2 flux across our study area.
[CO 2 ] estuarine , an indicator of estuarine biogeochemical processes, played a crucial role in determining CO 2 flux of this region (Fig. 7i-l). For example, CO 2 flux in MAE was primarily controlled by [CO 2 ] estuarine in both dry and wet conditions (R 2 = 0.801 and 0.960, respectively; p < 0.001 for both, Fig. 7k). Likewise, [CO 2 ] estuarine contributed significantly to CO 2 flux in all of our studied estuaries (p < 0.05 for all cases; Fig. 7i,j,l), especially during wet periods. Due to long residence times (40-360 d; Bianchi et al. 1999), pulses of freshwater inflow can generate profound influences in nwGOM estuarine systems, that is, CO 2 flux was elevated during flooding, but on an annual scale CO 2 flux was relatively depressed, potentially because nutrient input from floods stimulated phytoplankton growth. Both DIC estuarine and categorized [CO 2 ] revealed the strongest autotrophy during a mild drought after a major storm (DIC estuarine ranged −2000 to 1000 μmolÁkg −1 in salinity 10-20, Fig. 4c, and area-weighted [CO 2 ] estuarine ranged −15 to 45 μmolÁkg −1 in Fall 2015, Fig. 6). The autotrophy-dominated drought period yielded a similar CO 2 flux among LCE, MAE, and NE ($ 18.0 mmolÁCÁm −2 Ád −1 ), and it yielded a CO 2 sink (−0.7 AE 42.5 mmolÁCÁm − 2Ád − 1) in GE despite GE's high river discharge relative to MAE and NE.
Even though ocean-dominated scenarios were identified in LCE, MAE, and NE (Fig. 6a,c,d), the lack of correlation between [CO 2 ] ocean and CO 2 flux suggested that [CO 2 ] ocean exerted the weakest control on CO 2 flux in this area (Fig. 7e-h). The only exception was in LCE, where [CO 2 ] ocean was significantly correlated with CO 2 flux (p = 0.014 for dry and p = 0.002 for wet period; Fig. 7e). This exception may be caused by the decrease in organic matter with decreasing riverine input (seen as lower [CO 2 ] river ; Fig. 6a), which subsequently led to lower [CO 2 ] estuarine (Fig. 6a). [CO 2 ] estuarine was the most important contributor to CO 2 flux in LCE (R 2 = 0.646 and 0.952 for dry and wet conditions, respectively; Fig. 7i).
Annual and seasonal CO 2 flux During the 3 years of this study, the nwGOM estuaries acted as an overall CO 2 source to the atmosphere (Table 4; Fig. 2). By applying an area-weighted average (Eq. 4), annual air-water CO 2 flux from the entire system was Fig. 8. Seasonality of air-water CO 2 flux and categorized aqueous CO 2 in studied estuaries (unit: μmolÁkg −1 for [CO 2 ] x ; mmolÁCÁm −2 Ád −1 for air-water CO 2 flux).

Yao et al.
Estuarine CO 2 chemistry and flux 16.6 AE 17.1 molÁCÁm −2 Áyr −1 , with approximately an order of magnitude decline from the northeast to southwest (i.e., from 25.6 molÁCÁm −2 Áyr −1 in LCE to 2.7 molÁCÁm −2 Áyr −1 in NE, Table 4). This region-wide flux estimate agrees with an earlier estimation of CO 2 flux of global lagoonal estuaries (17.3 AE 16.6 molÁCÁm −2 Áyr −1 ; Laruelle et al. 2010). On the other hand, Chen et al. (2013) estimated that average CO 2 flux in North American estuaries is 2.2 molÁCÁm −2 Áyr −1 , although their study was mostly based on U.S. east coast estuaries. Furthermore, annual CO 2 flux in MAE during these 3 years was 6.9 AE 6.5 molÁCÁm −2 Áyr −1 ; our previous study estimated a much higher CO 2 flux in MAE (12.4 AE 3.3 molÁCÁm −2 Áyr −1 ), although that study only spanned May 2014 to April 2015, a period of drought that preceded a large freshwater discharge event in April 2015 (Yao and Hu 2017). In the longer time series that we present in this study, it is likely that elevated freshwater inflow delivered nutrients and subsequently enhanced autotrophic production and reduced CO 2 emission, especially when MAE returned to an ocean-dominated condition following the flooding. Moreover, wind speeds can play a crucial role in estuarine CO 2 flux. This windy and shallow environment seems to favor fast air-water gas exchange even though average estuarine pCO 2 was not much greater than atmospheric pCO 2 (Table 1). We want to note that we used a wind-dependent function that was originally derived from the open ocean to calculate gas transfer velocity (Eq. 3) (Wanninkhof 1992). There is not yet a consensus on the parameterization of estuarine gas transfer velocity, nor have there been direct measurements of this velocity in our study area; therefore, a better gas transfer velocity measurement is crucial to developing more accurate CO 2 flux estimates in future estuarine research (Raymond and Cole 2001;Jiang et al. 2008;Rosentreter et al. 2017). Regardless, the comparison between the two data series with different temporal coverages highlights the dynamic nature of estuarine environments. Not only did estuarine CO 2 flux in the nwGOM vary spatially, it also exhibited strong seasonality (Table 4; Fig. 8a). In this study, strong CO 2 efflux in spring and summer months weakened and transitioned to influx in fall and winter (Table 4; Fig. 8a). Large riverine inflows from three major freshwater discharge events and the simultaneous high wind speeds (data not shown) resulted in the strong CO 2 efflux in spring and summer. After the flood influence subsided, net autotrophy (as demonstrated by widespread high pH, low pCO 2 , and negative [CO 2 ] estuarine (Table 1; Fig. 6) significantly reduced the fall CO 2 efflux (Fig. 8a). Water temperature decreased from fall to winter along with pCO 2 and CO 2 flux. A temperature decline of $ 10 C could result in $ 200 μatm decrease in pCO 2 in MAE (Yao and Hu 2017). Together, the lowered pCO 2 and the lowest seasonal wind speed resulted in a moderate to weak CO 2 sink in these estuaries in the winter (Fig. 8a). A larger decline in winter pCO 2 was expected in the more northern estuaries LCE and GE since this thermal effect was also salinity-related, that is, temperature-normalized pCO 2 would decrease in low-salinity region with smaller temperature dependence (∂lnpCO 2 /∂T) (Joesoef et al. 2015). Therefore, the CO 2 flux reversal in winter appears to be controlled more by weather (temperature and wind speed) than biological activities. For example, CO 2 flux dropped $ 20 mmol CÁm −2 Ád −1 in GE and MAE from fall to winter (Table 4), yet [CO 2 ] river and [CO 2 ] estuarine were relatively the same (Fig. 8c,d). While the thermal effect became important in winter months, hydrologic change remained the most important to the system; this was demonstrated in LCE, where (despite high variations) elevated freshwater inflow ([CO 2 ] river increased from 8.3 AE 3.1 to 14.3 AE 6.5 μmolÁkg −1 ; Fig. 8b) and caused the seasonal CO 2 efflux to double from fall to winter (from 24.4 AE 44.5 to 53.5 AE 178.7 mmol CÁm −2 Ád −1 , Table 4).

Conclusions
Our study covered an extreme range of hydrologic conditions, from drought to flooding, and estuaries that ranged from river-dominated to ocean-dominated. In general, nwGOM estuaries LCE, GE, MAE, and NE are a net CO 2 source to the atmosphere on an annual scale with large temporal and spatial variations. About an order of magnitude decline in annual average air-water CO 2 flux was observed from northeast to southwest along the coastline. Substantial CO 2 degassing due to large freshwater inflow events and high wind conditions occurred in spring and summer, and CO 2 flux switched to a weak CO 2 sink in fall and winter. Hydrologic and meteorological influences, such as river discharge, wind speed, and water temperature played important roles in driving CO 2 fluxes. Both negative DIC estuarine and [CO 2 ] estuarine implied overall autotrophy after summer flooding, despite high initial river-borne CO 2 enrichment. In addition, CO 2 emission was elevated by evaporation and high wind speeds, which further led to net DIC consumption in hypersaline water. Overall, our findings indicate that estuarine carbon cycle variability is highly dependent on estuarine hydrologic condition, and more comprehensive studies should be done to further assess this effect in a broader context.