Responses of carbonate system and CO2 flux to extended drought and intense flooding in a semiarid subtropical estuary

Globally, estuaries are considered important CO2 sources to the atmosphere. However, estuarine water carbonate chemistry and CO2 flux studies have focused on temperate and high latitude regions, leaving a significant data gap in subtropical estuaries. In this study, we examined water column carbonate system and air–water CO2 flux in the Mission‐Aransas Estuary, a subtropical semiarid estuary in the northwestern Gulf of Mexico, by collecting samples at five System Wide Monitoring Program stations from 05/2014 to 04/2015. The carbonate system parameters (total alkalinity [TA], dissolved inorganic carbon [DIC], pH, CO2 partial pressure [pCO2], and carbonate saturation state with respect to aragonite [ΩAr]) and air–water CO2 flux all displayed substantial seasonal and spatial variations. Based on freshwater inflow conditions, a drought period occurred between 05/2014 and 02/2015, while a flooding period occurred from 03/2015 to 04/2015. Average DIC was 2194.7 ± 156.8 μmol kg−1 and 2132.5 ± 256.8 μmol kg−1, TA was 2497.6 ± 172.1 μmol·kg−1 and 2333.4 ± 283.1 μmol kg−1, pCO2 was 477 ± 94 μatm and 529 ± 251 μatm, and CO2 flux was 28.3 ± 18.0 mmol C·m−2·d−1 and 51.6 ± 83.9 mmol·C·m−2·d−1 in the drought and flooding period, respectively. Integrated annual air–water CO2 flux during our studied period was estimated to be 12.4 ± 3.3 mol·C·m−2·yr−1, indicating that this estuary was a net CO2 source. High wind speed, warm climate, riverine input, and estuarine biogeochemical processes all contributed to the high CO2 efflux despite the modest pCO2 levels year round.

Human activities have significantly increased atmospheric CO 2 concentration since the Industrial Revolution. Although occupying a small area (approximately 0.3%) of the global ocean, estuaries play a disproportionately important role in the global CO 2 budget (Bauer et al. 2013). In general, estuaries are a net CO 2 source due to net heterotrophy. For example, Frankignoulle (1998) suggested that CO 2 efflux from European estuaries represents 5-10% of anthropogenic CO 2 emissions throughout Europe. Global estuarine CO 2 emissions could reach an approximate rate of 0.1-0.25 PgÁCÁyr 21 , which is on the same order of magnitude as continental shelf CO 2 uptake and equivalent to as much as 30% of total riverine carbon export (Zhai et al. 2007;Cai 2011;Bauer et al. 2013;Chen et al. 2013;Regnier et al. 2013;Laruelle et al. 2015). Large CO 2 release from estuaries could be attributed to hydrologic conditions (i.e., due to higher dissolved inorganic carbon [DIC] to alkalinity ratio in river waters than the receiving seawater) and intensive biological activities. For example, Joesoef et al. (2015) found that more CO 2 is released into the atmosphere in the upper Delaware estuary than the lower estuary, and Guo et al. (2009) reported that CO 2 flux in the Pear River Estuary is dominated by aerobic remineralization of organic matter.
There are many uncertainties in estimating CO 2 flux in an estuary. A major reason for such uncertainty is temporal variation of riverine fluxes (Abril et al. 2004;Crosswell et al. 2014) that are not easily captured in estuarine studies, which tend to have low temporal sampling resolution. On the other hand, despite that lagoonal estuaries are generally recognized as important CO 2 sources particularly in tropical and temperate areas (Laruelle et al. 2010), whether the existing global estimate in those studied areas can represent all estuarine types is unknown given the drastically different hydrologic conditions of these estuaries. As one of the world's largest subtropical lagoonal estuary systems (D€ urr et al. 2011), the northwestern Gulf of Mexico (GOM) estuaries lack data for studying the CO 2 source/sink issue. Located in a semiarid subtropical region, south Texas has been experiencing prolonged drought with intense flooding occurring intermittently, thus estuaries in this area receive generally low riverine inflows, punctuated by large storms (Milliman et al. 2008;Mooney and McClelland 2012). How the changing hydrologic state alters CO 2 flux and its magnitude is unclear.
In an estuary, carbonate system parameters are usually controlled by mixing and biogeochemical processes. An estuary that receives variable nutrient and organic matter input from a river will have altered metabolic processes, which will affect the carbonate system (Doney et al. 2009;Feely et al. 2010). For example, in the Chwaka Bay of Tanzania, due to the presence of seagrass, pH increases and total DIC decreases through enhanced photosynthesis and calcification (Semesi et al. 2009); while in the Long Island Sound, significant pH reduction could occur when enhanced respiration in subsurface water is coupled with CO 2 production ). In estuaries with significant freshwater influence, carbonate saturation state with respect to aragonite (X Ar ) is strongly correlated with salinity. For example, in Glacier Bay in the eastern Gulf of Alaska, low total alkalinity (TA) concentration, resulting from glacier discharge, decreases X Ar substantially to below aragonite undersaturation (Reisdorph and Mathis 2014). The carbonate system parameters in high latitude estuaries have large fluctuations due to significant seasonal changes in freshwater flux, nutrient delivery, and light intensity (Cross et al. 2013;Reisdorph and Mathis 2014). Even in a single estuary, carbonate system parameters could change quickly because of varying riverine input and strength of biological activities, as well as weather conditions (Mooney and McClelland 2012).
In this study, we characterized carbonate chemistry and CO 2 flux in the Mission-Aransas Estuary (MAE). Our primary goals were to understand how the carbonate system in MAE responded to freshwater input, and to study the air-water CO 2 flux in this estuary and understand its control(s).

Field sampling
MAE is a shallow lagoonal estuary located along the south Texas coast in the northwestern GOM (Fig. 1). It consists of three connected water bodies: Aransas Bay is a defined as a primary bay and is connected to the GOM via the Aransas ship channel; Copano and Mesquite are secondary bays closer to riverine input (Kim and Montagna 2012). Copano Bay receives freshwater inflows directly from the Mission and Aransas Rivers, the two major freshwater sources for the entire estuary. Mesquite Bay receives inflow from adjacent San Antonio Bay during flooding periods.
Five long-term System Wide Monitoring Program (SWMP) stations were established in the MAE by the Mission-Aransas National Estuarine Research Reserve (MA-NERR) in 2007 ( Fig.  1, https://sites.cns.utexas.edu/manerr). The SWMP program is a nationally coordinated and standardized program and is used for tracking short-term variability and long-term changes in a host of biological, physical, and chemical parameters. The five stations were designed to represent key estuarine conditions that range from freshwater inflow (station CW in Fig. 1) to hydrologic connection between the estuary and the GOM (SC; Evans et al. 2015). Water depths at these sampling stations ranged from 1.2 m (MB) to 6.2 m (SC).
Surface and bottom water samples were taken by a Van Dorn water sampler during the 05/2014-04/2015 period following the standard protocol for ocean CO 2 studies (Dickson et al. 2007). From 11/2014 to 03/2015, samples were taken monthly, and sampling occurred biweekly during the warmer months in the rest of our sampling year. All field samplings were done between 09:00 h and 14:00 h. Briefly, 250 mL narrow-neck borosilicate glass bottles were used to collect water samples for TA, DIC, and pH analyses. 100 lL of saturated HgCl 2 was added to the water sample to arrest biological activity. The samples were stored at 48C in the dark until analysis, usually within 2-3 weeks of sample collection. One hundred twenty-five milliliters polypropylene bottles were used to collect Ca 21 samples. In the study by Bockmon and Dickson (2014), filtration for coastal water carbonate system characterization was recommended. However, we did not find significant difference between filtered and unfiltered samples in this estuary (also see Hu et al. 2015), thus we used unfiltered samples for this study. A calibrated YSI 6600 V2 data sonde was used to obtain in situ temperature, pressure, and dissolved oxygen (DO) concentration at each station.
To examine the effect of freshwater inflow on the estuarine carbonate system, we imported daily discharge data from the United States Geological Survey (USGS) real-time stream flow record (http://waterdata.usgs.gov/tx/nwis/rt) for Mission (USGS #08189500) and Aransas Rivers (USGS #08189700), then calculated monthly riverine discharge. Wind speed and barometric pressure were obtained every 30 min from the weather station at Copano East (CE) (http:// lighthouse.tamucc.edu/pq), and daily mean wind speed was applied to sampling days. Wind speed data collected from the weather station (3 m) was converted to 10 m above the water surface using the wind profile power law, here u 2 is the wind speed at height z 2 5 10 m, u 1 is the collected wind speed data at height z 1 5 3 m, the exponent P (0.11) around GOM area is extracted by Hsu et al. (1994; p 5 0.11).

Chemical analyses
All water samples were analyzed for DIC, TA, pH, and salinity. Ca 21 was analyzed for surface water only. For DIC analysis, 0.5 mL water sample was acidified by 0.5 mL 10% H 3 PO 4 using a 2.5 mL syringe pump. The released CO 2 was analyzed on an AS-C3 DIC analyzer (Apollo SciTech). To analyze TA, 25 mL water sample was titrated with a 0.1 M HCl solution (in 0.5 M NaCl) using an AS-ALK2 alkalinity titrator (Apollo SciTech). Temperature of the titration vessel was maintained at 22 6 0.18C using a water-jacketed circulation system. Certified Reference Material or CRM (Dickson et al. 2003) was used to construct the standard curve for the DIC analysis and to calibrate the acid used for TA titration. DIC and TA analyses had a precision of 6 0.1%.
A spectrophotometric method (Carter et al. 2013) using purified m-cresol purple (mCP) obtained from Dr. Robert Byrne's lab (University of South Florida) was used for pH analysis on total scale (Liu et al. 2011). The indicator was adjusted to pH 7.92 6 0.01 every time before sample analysis with the aid of a calibrated Orion TM Ross TM glass electrode. A 10-cm water-jacketed absorbance cell used for pH measurements (Carter et al. 2013) was kept at 25 6 0.018C. The dye effect was corrected via duplicate runs of each sample by adding two volumes (30 lL and 60 lL) of mCP following the procedure in Clayton and Byrne (1993). This method had a precision of 6 0.0004 pH unit. Because of salinity limitations of the spectrophotometric method (20-40, Liu et al. 2011), for lower salinity samples we used the calibrated pH electrode to measure pH at 258C. All pH values obtained using the potentiometric method were converted to total scale using temperature and salinity (Millero 2001).
Salinity was measured using a benchtop salinometer (Orion Star TM A12, Thermo Scientific), which was calibrated using MilliQ water and known salinity CRM seawater each time before sample analysis. The salinometer was also regularly calibrated with 0.5 M KCl at 258C. Ca 21 concentration was determined by potentiometric titration (Kanamori and Ikegami 1980) using EGTA as the titrant. The end-point was detected using a Metrohm V R calcium-selective electrode on a

Yao and Hu
Carbonate chemistry and CO 2 flux in an estuary semi-automated titration system. This method had a precision of 6 0.5% for estuarine waters.
Water pCO 2 and X ar at field conditions were calculated using the program CO2SYS (Lewis and Wallace, 1998) based on DIC and lab measured pH. We used measured pH and DIC as input variables instead of the pH/TA combination to calculate pCO 2,water to avoid possible errors that were resulted from organic alkalinity component (Abril et al. 2015). In fact, we found that the difference between titration alkalinity and calculated alkalinity (with pH and DIC) showed strong salinity dependence in the MAE (Fig. 2). Although calculated pCO 2,water was subject to some degree of uncertainty, previous studies showed an 1 : 1 linear relationship between measured and calculated pCO 2,water over a range of 300-4000 latm in estuarine and coastal waters (Frankignoulle and Borges 2001).
Calculated X ar from CO2SYS output was corrected using measured Ca 21 concentration, which was a near perfect linear function of salinity throughout the year (r 2 5 0.99, data not shown). Carbonic acid dissociation constants (K 1 , K 2 ) in Millero (2010), bisulfate dissociation constant in Dickson (1990) were used in this calculation. In the MAE (http://missionaransas.org/science/download-data), the soluble reactive phosphorous concentration was on the order of $ 2.0 lmol kg 21 or less, and silicate concentration was on the order of $ 200 lmol kg 21 . The effect of nutrients on calculated X ar and pCO 2 was minimal. For example, pCO 2 would only differ by 0.2 latm with or without nutrients in the CO2SYS program. Because not all water samples that were characterized for carbonate chemistry had concurrent nutrient information, nutrients were omitted in all carbonate speciation calculations.

Air-water CO 2 flux calculation
The air-water flux of CO 2 was calculated using the following equation: F5kK 0 pCO 2;water 2pCO 2;air À Á where k (mÁd 21 ) is the gas transfer velocity calculated from wind speed, K 0 (molÁm 23 Áatm 21 ) is the gas solubility at measured in situ temperature and salinity (Weiss 1974), pCO 2,water and pCO 2,air are partial pressure of CO 2 in surface water and the atmosphere, respectively. Positive F value means CO 2 degassing to the atmosphere. pCO 2,air can be calculated from: Here P b (atm) is the barometric pressure, which was downloaded from the weather station at Copano East (CE), 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 atmospheric CO 2 in dry air. We did not measure air xCO 2 directly but chose to download monthly averaged xCO 2 data from http://www.esrl.noaa.gov/gmd/ ccgg/trends. We recognized the spatial difference in xCO 2 on the global scale. However, compared to the xCO 2 data from a coastal CO 2 monitoring site in eastern GOM (https://www. pmel.noaa.gov/co2/story/Coastal 1 MS), monthly average xCO 2 is no more than 15 ppm greater in the GOM than the Mauna Loa site in winter and the two xCO 2 records generally agree with each other to within 6 2 ppm in summer. Therefore the subsequently calculated CO 2 flux would differ by a few percent using either xCO 2 record (or the actual xCO 2 at MAE). Given that the GOM site does not have a continuous dataset, we used the Mauna Loa data for our calculation.
Gas transfer velocity (k) would differ depending on wind speed, tidal current, and bottom topography (Wanninkhof 1992;Raymond and Cole 2001). Unfortunately, there is no widely accepted k formulation in shallow estuary regions, and the carbon cycle community still has to rely on wind speed dependence to estimate gas exchange rates. Because the equation in Raymond and Cole (2001) was mainly based on relatively low wind speed (< 7 mÁs 21 ) estuarine data, here we used the equation in Jiang et al. (2008), which covered larger amount of high wind speed data (up to 12 mÁs 21 ) and was more appropriate to our studied area: where U is the wind speed at 10-m height (mÁs 21 ), Sc SST indicates Schmidt number of CO 2 at the in situ temperature from the freshwater (flooding period) and seawater (drought period) equations, respectively (Wanninkhof 1992). This equation has also been adopted in recent estuarine studies (Bozec et al. 2012;Crosswell et al. 2012).
The total surface area of Aransas and Copano bays is 452 km 2 based on Texas Water Development Board (TWDB) record. Aransas Bay occupies 181 km 2 with estimated 10% of its area in the Aransas Ship Channel where the station SC is located, and Copano Bay is 271 km 2 . There is no published record on the area of Mesquite Bay. We estimated its area using a closet rectangle on Google Earth V R and calculated its area as 75 km 2 . Due to relatively limited sampling stations, we decided to use area-weighted average method to calculate CO 2 flux in the MAE instead of taking arithmetic mean of the individual stations. Our results showed that the difference between these two methods for each trip was 2.5 6 7.0 mmolÁCÁm 22 Ád 21 , and integrated annual CO 2 flux is 12.4 molÁCÁm 22 Áyr 21 , compared with 11.4 molÁCÁm 22 Áyr 21 using the arithmetic mean. We first calculated area-weighted CO 2 flux using CO 2 flux values at the five sampling stations and the respective areas above. Then to estimate seasonally (drought and flooding periods) or annually averaged CO 2 flux, we used the following equation: F avg is area-weighted CO 2 flux and has a unit of mmolÁC m 22 Ád 21 or molÁC m 22 Áyr 21 , F i is air-water CO 2 flux of each sampling trip, d i indicates days in between two consecutive trips. Note in our CO 2 flux calculations we did not consider diel effect. Given that our sampling always took place during the middle of the day (09:00-14:00 h) when pCO 2 was likely the lowest due to photosynthesis, calculated fluxes may represent lower estimates of the actual values ).
Thermal and non-thermal effects on pCO 2 variations Thermal and non-thermal effects on pCO 2,water variations were evaluated using the method in Takahashi et al. (2002, Eqs. 5-9). To remove the temperature effect, we normalized pCO 2,water to annual mean temperature of 23.08C (T mean , Eq. 6). In surface seawater pCO 2 doubles for every 168C increase in the oceanic waters (@lnpCO 2 /@T 5 0.0423; Takahashi et al. 1993). However, the average @lnpCO 2 /@T of MAE water was calculated to be slightly lower (0.0411 6 0.0068) than open ocean water. To evaluate the effect of temperature change on pCO 2 , we used Eq. 7.
where d is @lnpCO 2 /@T, subscripts mean and obs stand for annual mean and observed values, respectively. To understand relative contributions of thermal and nonthermal effects on temporal pCO 2 changes, we applied the following equations: T=B5DpCO 2;thermal =DpCO 2;non2thermal (10) Both DpCO 2,thermal and DpCO 2,non-thermal were calculated from the difference between the maximum and the minimum pCO 2 effects during either the drought or the flooding period. The T/B ratio illustrated the relative importance of thermal vs. non-thermal effects. In MAE, non-thermal effect indicated the combination of biogeochemical processes and physical mixing. Thermal effect on surface pCO 2,water would exceed non-thermal if T/B ratio was greater than 1; conversely, non-thermal effect would dominate if T/B ratio was less than 1.

Data analysis
A two-factor analysis of variance (ANOVA) was used to test the effects of hydrologic condition and sampling locations on carbonate variables and CO 2 fluxes. Probabilities (p) of < 0.05 were considered as significant. Normality and homogeneity of variance were ensured before ANOVA was conducted and there was no need to transform the data.

Hydrologic conditions
During our sampling year, the drought period (05/2014-02/2015,18 trips in total) was much longer than the flooding period (03/2015-04/2015, four trips in total; Fig. 3). Annual discharge for Mission River was 8.29 310 7 m 3 Áyr 21 while discharge for Aransas River was 1.88 310 7 m 3 Áyr 21 . Aransas River had more discharge during the drought period (07/ 2014-02/2015), in which the month with the least discharge was recorded in August 2014 for Aransas and Mission Rivers (97.40 3 10 3 and 0 m 3 month 21 , respectively). However, Mission River had more freshwater input during the flooding period (03/2015-04/2015). The MAE was affected by two strong storms in late 03/2015 and 04/2015 (Fig. 3). During this period, Mission River discharge was almost twice as much as Aransas River. The highest discharge from Aransas River (51.25 m 3 s 21 ) was recorded in 03/2015 and discharge from Mission River was the highest (118.93 m 3 s 21 ) in 04/ 2015.
The average values of all physical and chemical parameters for the five stations during the drought and the flooding periods are listed in Table 1 and Fig. 4. The average water temperature during our study period was 23.0 6 6.38C (N 5 216), with the highest temperature (average 29.  N 5 9), which was the highest during our study period. On the contrary, during the flooding period, there was a salinity gradient (surface and bottom average) across the estuary from the lowest (18.9 6 10.4, CW; N 5 8) to the highest (29.0 6 3.2, SC; N 5 8). Salinity stratification was only observed in Copano Bay. Salinity decreased sharply at CW after the flood (from 30.1 to 5.1 for surface water).

Carbonate chemistry
All carbonate system parameters displayed significant spatial (station) and temporal differences (between the drought    (N 5 216, Fig. 5). Average DIC concentration was 2194.7 6 156.8 lmolÁkg 21 (N 5 176) during the drought period, and was lower from summer to early fall (2106.8 6 100.6 lmolÁkg 21 , 06/2014-10/2014; N 5 109) but higher in winter (2230.6 6 96.0 lmolÁkg 21 , 11/2014-02/ 2015; N 5 40). CW had the highest DIC concentration and SC had the lowest in drought condition. DIC concentrations during the flooding period (with an average of 2132.5 6 256.8 lmolÁkg 21 ; N 5 40) were lower than the drought period. There was a large decrease in the surface waters at CW after the first storm event in late 03/2015 with DIC concentrations decreasing from 2239.4 lmolÁkg 21 (prior to the storm) to 1227.3 lmolÁkg 21 (after the storm). During the flooding period, MB had the highest DIC concentration, CW the lowest. Along with salinity stratification in Copano Bay during the flooding period, DIC concentration was lower in surface waters than the bottom water.
For the entire estuary, average pCO 2,water was 487 6 138 latm (N 5 216, Fig. 5). It was 477 6 94 latm (N 5 176) during the drought period, and 529 6 251 latm (N 5 40) during the flooding period. During the drought period, the highest average pCO 2,water occurred at CE (542 6 105 latm, N 5 36) and lowest at SC (423 6 70 latm, N 5 35). During the flooding period, pCO 2,water in most stations increased, and CW had the highest average pCO 2,water , and SC had the lowest. At the same time, calculated pCO 2,air in MAE area during our study period averaged at 391 6 6 latm using observed monthly xCO 2,air at the Mauna Loa site.
Overall, annual average air-water CO 2 flux was 12.4 6 3.3 molÁCÁm 22 Áyr 21 (N 5 110, Fig. 6), and the entire estuary released CO 2 during most of our sampled period. Further, the air-water CO 2 flux pattern displayed significant seasonal and spatial variations (p < 0.05 tested by two-way ANOVA; Fig. 6). Throughout the drought period, air-water CO 2 flux varied between 213.5 6 7.6 mmolÁCÁm 22 Ád 21 and 109.4 6 71.3 mmolÁCÁm 22 Ád 21 (02/2015 and 08/2014, respectively; N 5 5 for each average value). The estuary was a CO 2 source in the drought period during warm months (05/ 2014-10/2014) and was a sink during cold months (11/ 2014-02/2015). Moreover, there was an overall increase in CO 2 efflux from the drought to the flooding period mostly due to large increase of CO 2 emission at CW. In particular, there was a drastic increase (from 3.4 mmolÁCÁm 22 Ád 21 to 380.3 mmolÁCÁm 22 Ád 21 ) at CW in late 03/2015 right after the first storm event.

DIC and TA dynamics during the dry-wet cycle
Many Texas rivers have high levels of bicarbonate ion (HCO 2 3 ) as a result of high weathering rates of the drainage basins (Zeng et al. 2010) and generally high evaporation Table 2. Temporal and spatial analysis of the carbonate system tested by two-way ANOVA. "ͱ" denotes a significant main factor interaction between sampling time (i.e., drought-flooding cycle) and sampling location. rates in the area (http://www.twdb.texas.gov/surfacewater/ conditions/evaporation/). To unravel the factors that contributed to estuarine DIC and TA variations, we used the water chemistry data to construct two-endmember mixing diagrams for DIC and TA in the MAE during the drought and flooding periods, respectively (Figs. 7,8); we also investigated the effect of precipitation and evaporation. Note for the river endmembers, we did not take river endmember during the 2014-2015 boat trips. Since we did not have concurrent river chemistry data during our estuarine sampling period, to best estimate the river chemistry during the flooding and drought conditions, we collected river samples at Mission and Aransas rivers bimonthly between December 2015 and December 2016 (Table 3). Six trips were conducted during a drought period, and one trip was conducted at the end of May 2016 after significant flooding in south Texas. As discussed in Hu et al. (2015), along a river-ocean mixing line, the lowest solute to salinity ratio (i.e., the slope of the evaporation line) would be at the ocean endmember if there were no reactions that consumed this solute. Clearly, all the data points to the right of the dotted line (i.e., the precipitation-evaporation line for seawater endmember, Fig.  7) during the drought period reflected net removal of DIC (panels a, b, and c) and TA (panels d, e, and f), and the Copano Bay stations (CE and CW) also showed larger extent of DIC and TA consumption under hypersaline conditions (panels a and d). Furthermore, because the Mission River had higher DIC and TA concentrations than the Aransas River (Table 3), the mixing line between the Mission River (the dashed line in Fig. 7) and the ocean water endmember (values in Table 3) should have a slightly steeper slope than that between the Aransas River and the ocean endmember (the solid line in Fig. 7). Therefore, all the data points that are bracketed by the river-ocean mixing lines and the left of the precipitation-evaporation line (to the lower of the plot) also indicated removal.

Parameter
Normally in an estuarine mixing zone, data points that appear above a mixing line would indicate in situ production as many estuarine studies have indicated (Raymond et al. 2000). However, for lagoonal estuaries with prolonged residence time under drought conditions ), significant evaporation could increase solute concentration and salinity simultaneously based on the original river-ocean mixing line (see Hu et al. 2015 for a detailed discussion). Therefore, it is likely that evaporation may have played an important role in the water chemistry during the drought period (Fig. 7a,d). Our results could not rule out in situ production of DIC and TA in Copano Bay, although as one moves to Aransas Bay, consumption clearly dominated DIC and TA changes (Fig. 7b,e).
In the literature, alkalinity reduction has not been commonly reported in coastal estuaries (other than in coral reefs). Instead, alkalinity production due to net anaerobic reactions (pyrite burial and net denitrification) is often suggested as an important process Smith et al. 1991;Hu and Cai 2011). Alkalinity consumption due to calcification clearly is a possible mechanism as the MAE has abundant oyster reefs, representing the southernmost location for viable commercial production . While not directly observed in the MAE, sulfate concentration in the adjacent Corpus Christi Bay immediately to the south of the MAE appeared to show excess relative to a conservative mixing during the drought period (D. Murgulet pers. comm.), indicating possible external reduced sulfur contribution that was oxidized to sulfuric acid. In addition, Benoit et al. (1994) also observed high levels suspended matter that contains Fe in San Antonio Bay and Corpus Christi Bay, both of which "bracket" the MAE. It is known that oxidation of reduced sulfur and iron would produce H 1 that titrates TA (Chen 2002). Given the shallow depth of the MAE and windy conditions in this area, significant benthic contribution (i.e., reoxidation of reduced sedimentary compounds) is likely, especially during high wind conditions that could cause abundant sediment resuspension. More detailed studies, such as examining sulfur and metal dynamics under different hydrologic conditions, are needed.
During the flooding period, DIC and TA concentrations in the river endmembers decreased (Table 3) presumably due to the dilution effect. Apparently, it became more difficult to explain the DIC and TA variations using the simple mixing lines derived from the local rivers. However, a closer examination of the CW station data suggested that both surface (TA and DIC) and bottom waters (TA only) showed excellent mixing behavior while DIC in the bottom water showed slight respiration signal during the 2-month floodwater-dominated period. For example, intercepts of TA-salinity regressions were 1207 lmol kg 21 (r 2 5 0.95) for surface and 1449 lmol kg 21 (r 2 5 0.99) for bottom water respectively; those of DIC-salinity regressions were 1123 lmol kg 21 (r 2 5 0.98) for surface and 1374 lmol kg 21 (r 2 5 0.89) for  bottom water, respectively (Fig. 8a,d, regression lines not shown). Therefore, the stagnant endmember values (Table 3) were likely not accurate during the high freshwater inflow season, and the different regression intercepts between the surface and bottom waters may reflect changing freshwater endmember compositions (Cifuentes et al. 1990). Furthermore, Mesquite Bay (MB) (surface and bottom) and the Aransas Bay (bottom) showed "excess" DIC and TA compared to the river mixing line (Fig. 8b,c,e,f). This was likely caused by overflow of high alkalinity San Antonio River (up to 5500 lmol kg 21 , Hu, unpublished data) water through Mesquite Bay (Evans et al. 2015) during the flooding period. This high alkalinity water influence decreased from Mesquite Bay to Aransas Bay.
Controlling factors on the pCO 2,water variations during the dry-wet cycle Carbonate system speciation in estuarine waters is controlled by various factors such as temperature, biological processes (primary production, calcification, etc.), and riverocean mixing (Bozec et al. 2012;Hu and Cai 2013;Hunt et al. 2014;Joesoef et al. 2015). Here we evaluated the effect of temperature and riverine inputs during the drought and flooding periods.
During the drought period, the mean amplitude of thermal effect on water pCO 2 (pCO 2,thermal -pCO 2,water ) was $ 190 latm, which was partially compensated ($ 100 latm) by counteractive non-thermal effect (pCO 2,non-thermal -pCO 2,water ), thus a net $ 100 latm of seasonal amplitude on pCO 2,water was observed when there was a warming effect from spring to summer or a cooling effect from summer to winter (Fig. 9). During the flooding period however, mean thermal effect declined to $ 90 latm as water temperature increased slightly in 03/2015, whereas there was a large fluctuation of nonthermal effect in 03/2015-04/2015. In particular, sharp increases in pCO 2,water was observed in Copano Bay (CW and CE, respectively, Fig. 9), and pCO 2,non-thermal reached its highest concentration (1211 latm) when the first storm came in late 03/2015. These increases indicated that the storm event caused a dramatic increase in non-thermal pCO 2 (biological and/or mixing). According to Bruesewitz et al. (2013), community respiration in Copano Bay would greatly increase following a storm (thus increasing pCO 2 ). In addition, high river inflow resulted from storm event also brought in high pCO 2 water. The appearance of peak pCO 2,non-thermal (715 latm, early 04/2015) was delayed at CE compared to CW (Fig.  9) and the value was lower, indicating that the non-thermal effect became less pronounced with freshwater inflow propagating along the estuary.
Overall, thermal effect dominated the drought period from 05/2014 to 02/2015 for all stations except in AB (Table  4), where thermal and non-thermal effects were about the same. However, during the flooding period, Copano and Mesquite Bay had much lower T/B ratio, indicating that non-thermal effect played a dominant role in controlling pCO 2,water variation, consistent with the discussion above. In addition, AB and SC exhibited similar T/B ratios across the dry/wet cycle, indicating that hydrologic state probably did not have a significant effect on the relative importance of thermal vs. non-thermal effect, at least during our sampled time. Considering that station AB is located in the primary bay that has direct exchange with the GOM through the ship channel, and that both AB and SC had smaller salinity variations during the annual cycle (Table 1), these stations clearly were less influenced by river inflow but more by exchange with the GOM.
Air-water CO 2 flux dynamics and controlling factors during the dry-wet cycle Although accounting for 41% of world's estuarine area, average CO 2 flux in North American estuaries accounts for only 12% of global estuarine CO 2 emission, at a moderate value of 2.2 molÁCÁm 22 Áyr 21 , according to a recent synthesis (Chen et al. 2013). However, the average annual air-water CO 2 flux from MAE reached 12.4 6 3.3 molÁCÁm 22 Áyr 21 based on our study. Even though this flux was consistent with results obtained in a limited number of tropical lagoons   (Laruelle et al. 2010), it was significantly higher than many other North American estuaries and even some macrotidal estuaries in Europe (Table 5).
As an important factor that determines gas transfer velocity (Eq. 4), wind speed plays an important role in estuarine CO 2 flux. According to Chen et al. (2013), the mean airwater pCO 2 gradient of European estuaries is only about one third of those in Asian estuaries, whereas the mean air-water CO 2 flux from European estuaries doubled that in Asia due to much higher wind speed. Similarly, despite a relatively low average pCO 2,water (487 6 138 latm) in MAE compared with much higher pCO 2,water in European and Asian estuaries (1600 latm and 4000 latm, respectively; Chen et al. 2013), high mean wind speed (5.4 6 2.3 mÁs 21 , compared with approximately 4 mÁs 21 and 1.6 mÁs 21 on European and Asian coasts, respectively; Chen et al. 2013) contributed to a relatively high CO 2 efflux in this estuary. Therefore, it is desirable to further investigate the role of other low latitude regions, not only within the estuaries, but coastal ocean, in CO 2 budget calculations.
Similar to many other estuaries, air-water CO 2 flux in MAE displayed strong temporal changes during our studied period (Tables 2 and 6). During the drought period, CO 2 emission from the estuary was the highest in summer and early fall (06/2014-09/2014, Fig. 6). This is likely a result of high wind speed (average 6.0 6 2.4 mÁs 21 ) and high water temperature, and the later could enhance community respiration and evaporation, as more concentrated seawater holds less CO 2 . For example, a simulation using CO2SYS suggests that for seawater with salinity 35, total alkalinity 2270 lmol kg 21 , total DIC 1977 lmol kg 21 , pCO 2 of this seawater would be 400 latm at 258C. Then allowing the water evaporate to salinity 40 (hypersaline condition in Copano Bay), pCO 2 would increase to 494 latm (Fig. 3). If allowing all * Two-way ANOVA suggests that both "dry-flooding cycle" and "stations" had significant effect on CO 2 flux, and the two main effects also had significant interaction (df 5 4, F 5 2.538, p < 0.05). "excess" CO 2 to degas and reach equilibrium with the 400 latm atmosphere, a degassing process that may take only a few days given the shallow water depth here, this evaporation-concentrated water would hold $ 50 lmol kg 21 less DIC than the "concentrated" original seawater. With additional alkalinity reduction, regardless its cause, prolonged water residence time would lead to even more CO 2 loss to the atmosphere. Indeed, TA/DIC ratio in all stations increased with increasing salinity (i.e., average TA/DIC increased from 1.116 6 0.012 to 1.172 6 0.024, along with the increase of average salinity from 31.6 6 1.0 to 38.7 6 0.8, N 5 10; Fig. 10). Part of the TA/DIC increase can be attributed to higher TA/DIC ratio in the ocean endmember, and then further depletion of DIC relative to TA toward hypersaline conditions indicated that evaporation contributed to a net CO 2 loss. Toward the end of the drought period, pCO 2,water started to decline with decreasing temperature, and the entire estuary became a weak CO 2 sink from late fall to winter (11/ 2014-02/2015). Particularly in 02/2015, at all five stations surface water was undersaturated for CO 2 (average pCO 2,water 5 341 6 49 latm). This could be attributed to low riverine input and lower water temperature (Hunt et al. 2014). Note this period also had lower wind speed (average 3.7 6 0.9 mÁs 21 ) thus CO 2 uptake was modest.
CO 2 efflux also displayed spatial variations during our sampling period (Table 6). There was a decreasing trend of average CO 2 emission (46.2-8.1 mmolÁCÁm 22 Ád 21 ) from the secondary bays (Copano and Mesquite) to the primary bay (Aransas). Mean air-water CO 2 flux at the ship channel (SC) was the lowest (8.1 6 34.2 mmolÁCÁm 22 Ád 21 ; N 5 22). This spatial distribution agrees with other estuarine studies as CO 2 efflux typically decreases toward the ocean (Crosswell et al. 2012;Hunt et al. 2014;Joesoef et al. 2015). However, higher CO 2 emission was observed at mid-estuary (CE and MB) especially in the drought period. This was possibly due to more intense remineralization reactions (Table. 1; Fig. 5), which was in favor of CO 2 production. In flooding season, non-thermal effect dramatically increased (Table 4) air-water CO 2 flux at CW, whereas CO 2 degassing was maintained at similar level or even decreased in other part of this estuary, due presumably to nutrient-enhanced primary production (Reyna et al. 2017).

Conclusions
Both carbonate chemistry and CO 2 flux demonstrated substantial temporal and spatial variations in the subtropical MAE, which was affected by strong interannual changes in estuarine hydrologic states. There was a gradient for carbonate variables and CO 2 flux from the secondary bays to the primary bay. We observed significant TA and DIC removal during the drought period and mixing dominated distribution in the flooding season, although detailed mechanisms for the alkalinity loss still await further investigation.
Overall, the MAE was a CO 2 source with an annual average air-water CO 2 flux 12.4 6 3.3 molÁCÁm 22 Áyr 21 . High wind speed played an important role for this high CO 2 efflux despite the relatively small air-water pCO 2 gradient. This estimate is much higher than existing, yet scarce, estimates in other subtropical estuaries. During the drought and warm period, CO 2 emission was enhanced by increased temperature (hence increased respiration) and evaporation, and highest CO 2 emission (74.5 6 41.1; N 5 30 mmolÁm 22 Ád 21 ) was found in summer in the drought period. However, lower temperature in winter would turn the entire estuary into a weak CO 2 sink (28.9 6 5.1 mmolÁm 22 Ád 21 ; N 5 15). In the flooding period, storm events briefly yet significantly enhanced air-water CO 2 flux at CW due to much elevated pCO 2 level there. Overall, our work suggests that global estuarine CO 2 flux estimates need to be improved by incorporating new studies that focus on subtropical and/or windy estuaries.