Eutrophication Amplifies the Diel Variability of Carbonate Chemistry in an Equatorial, Semi-Arid, and Negative Estuary

Introduction

Monthly average global atmospheric CO₂ concentration reached 418.9 ppm in May-2021 (UCSD-SIO, 2021), which is the highest concentration in the past 3 million years (Willeit et al., 2019). Considering this dramatic increase of atmospheric CO₂, it is essential to improve carbon budget estimations at local/regional and global scales, identifying sources and sinks of this major anthropogenic greenhouse gas (IPCC, 2021). Estuarine ecosystems play a disproportional role in coastal carbon budget. They occupy a modest global area (0.2% of global ocean), but contribute with a significant global CO₂ emission on the order of 0.10 Pg C yr^–1 (Chen et al., 2013); these emissions respond to ∼5% of the annual ocean CO₂ sink (Friedlingstein et al., 2019). However, the controls of carbonate chemistry and air-water CO₂ flux in estuaries are complex and still not well understood and quantified (Borges, 2005). Estuarine ecosystems are highly dynamic transitional areas, with diverse “coastal typologies” and high spatial-temporal variabilities of biogeochemical properties (Dürr et al., 2011).

Estuaries are generally considered sources of CO₂ to the atmosphere because they exhibit net heterotrophic metabolism. It means the rates of community respiration (both autotrophic and heterotrophic) are higher than the rates of gross primary production (Borges and Abril, 2011). Furthermore, estuaries receive CO₂-rich waters from riverine discharges (Jiang et al., 2008; Van Dam et al., 2018; Cotovicz et al., 2020a), and lateral inputs of DIC and TA from vegetated coastal ecosystems (mangroves, saltmarshes) (Ovalle et al., 1990; Santos et al., 2021). Most studies concerning air-water CO₂ exchanges and carbonate chemistry were conducted in temperate and subtropical regions, and generally in turbid estuaries dominated by significant river discharges (Frankignoulle et al., 1998; Borges and Abril, 2011). Studies of carbonate chemistry are overlooked in tropical ecosystems, particularly in semi-arid regions on and near the equator. In such regions, the river flow changes drastically during the year depending on the balance between evaporation and precipitation (Lavín et al., 1998; Dias et al., 2009). Furthermore, the significant water retention upstream in artificial reservoirs has been described in several estuaries impacted by dams (Dias et al., 2009; Mulligan et al., 2020). Consequently, the turbidity maximum zone, rarely reported in studies of estuarine systems in northeastern Brazil, can be retained by tidal forcing and favor oxygen consumption, with modifications in organic matter processing and the export of carbon to the ocean (Dias et al., 2016; Cavalcante et al., 2021). When the evaporation exceeds the freshwater supply from rivers and forms hypersaline waters the estuaries exhibit negative or inverse circulation (Lavín et al., 1998; Dias et al., 2009). Hypersalinity occurs under prolonged drought (low humidity) and warm climate, leading to changes in the physico-chemical properties and carbonate chemistry of estuaries with implications considering the magnitude of air-water CO₂ flux (McCutcheon et al., 2019; Yao et al., 2020). Hypersalinity occurs in many coastal ecosystems; however, investigations of the carbonate chemistry in inverse estuaries are almost non-existent.

In addition, there is a lack of information regarding the temporal variability of carbonate chemistry and associated CO₂ fluxes in estuaries, mostly at diel time scales (day—night changes). The diel cycles in solar radiation exert a periodicity on biogeochemical processes, creating diel (that is, 24-h) patterns of aquatic CO₂ concentrations (Gómez-Gener et al., 2021). The diel-CO₂ variability in general is governed by the photosynthetic activity occurring at daytime with assimilation of DIC from the water (mainly CO₂); however, during night-time the photosynthesis is interrupted, and the aquatic CO₂ concentrations increases due to the microbial respiration of organic matter. Despite this obvious process, there is a lack of CO₂ measurements especially at night-time. The diel variably has showed to be important in diverse estuarine typologies and climate domains, including subarctic, temperate, subtropical, and tropical coastal regions (Yates et al., 2007; Dai et al., 2009; Bozec et al., 2011; Cotovicz et al., 2015; Fairchild and Hales, 2021; Miller and Kelley, 2021). The diel variability affects all carbonate chemistry parameters, including concentrations of DIC, TA, CO₃^2–, and values of pH, Ω_cal, and Ω_ara (Cyronak et al., 2018; Fairchild and Hales, 2021; Miller and Kelley, 2021). The diel patterns of CO₂ system in equatorial coastal regions is particularly neglected.

Tropical coastal ecosystems have been suffering with increasing nutrient pollution and eutrophication due to the insufficiency of wastewater treatment facilities. The eutrophication has been associated with changes in coastal carbon budgets (Borges and Gypens, 2010; Cai et al., 2011; Cotovicz et al., 2015). However, the response of coastal ecosystems to eutrophication is strongly site-specific, in which some ecosystems develop acidification of subsurface waters (high pCO₂, low pH, CO₃^2–, Ω_cal, and Ω_ara), whereas others are the opposite and exhibit an increase in pH (low pCO₂, high pH, CO₃^2–, Ω_cal, and Ω_ara) (Borges and Gypens, 2010; Cai et al., 2011; Cotovicz et al., 2015, 2021).

The Jaguaribe River Estuary (JRE), located at an equatorial region in the northeastern coast of Brazil, is classified as a well-mixed estuary (Dias et al., 2009). The reduced rains, in addition to the constructions of large dams in the river basin and high average atmospheric temperatures (∼28°C), have decreased the freshwater flux into the estuary. The very low freshwater supply and the high rates of evaporation contribute to the insignificant riverine discharge during dry periods. For this reason, the salinities during dry periods are higher or similar to those of adjacent coastal waters (Dias et al., 2009, 2016). The hydrochemistry in the JRE reflects the seasonal variability of the semi-arid climate, and the estuary behaves as a retainer of DIC in the dry season (Cavalcante et al., 2021). However, the carbonate chemistry dynamics in the estuary was not addressed. The JRE suffers with increasing eutrophication, particularly through discharges from shrimp aquaculture ponds and domestic effluents that impact local mangroves and the main estuarine channel (Eschrique et al., 2014; Marins et al., 2020; Lacerda et al., 2021). This tropical semi-arid estuary shows rapid environmental change related to anthropogenic disturbances, including eutrophication, river damming, decreasing of freshwater discharge and hypersalinity. The main objective of the present study was to investigate the diel variability of carbonate chemistry in the estuary, with emphasis on the quantification of air-water CO₂ exchanges. Our main hypotheses are that: (1) the concentrations of TA and DIC in the estuary will be higher compared to the adjacent coastal ocean due to the high rates of evaporation, the heterotrophic metabolism, and the supply of DIC and TA by mangrove forests and anthropogenic effluents; (2) The diel variability of carbonate chemistry parameters will present both semi-diurnal (tidally-driven) and diel (biologically driven) tendencies, with minimum buffering capacity at low tide and night-time conditions.

Materials and Methods

Study Area

The JRE basin (4° 23′ S and 37° 43′ W; 4° 36′ S and 37° 43′ W) has a surface area of about 1,350 km², located in a tropical region near the equator (Dias et al., 2009; Figure 1). The watershed of the Jaguaribe River that covers 74,327 km² with an extension of about 610 km contributing to the Western Equatorial Atlantic Ocean (Dias et al., 2009; Cavalcante and Cunha, 2012). Historically, the freshwater discharge into the estuary ranged from 0 to 7,000 m³ s^–1 (Campos et al., 2000). However, the freshwater inputs have decreased due to the construction of a series of dams at the Jaguaribe River watershed (Dias et al., 2009). Historical monthly rainfall shows a marked seasonal behavior, with a wet period between March and April (200–400 mm), and a dry period between August and November, when the precipitation can frequently be zero. Between 2012 and 2017, a persistent extended drought (<700 mm.yr^–1) dominated the river basin climate, and in 2018 annual rainfall started to return to the historical average (FUNCEME, 2021). Indeed, there is a decreasing continental runoff in the semi-arid northeastern region of Brazil as the result of decreasing annual rainfall, and this alarming scenario of decreasing freshwater supply have worsened due to river damming (Lacerda et al., 2020). The water residence time during the dry season is longer in the higher/middle estuary with average of 3 days, and up to 13 days (Lacerda et al., 2013). The estuary has a semi-diurnal and meso-tidal regime with average tidal height of 2.8 m, with peaks during spring tide reaching up to 3.4 m (Dias et al., 2009). The human population in the estuarine basin is estimated at about 100,000 inhabitants, whereas the wastewater treatment services comprise less than 40% of the households (IPECE, 2017). This has contributed to increasing levels of nutrients in the estuary, and occasional occurrence of hypoxia particularly in mangrove channels of the upper/middle estuarine regions (Eschrique et al., 2014; Marins et al., 2020). The trophic status in the estuary, using the trophic state index (Lamparelli, 2004; Silva, 2019), was classified as mesotrophic to eutrophic during the sampling period of this study. In addition to urban wastewaters, discharges from about 3,640 ha of shrimp farms significantly increases nutrient loads into the estuary (Lacerda et al., 2021). The estuary is surrounded by mangroves, which occupy about 13,000 ha (Godoy et al., 2018).

Sampling Surveys and Analytical Procedures

Two eulerian time series were performed, one in October-2017 (33 h of sampling), and other in September-2018 (44 h of sampling) (Figure 1). Both samplings were performed during spring tides and dry season, and at the same moored sampling station. The monthly averaged precipitation during these two surveys were null, and the water budget was negative, i.e., evaporation exceed rainfall rates (Figures 1B,C). The moored station was chosen because it represents a transitional region in the central estuary, which receives influences from continental and marine intrusion and in dry season often represents a turbidity maximum zone, where organic matter interactions and DIC retention occur (Dias et al., 2016; Cavalcante et al., 2021).

Continuous Measurements

The collection campaign included continuous and discrete sampling. The continuous and real-time measurements were performed for the partial pressure of CO₂ (pCO₂), salinity, temperature, and wind velocity. The continuous measurement system is based on Pierrot et al. (2009), and well-described in Carvalho et al. (2017) and Cotovicz et al. (2020b). Briefly, a water pump was placed at a depth of ∼0.5 m and provided continuous water flow (∼2.5 L min^–1) to the boat. The water flow was directed to a thermosalinograph (SeaBird Electronics ^®) to record the surface temperature and salinity. After, passing through the thermosalinograph, this water flux was directed to two equilibrators (shower head type) and then discharged. These showerhead equilibrators promote fast equilibration between air and water phases inside the system. The gas, free from humidity, passes through a Non-dispersive InfraRed gas analyzer (NDIR) for CO₂ quantification (Licor-7000 ^®CO2/H2O gas analyzer). A data acquisition system determined the following parameters every 5 min: date and time, position of the ship, velocity of the ship, molar fraction of CO₂ in the equilibrator (xCO₂), water content in the detector, sea surface temperature (SST) and sea surface salinity (SSS). The molar fraction of CO₂ (xCO₂ ppm) measured in the equilibrator in dry gas was computed using the equation below:

where Peq is the pressure in the equilibrator (assumed to be the same of atmosphere), and P_weq is the pressure of water vapor (atm), according to Weiss and Price (1980). The temperature measured in the surface water and in the equilibrator were slightly different, and then a correction was applied to compensate such difference, according (Takahashi et al., 1993):

where pCO₂ represents the seawater pCO₂ at in situ conditions, SST is the temperature measured in situ and Teq is the temperature measured in the equilibrator.

The NDIR was calibrated prior the first analyses and after every 6 h of seawater molar fraction records using air free CO₂ by passing N₂ standard followed by standard gases with nominal concentrations of 360, 1,009, and 2,009 ppmv (99.9% purity, supplied by White Martins Certified Gases).

Atmospheric pCO₂ measurements (pCO_2air) were performed every 6 h with air taken from the top of the vessel at ∼10 m high. The accuracy of the pCO₂ measurements was estimated at ± 2 μatm. One anemometer model Davis S-WCF-M003 was used to measure the wind velocity and placed at about 10 m in height.

Discrete Water Sampling

Discrete water samples, at a depth of ∼0.5 m, were collected hourly using a 3-L Niskin bottle. The water samples were filtered in the boat using Whatman GF/F filters (diameter 0.47 mm, pore size 0.7 μm), which were used for chlorophyll a (Chl a) analysis, and the filtrate used for nutrients (phosphate) and TA analysis. All filters were pre-combusted (at 500°C for 6 h). The filters and the filtered water samples were conditioned (fixed and/or maintained in ice in the dark) for further analysis in the laboratory. At the laboratory, the Chl a was extracted in 90% acetone and determined by spectrophotometry following the procedures described by Jeffrey and Humphrey (1975). The phosphate (PO₄^3–) was determined by colorimetric method according to Hansen and Koroleff (1983). TA was determined on 60 mL of filtrate using the Gran (1952) electro-titration method with an automated titration system (Mettler Toledo model T50). The reproducibility of TA was about 3 μmol kg^–1 (n = 7). Measurements were compared to certified reference material (CRM, provided by A. G. Dickson from Scripps Institution of Oceanography) and consistent at an accuracy level of ± 5 μmol kg^–1.

Carbonate Chemistry Calculations

The pH (at the National Bureau of Standards scale, NBS), DIC, Ω_cal, and Ω_ara were calculated from pCO₂, TA, seawater temperature, and salinity using the CO2calc 1.2.9 program (Robbins et al., 2011). The dissociation constants for carbonic acid were those proposed by Mehrbach et al. (1973) refitted by Dickson and Millero (1987), the borate acidity constant from Lee et al. (2010), the dissociation constant for the HSO₄^– ion from Dickson (1990) and the CO₂ solubility coefficient of Weiss (1974). The Ksp values for aragonite and calcite were taken from Mucci (1983) and the concentrations of calcium (Ca²⁺) were assumed proportional to the salinity variations according to Millero (1979).

Mixing Model

According to Jiang et al. (2008), in estuaries and lagoons with low or negligible freshwater inputs, the concentrations of DIC and TA during estuarine mixing can be calculated as:

Where SSS_measured is the measured surface salinity, SSS_ocean the surface salinity of the ocean endmember, DIC_ocean and TA_ocean the DIC and TA concentrations of the ocean endmember. This model assumes that the DIC_conservative is the DIC concentration after the ocean endmember is diluted by a zero DIC freshwater. Here, we applied this same approach; however, we assumed that the ocean endmember is linearly concentrated by the evaporation considering the salt conservation. This same approach was recently applied by Cotovicz et al. (2021) studying an evaporative coastal lagoon. The salinities in the inner shelf off the Jaguaribe River during dry season are relatively constant (∼36.7); typical variabilities are < 1 in the ocean endmember (Dias et al., 2013).

The deviation from conservative mixing (ΔDIC) is defined as the DIC addition or loss relative to the theoretical DIC concentration during mixing:

Where DIC_measured is the measured DIC. In the same way, the deviation of TA from the evaporation path (ΔTA) can be calculated. As we did not measure the ocean endmember, we use the value from Cotovicz et al. (2020b), which measured TA and DIC in an adjacent coastal region with similar salinity and water temperature and during the dry season.

Air-Water CO₂ Fluxes

The air-water CO₂ fluxes (FCO₂) was calculated according:

where k is the gas transfer velocity coefficient (cm.d^–1), K₀ (mol.cm^–3.atm^–1), is the solubility coefficient of CO₂ at in situ temperature and salinity (Weiss, 1974), and pCO_2water and pCO_2air (pressure) are the partial pressures of CO₂ in equilibrium with surface water and in the overlaying air, respectively.

The gas transfer velocity was parameterized as a function of wind speed taking account two parameterizations available for estuaries (Raymond and Cole, 2001, RC01; Jiang et al., 2008, J08), and two specific parameterizations for oceanic waters (McGillis et al., 2001, M01; Wanninkhof, 2014, W14). The parameterizations of RC01, M01, J08, and W14 can be calculated as follows:

Where U represents the measured wind velocity (m s^–1), Sc is the Schmidt number for CO₂, and 660 is the Schmidt number of CO₂ in seawater at 20°C.

Net Community Production

The NCP was calculated considering the changes in dissolved inorganic carbon (DIC) with time, comparing consecutive peaks of maximal and minimal tidal heights. For highest tidal heights, the NCP_HT was calculated according:

Where NCP_HT is in (mmol m^–2 h^–1), ρ is the seawater density (kg m^–3), d is the average depth (m) of the area, t represents the time interval (12 h), and FCO₂ is the averaged carbon dioxide flux (mmol m^–2 h^–1) across the water–atmosphere interface during this period. DIC_1HT and DIC_2HT represent the salinity-normalized concentration of DIC (mmol kg^–1) during two consecutive measurements at the highest tide. In this manner, we can compare consecutive DIC measurements at the maximal tidal level and in different diel period (daytime and night-time). Furthermore, we calculated the NCP for the lowest tidal heights (NCP_LT) according:

Where DIC_1LT and DIC_2LT represent the salinity-normalized concentration of DIC (mmol kg^–1) during two consecutive measurements at the lowest tide. The averaged values of NCP_LT and NCP_HT provide the NCP of the estuary.

Statistics

The Shapiro-Wilk test was applied to check whether a given variable follows a parametric or a non-parametric distribution. All investigated variables in this study were not normally distributed; therefore, we only applied non-parametric statistics. To compare two unpaired groups, we calculated the Mann-Whitney test. To compare three or more unmatched groups, we calculated the Kruskal-Wallis test. Spearman rank coefficient was calculated to assess the statistical correlation between the rankings of two variables. All statistical analysis were based on significance level of 0.05. We used the GraphPad Prism 8 program (GraphPad Software, Inc., La Jolla, California) to perform graphs and statistical tests.

Results

Main physico-chemical parameters analyzed in this study are shown in Table 1. The carbonate chemistry parameters presented marked variability at the semi-diurnal and diel time scales (Figures 2, 3). Both sampling campaigns presented similar tidal heights and amplitudes, with maximal tidal height of 3.60 m and minimum tidal height of 0.10 m. The averaged tidal amplitude was 3.3 m in Oct-2017, and 3.0 m in Sep-2018. The salinity showed a strong and negative trend with tidal height, showing characteristics of an inverse estuary (Figure 3). Higher salinities (max = 40.8) were measured during low tides (p < 0.01), whereas lower salinities (min = 38.2) were measured during high tides. The averaged salinity was 39.3 ± 0.7 in Oct-2017 and 39.4 ± 0.2 in Sep-2018, which are higher than the values found in the adjacent coastal region. The water temperature followed a clear diel pattern, with warming from dawn to midday (highest temperature of 29.1°C) and cooling from midday to dusk/night (lowest temperature of 26.7°C) (Figure 3). The averaged water temperature was 28.2 ± 0.6 and 28.0 ± 0.6°C for Oct-2017 and Sep-2018, respectively, with no significant differences considering the two time-series.

The TA variability followed a clear semi-diurnal pattern, but with an inverse trend compared to the tidal height (Figures 2A,B). The highest concentrations were measured during lowest tidal heights (up to 2,724 μmol kg^–1), whereas lowest concentrations were measured during highest tides (down to 2,440 μmol kg^–1) (p < 0.01). TA averaged concentrations were not different considering daytime and night-time periods (p > 0.05). DIC variability also followed a semi-diurnal tendency, but with diel influences (Figures 2C,D). Highest concentrations were verified during lowest tidal heights and night-time (max = 2,450 μmol kg^–1), whereas lowest concentrations were verified during highest tidal heights and daytime (min = 2,145 μmol kg^–1) (p < 0.01). The red arrows in Figures 2C,D show the difference of DIC concentrations between consecutive peaks of maximal tidal heights occurring in daytime and nigh-time conditions, illustrating the influence of diel variability (Figure 4). On Oct-2017 presented DIC concentration averaging 2,272 ± 76 μmol kg^–1 that is about 50 μmol kg^–1 lower than the averaged DIC concentration on Sep-2018. DIC and TA concentrations presented positive correlation with salinity (Figure 5). The Spearman correlation coefficient between TA and salinity was 0.98 and 0.97 on Oct-2017 and Sep-2018. Considering the correlation between DIC and salinity these coefficients were 0.97 and 0.95, respectively, on Oct-2017 and Sep-2018. This positive correlation with salinity was also verified for PO₄^3– concentrations, with minimum concentration of 0.09 mg L^–1 (high tide and daytime) and maximum concentration of 0.40 mg L^–1 (low tide and night-time). The correlation coefficient (Spearman rank) between PO₄^3– and salinity was 0.96 for both sampling campaigns. Measured TA and DIC concentrations are overall higher than those predicted by the evaporation model, creating positive deviations in the values of ΔTA and ΔDIC (Figures 5A,B). Considering the mean value of ΔTA = 48 μmol kg^–1, ΔDIC = 107 μmol kg^–1, the surface area in the middle estuarine portion = 450 km², the mean water depth = 5 m, and considering the residence time of 3 days, the estimated net gains of TA and DIC are about 80.8 and 179.2 mmolC m^–2 d^–1, respectively. Figure 5C present unitless directional vectors representing the slopes of the main processes affecting TA and DIC, showing that data points are close to the slopes of sulfate reduction and denitrification.

Concentrations of CO₃^2–, Ω_cal, and Ω_ara presented a similar tendency, exhibiting a combination of semi-diurnal and diel influences (Figure 2). Highest values of CO₃^2–, Ω_cal, and Ω_ara were verified during low tide at daytime, with values reaching up 251 μmol kg^–1, 5.79 and 3.88, respectively, during sampling in Oct-2017. The lowest values of CO₃^2–, Ω_cal, and Ω_ara were verified during high tide at night-time, with values of, respectively, 185 μmol kg^–1, 4.30 and 2.86. Overall, on Oct-2017 the values of CO₃^2–, Ω_cal, and Ω_ara were higher than those on Sep-2018. The red arrows in Figure 2 shows the difference of CO₃^2–, Ω_cal, and Ω_ara values comparing consecutive peaks of maximal tidal heights occurring at different time of the day (daytime and night-time). The values of CO₃^2–, Ω_cal, and Ω_ara were higher during daytime compared to night-time (Figure 4) (p < 0.001). Chl a concentrations exhibited values averaging 2.55 ± 0.97 μg L^–1 in Oct-2017, and 1.98 ± 0.98 μg L^–1 in Sep-2018. No significant trend was verified considering flood and ebb tides; however, a significant difference was verified considering the diel variability: highest concentrations were measured at daytime.

Following the tendencies described above for TA, DIC, CO₃^2–, Ω_cal, and Ω_ara, the values of pCO₂ also exhibited marked differences considering the hour of sampling (diel variability) and tidal height (semi-diurnal) (Figures 2, 4). On Oct-2017 high pCO₂ values occurred at low tide and night-time with pCO₂ values reaching up to 731 μatm, and a maximal diel amplitude of 219 μatm. On Sep-2018 we also found this same tendency, with values reaching the maximum value of 860 μatm during low tide and night-time, and maximal diel amplitude of 231 μatm. On Oct-2017 the averaged pCO₂ value (716 ± 45 μatm) was slightly above than the average measured in Sep-2018 (632 ± 49 μatm). The values of pCO₂ measured during daytime were lower than those measured during night-time, in both sampling campaigns (Figure 4) (p < 0.001). The red arrows in Figure 2E, f illustrate the diel variability of pCO₂ values considering consecutive peaks of maximal pCO₂ occurring during low tide, but with different magnitude depending on the hour of sampling (daytime × night-time). As expected, the pH followed a significant and inverse trend compared to pCO₂. Highest values of pH were verified during daytime and high tide (7.96), whereas lowest values were verified during night-time and low tide (7.83), with diel pH amplitude reaching a maximum of 0.10 units (NBS scale).

Table 2 shows the main parameters used to calculate the air-water CO₂ fluxes (FCO₂). The wind velocities were similar for both sampling campaigns, however, with some particularities. Oct-2017 presented wind velocities averaging 4.15 ± 0.99 and 3.96 ± 1.59 m s^–1 for daytime and night-time, respectively, with no statistical difference between these periods. Sep-2018 presented wind velocities averaging 4.73 ± 0.74 m s^–1 for daytime that was higher than the average of 3.47 ± 1.28 m s^–1 measured during night-time (p < 0.0001). This diel pattern verified for wind velocities was the same of that verified for the values of gas transfer velocities (k₆₆₀). The parameterization of RC01, which is specific for estuaries, provided the highest values of k₆₆₀ that averaged 10.15 cm h^–1 on Oct-2017, and 10.21 cm h^–1 on Sep-2018 (Figures 6C,D). The parameterization of W14, which is specific for oceanic waters, provided the lowest values of k₆₆₀ with a mean of 5.58 cm h^–1 in Oct-2017, and 5.65 cm h^–1 in Sep-2018. Considering the averaged values, the wind velocities and gas transfer velocities did not present differences between Oct-2017 and Sep-2018. The FCO₂ were always positive, indicating that the estuary is a source of CO₂ to the atmosphere in all conditions (Table 2 and Figure 6). Following the tendency verified for k₆₆₀ values, the CO₂ emissions were highest applying the parameterization of RC01 (mean of 67.1 mmolC m^–2 d^–1), and lowest with the parameterization of W14 (mean of 37.3 mmolC m^–2 d^–1), whereas the parameterizations of M01 and J08 provided intermediate values (Table 2 and Figure 6). On Oct-2017, the diel difference of wind velocities and k₆₆₀ values were not statistically different. Therefore, the highest FCO₂ values were verified during low tide and night-time, following the trend verified for pCO₂ values. The emissions during night-time (48.1 mmolC m^–2 d^–1) were higher than during daytime (42.5 mmolC m^–2 d^–1). On Sep-2018, however, the wind velocities and k₆₆₀ were lower during night-time and pCO₂ values were higher. For this, despite the fact that night-time presented higher pCO₂ values, the values of FCO₂ were lower during night-time (54.1 mmolC m^–2 d^–1) compared to daytime (67.5 mmolC m^–2 d^–1). Considering all FCO₂ data and all parameterizations, the averaged CO₂ emission was estimated at 51.9 ± 26.7 mmolC m^–2 d^–1. As we could not account for current velocity, which can be an important driver of k₆₆₀ in shallow estuaries, our evasion rates should be considered conservative (Abril et al., 2009; Jeffrey et al., 2018).

The net heterotrophic metabolism in the JRE was estimeted at about −5.17 ± 7.39 mmol C m^–2 h^–1, varying from +6.86 to −14.66 mmol C m^–2 h^–1 (Table 2). Most NCP values were negative; however, there are few periods of autotrophy as indicated by the positive NCP. For each individual tidal cycle, we compared consecutive peaks of DIC concentrations from low-to-low tides occurring at different hour of the day (daytime and night-time). The same was calculated for consecutive peaks of high tides (high-to-high). The maximal heterotrophy is verified when we compared the DIC concentrations from daytime to night-time (low-to-low and high-to-high tides), considering that DIC concentrations (and pCO₂ values) were always higher at night-time. The two situations of autotrophy occurred when we compared the DIC concentrations from night-time to day time (low-to-low and high-to-high tides).

Discussion

Production of Dissolved Inorganic Carbon and Total Alkalinity in the Estuary

Changes in carbonate chemistry in semi-arid, hypersaline estuaries are still overlooked (McCutcheon et al., 2019; Yao et al., 2020). Hypersalinity occurs frequently in the JRE during dry months because evaporation exceeds water inputs (Dias et al., 2013). This favors water retention by tidal forcing which increases water residence time in the estuary and enables sediment trapping in the middle estuary (Dias et al., 2013). Water in the middle estuary is blocked for long periods before reaching the sea (Lacerda et al., 2020). Negative water balance is reflected in higher salinities in the estuary (averaging 39.5) compared to the adjacent coastal ocean (<36.5; Cotovicz et al., 2020b; 36.45 ± 0.24 along the equatorial coast; Carvalho et al., 2017). This corroborates previous results in this estuary that exhibits negative estuarine circulation (Dias et al., 2009). Indeed, the JRE behaves as a retainer of DIC during the dry season (Cavalcante et al., 2021). The evaporation of seawater increases the salinity in the remaining seawater (evaporation-concentrated seawater) in the same proportion to the amount of seawater that is evaporated, considering that precipitation and evaporation create no salt flux across the air-sea interface. Thus, it is possible to calculate the effect of evaporation in the concentrations of TA and DIC in the estuary. The values of ΔTA and ΔDIC were mostly positive (above the evaporation line), evidencing that estuary is producing TA and DIC in significant amounts by oxic and sub-oxic processes (Figure 5). The estimated net gains of TA and DIC are about 80.8 and 179.2 mmolC m^–2 d^–1, respectively. This gain of DIC is higher than the average of CO₂ emissions in the world’s estuaries (45.2 mmolC m^–2 d^–1; Chen et al., 2013) and higher than the NCP calculated for this estuary (section “pCO2 Variability Amplified by Eutrophication” of this manuscript).

The reactions of reduction and oxidation are coupled to proton production and consumption leading to changes in TA and DIC, particularly in estuaries and coastal regions enriched in organic matter (Abril and Frankignoulle, 2001; Hu and Cai, 2011). The stoichiometry of diagenetic reactions alters TA and DIC in specific ways (Rassmann et al., 2020). The highest deviations of TA and DIC concentrations from the expected evaporation line (Figures 5A,B) occurred in the samples with highest salinities from JRE, indicating that there are net gains of TA and DIC from middle/upper estuarine zones. The data from the JRE are between the slopes of sulfate reduction (vector v), denitrification (vector vi), and aerobic respiration (vector vii). Indeed, the biogeochemical processes of aerobic respiration, ammonification, denitrification, and sulfate reduction explained the gains of TA and DIC in other tropical mangrove-dominated estuaries (Dutta et al., 2019; Akhand et al., 2021). A previous study has suggested the influence of sulfate reduction in the JRE (Lacerda et al., 2013), which is particularly pronounced during the dry period due to the longer water residence time and greater influence of the waters trapped inside the mangrove area and when bidirectional flux can be observed at the middle of the estuary (Dias et al., 2016). Sulfate reduction is the main digenetic pathway of organic matter mineralization in mangroves (Borges et al., 2003; Bouillon et al., 2007). Furthermore, high levels of ammonium (NH₄⁺) and low levels of dissolved oxygen were reported at mangrove tidal creeks in the JRE, indicating processes of ammonification and denitrification (Eschrique et al., 2014; Marins et al., 2020). Increasing levels of NH₄⁺ and PO₄^3– in the middle estuarine portion suggest influences of wastes from shrimp aquaculture ponds and domestic sewage (Eschrique et al., 2014). The net heterotrophy explains the higher production of DIC compared to TA.

Controls of Carbonate Chemistry at Semi-Diurnal and Diel Variabilities

In classical river-dominated estuaries, there is a clear gradient of salinity in the main estuarine channel, with salinities decreasing from the sea to the higher estuary (freshwater domain) (Frankignoulle et al., 1998; Chen et al., 2013). Tropical rivers present lower TA concentrations than temperate/boreal rivers overall (Cai et al., 2008, 2013; Cotovicz et al., 2020a). Consequently, tropical river-dominated estuaries show lower values of pH, DIC, CO₃^2–, Ω_cal, and Ω_ara in the freshwater domain (Salisbury et al., 2008; Hu and Cai, 2011), and these carbonate chemistry parameters tend to increase seaward or during flooding tides. However, the JRE did not present this classical salinity gradient frequently, and principally during dry conditions; to the contrary, the estuary presents increasing salinity, pH, DIC, TA, CO₃^2–, Ω_cal, and Ω_ara during ebb tides (Figure 3). The high residence time of water in the estuary contributes to the high rates of evaporation, resulting in hypersalinity and creating this atypical pattern of carbonate chemistry at semi-diurnal (tidal-driven) time scale. Evaporation would increase solute concentrations (TA, DIC) and salinity simultaneously without changing the TA/salinity and DIC/salinity ratios (Hu et al., 2015; McCutcheon et al., 2019). Therefore, the residual water (after evaporation) shows higher values of DIC and TA concentrations than the initial seawater entering the estuary. Since evaporation also exerts control on major cations (like Ca²⁺) and consequently changes the solubility of calcium carbonate minerals, it may have the potential to change the CO₃^2–, Ω_cal, and Ω_ara (Millero, 2007; McCutcheon et al., 2019). Indeed, changing only salinity and keeping the other parameters constant, the CO₃^2–, Ω_cal, and Ω_ara will be lower in waters with salinity of 40 than in waters with salinity of 35, due to the changes in DIC speciation (Middelburg et al., 2020). This means that the proportion of dissolved CO₂ to the DIC pool increases, whereas the proportion of CO₃^2– decreases with increasing salinity (maintaining TA and DIC constant).

In addition to this atypical semi-diurnal pattern, the carbonate chemistry in JRE also exhibited a clear light-driven diel cycle due to the influence of biological metabolism. Comparing consecutive periods of lowest tidal height occurring at daytime and night-time, there is a marked difference of CO₃^2– concentrations, and values of pH, Ω_cal, and Ω_ara. The values of these parameters are always lower at night-time because of respiratory processes exceed photosynthetic production of organic matter. There is a variety of organic matter sources available in the estuary, with autochthonous (estuarine phytoplankton) and allochthonous (predominantly mangrove detritus) origins (Mounier et al., 2018; Cavalcante et al., 2021). The respiratory processes lead to a production of CO₂, reducing the CO₃^2– concentrations, and values of pH, Ω_cal, and Ω_ara. On the other hand, the concentrations of CO₃^2– and values of pH, Ω_cal, and Ω_ara were higher at daytime because of primary production, with uptake of CO₂ (Cotovicz et al., 2018; Cyronak et al., 2018). The phytoplankton biomass is highest at daytime and during ebb tide conditions (Chl a concentrations reach up to 2.5 μg L^–1), suggesting that the phytoplankton biomass is concentrated in middle/upper estuarine sections.

The DIC/TA ratio gives information about the buffering capacity of seawater; the waters have lowest buffering capacity when DIC/TA ratio are close to 1 (Egleston et al., 2010). Indeed, the CO₃^2–, pH, Ω_cal, and Ω_ara were lowest when the DIC/TA ratio was at highest levels approaching to 1 (Figure 7). As sulfate reduction and denitrification produce TA and DIC in a very close proportion (1:1), aerobic respiration can be considered the main driver of reduced buffering capacity in the estuary, increasing the DIC/TA. Indeed, the prevalence of respiratory processes over primary production reduces buffering capacity (Anthony et al., 2011; Cotovicz et al., 2018; Cyronak et al., 2018). Enhanced respiratory processes caused by coastal eutrophication has recently been described in estuaries, contributing to coastal acidification (Cai et al., 2011; Wallace et al., 2014). In the JRE, the buffering capacity was lower at night-time than at daytime. This diel pattern was verified in other productive coastal ecosystems, following the diel balance between photosynthesis and respiration. The concentrations of CO₃^2–, and values of pH, Ω_cal, and Ω_ara in the JRE were, in average, slightly lower than those found in adjacent shelf waters (Cotovicz et al., 2020b) suggesting a slight acidification driven by eutrophication in the estuary compared to shelf waters. Calcite minerals were observed in JRE sediments at 30 cm depth by X-ray microanalysis (SEM/EDS), suggesting that the trophic state that affects pH, Ωcal, and Ωara has changed through time (Miguens et al., 2010, 2011).

pCO₂ Variability Amplified by Eutrophication

Overall, estuarine waters are oversaturated in CO₂ with respect to atmospheric equilibrium because estuaries are net heterotrophic ecosystems (Borges and Abril, 2011). The JRE follows this pattern, with pCO₂ values always above the atmospheric pCO₂ (∼ 410 uatm). The variability of pCO₂ shows strong influences from semi-diurnal and diel time scales, as described for other carbonate chemistry parameters. The values of pCO₂ increases during ebb tides, and decreases during flood tides, a common pattern found in diverse estuaries worldwide (Borges and Abril, 2011). However, in the JRE the higher pCO₂ values are coincident with higher salinities, which is not commonly described in estuaries. The hypersalinity in the JRE was developed in more confined waters, coinciding the highest organic matter concentrations and nutrients, which comes from diverse sources, both natural and anthropogenic (Lacerda et al., 2013; Eschrique et al., 2014; Mounier et al., 2018; Marins et al., 2020; Cavalcante et al., 2021). These sources of organic matter contribute to maintain high rates of microbial respiration, sustaining a net heterotrophic metabolism. Indeed, there is a strong positive correlation between pCO₂ values and PO₄^3– concentrations, evidencing that microbial production of CO₂ and remineralization of nutrients by respiratory processes are coupled in the estuary. Other probable explanation is that part of the CO₂ and part of the PO₄^3– in the estuary are coming from same allochthonous sources (shrimp farm; domestic effluents) (Marins et al., 2011, 2020), as well in other eutrophic tropical estuaries (Cotovicz et al., 2015, 2018).

The highest levels of pCO₂ occurred at night-time due to the organic matter respiration and absence of photosynthesis (light-driven diel cycle) (Figures 4A,B). Most of the calculated NCP showed negative values (Table 2) indicating that on average the ecosystem is net heterotrophic, confirming an overall metabolism found in estuaries (Borges and Abril, 2011). Even taking into account the net heterotrophy, the estuarine phytoplankton consumes some CO₂ from the water column considering that pCO₂ values were always lower at daytime compared to night-time considering similar tidal heights. The uptake of CO₂ by estuarine phytoplankton is particularly exacerbated in shallow, stratified, and eutrophic ecosystems (Cotovicz et al., 2015). Reduction in pCO₂ over daylight corresponds to a diel increase in water temperature and solar radiation (Reiman and Xu, 2019). The net heterotrophy in the JRE seems to amplify the diel variability of pCO₂. The diel amplitude of pCO₂ values in the JRE was between 219 and 231 μatm. It is highly above than diel variations found in offshore and oligotrophic sites (Dai et al., 2009) but similar to productive coastal bays, nearshore macrophyte meadows, and hypersaline mangrove-dominated wetlands (Borges et al., 2003; Dai et al., 2009; Delille et al., 2009; Saderne et al., 2013), and lower than those found in coral reef ecosystems and mangrove creeks (Dai et al., 2009; Cotovicz et al., 2020a).

In addition to NCP, the contribution of DIC (and pCO₂) from other sources can be relevant in the main estuarine channel, particularly from mangrove ecosystems and effluent discharges. Carbon export from mangrove soils to the ocean occurs due to soil’s semi-diurnal inundation by tide, a process called “tidal pumping” (Becherer et al., 2016; Santos et al., 2021). Waters flowing from mangroves during ebb tide are enriched in DIC and pCO₂ compared to flood tides (Borges et al., 2003). This outwelling of mangrove-derived carbon can partially sustain the high pCO₂ values found in the main estuarine channel of the JRE, as the estuary has about 13,000 ha covered by mangrove forests (Godoy et al., 2018). Furthermore, there is increasing evidence of direct inputs of effluents from shrimp aquaculture ponds and urban discharges into the estuary (Eschrique et al., 2014; Marins et al., 2020; Lacerda et al., 2021). High pCO₂ values associate with the mineralization of domestic organic matter was described in urbanized estuaries (Frankignoulle et al., 1998; Borges and Abril, 2011). The relationship between E-DIC and AOU shows that JRE present values above the 1:1 line (Figure 8), which represents the quotient between photosynthesis and respiration (Borges and Abril, 2011). The deviation above this line confirms additional sources of sources of DIC (and pCO₂) other than NCP. Moreover, there is thermodynamic influences supporting high pCO₂ values. The CO₂ solubility in hypersaline and warm waters is low (Millero, 2007). Thermodynamic calculations using CO2Calc (Robbins et al., 2011) show that the increase of 1°C in water temperature and 1 unit in water salinity increase the pCO₂ in 25 and 16 μatm, respectively. Hypersaline and warm water holds less dissolved CO₂ increasing the pCO₂, a process described in another evaporative estuary (Yao and Hu, 2017; McCutcheon et al., 2019). However, at the diel time scale, the changes of pCO₂ driven by thermodynamics are minor compared to influences driven by ecosystem metabolism and allochthonous sources.

Emissions of CO₂ (Air-Water CO₂ Fluxes) and Implications Considering the Undergoing Processes of Climate Change and Eutrophication Worldwide

Estuaries show large range of variability in the values of gas transfer velocity (k₆₆₀) as a function of wind speed compared to lakes, streams, rivers, and open ocean environments (Borges and Abril, 2011). The gas transfer velocity is highly site-specific and one of the major uncertainties in estimating flux rates for CO₂ (Borges and Abril, 2011; Wanninkhof, 2014; Jeffrey et al., 2018). For that reason, when the gas transfer velocity is not directly measured, it is recommended to quantify the air-water CO₂ fluxes using different parameterizations to provide a wider range of air-water CO₂ flux estimates (Jeffrey et al., 2018; Akhand et al., 2021). That is the reason we have applied four commonly used k₆₆₀ parameterizations (McGillis et al., 2001; Raymond and Cole, 2001; Jiang et al., 2008; Wanninkhof, 2014). The FCO₂ in the JRE showed a strong diel influence driven by wind velocity and pCO₂ disequilibria. Indeed, some estuaries show daytime decrease in pCO₂ in coincidence with increase in average wind speed (Reiman and Xu, 2019). In others, the increased evasion during the day is explained by usually higher wind velocities during this period (Cotovicz et al., 2015; Maher et al., 2015). Summarizing, the highest CO₂ effluxes in the JRE are verified when the gradient of pCO₂ at the air-water interface and wind velocities are greatest.

A global estuarine compilation of air-water CO₂ emissions showed that lower estuaries (with salinities higher than 25) are weak sources of CO₂ (23.0 ± 38.3 mmol C m^–2 d^–1) (Chen et al., 2013). With respect to latitude, lower latitude estuaries exhibited moderate emissions (23.5–0°S: 44 ± 29 mmol C m^–2 d^–1; 0–23.5°N: 39 ± 55 mmol C m^–2 d^–1) (Chen et al., 2013). The average of CO₂ emissions in the middle JRE estuarine portion is 51.9 ± 26.7 mmol C m^–2 d^–1 and higher than values verified in other estuaries with similar characteristics. Our values are close to those found in tropical waters surrounding mangrove forests, with air-water CO₂ fluxes of about 50 mmolC m^–2 d^–1 (Borges et al., 2003; Maher et al., 2015). Overall, the hypersaline estuarine waters hold less dissolved CO₂ contributing to enhanced CO₂ efflux, a pattern found in estuaries located at the northwestern Gulf of Mexico (Yao and Hu, 2017; Yao et al., 2020). Our averaged values of CO₂ outgassing exceeded those found in hypersaline waters in the northwestern Gulf of Mexico estuaries, suggesting that eutrophication and allochthonous contributions from mangrove and effluent discharges in the JRE can amplify the CO₂ efflux in these hypersaline waters.

Figure 9 shows a conceptual model with interrelations between semi-arid climate, river damming, eutrophication, and carbonate chemistry (with emphasis on CO₂ emissions) in the JRE. The presence of organic matter retained by the negative estuary enhances the respiratory processes, decreasing the pH values and DIC/TA ratio (buffering capacity), but increasing the pCO₂ and CO₂ emissions. Consequently, CO₂ emissions by the estuarine waters are about 10 times higher compared to nearshore waters (Carvalho et al., 2017; Cotovicz et al., 2020b). Observations and modeling presume the increase of evaporation compared to precipitation in arid/semi-arid regions as the result of climate change (Chou et al., 2009; Huang et al., 2017) and in JRE this process is potentiated by river damming (Dias et al., 2009; Lacerda et al., 2020). Currently, semi-arid regions cover ∼15% of the Earth’s continental surface (Safriel and Adeel, 2005), and the global area of drylands is estimated to expand ∼10% by 2,100 (Feng and Fu, 2013). Indeed, Godoy et al. (2018) reported a 5.6 mm yr^–1 decrease in annual rainfall over the Jaguaribe River basin from the late 1960s that has accelerated in the past 50 years. The consequences of this decrease in rainfall were potentialized by the construction of five large dams. Concomitantly, the eutrophication of estuaries is a widespread environmental problem and particularly enhanced in densely populated regions (Cloern et al., 2014). Despite the recent efforts to implement wastewater treatment plans in developing countries (Tong et al., 2020), the delivery of nitrogen and phosphorus is continuing to grow with urbanization and the associated population increase (Larsen et al., 2016). The water quality of the Jaguaribe River reflects this tropical scenario of environmental degradation, sustaining an increasing process of eutrophication (Eschrique et al., 2014; Marins et al., 2020). The current levels of Ω_ara and pCO₂ are still higher than 2 (supersaturation) and less than 1,000 μatm (hypercapnia), respectively, and similar to other tropical evaporative estuaries (McCutcheon et al., 2019). However, eutrophication and warming tend to increase the values of pCO₂ and CO₂ efflux in the estuary, which could be critical in the coming years considering the actual scenario of climate change and environmental degradation. The enhanced respiratory processes driven by eutrophication tend to decrease the pH, concentrations of CO₃^2–, and vales of Ω_ara and Ω_cal. The critical levels of Ω_cal and Ω_ara (<1) will probably appear early in the JRE compared to adjacent coastal waters due to the ongoing process of eutrophication. Indeed, the coastal eutrophication could increase the susceptibility of coastal waters to ocean acidification (Cai et al., 2011). The diel variability of pH exceeded 0.1 in the JRE, which is the magnitude comparable to the change in the mean ocean pH during the industrial era (Orr et al., 2005). Finally, the amplitude of diel variations is expected to increase with increasing aquatic pCO₂ values in the future, because of simultaneously decreasing buffer capacity driven by global (ocean acidification) (Schulz and Riebesell, 2013), and local (enhanced respiration) processes (Cai et al., 2011).

Conclusion

The diel variability of carbonate chemistry was investigated in the JRE under severe drought conditions. The estuary shows higher salinity, TA and DIC concentrations compared to adjacent coastal waters, confirming an evaporative and inverse estuary. The salinity, TA and DIC concentrations increase during ebb tides, a different trend compared to most estuaries. However, the concentrations of TA and DIC are above than those calculated by the evaporation model, evidencing additional sources in the estuary (biological metabolism, mangroves, effluent discharges). Overall, the carbonate chemistry parameters were governed by a combination of semi-diurnal and diel tendencies. Semi-diurnal variability was driven by the tidal movement, with increasing pCO₂, and decreasing pH, CO₃^2–, Ω_ara, and Ω_cal during ebbing tides, and an inverse pattern verified during flooding tides. The diel variability was governed by the light-driven biological cycle. Comparing pCO₂ values occurring in similar conditions of tidal heights but at different hour of the day, the highest values of pCO₂ were always verified during night-time due to respiratory processes and absence of photosynthesis. The aquatic pCO₂ values were always above the atmospheric values. The estuary behaved as a CO₂ source, averaging 51.9 ± 26.7 mmol C m^–2 d^–1, which is higher than verified in other similar environments due to enhanced heterotrophy driven by eutrophication. Considering that human population is increasing in coastal regions, and many coastal regions of the world are experiencing increases in aridity as the result of climate change (IPCC, 2021), the present study may provide important insights of future conditions regarding carbon cycling and carbon budget in other aquatic coastal ecosystems worldwide.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

LC: data analysis, writing, reviewing, and editing. RM: conceptualization, experiment design, sampling, funding acquisition, writing, and reviewing. AS: experiment design, sampling, and analytical procedures. All authors contributed to the article and approved the submitted version.

Funding

This study was financed by the Fundação Cearense de Apoio ao Desenvolvimento Científico e Tecnologico (FUNCAP; Proc. No. INT-00159-00009.01.00/19 and FUNCAP PV2-00125-00405.01.00/21), and the Programa de Apoio a Núcleos de Excelência (PRONEX; Proc. No. PR2-0101-0052.01.00/2015).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

LC thanks the Fundação Cearense de Apoio ao Desenvolvimento Científico e Tecnológico (FUNCAP) for a visiting professor grant at the Marine Sciences Institute (LABOMAR)—UFC. Thanks are also due to Francisco ATF Silva, from INPE- Instituto Nacional de Pesquisa Espacial, Eusébio, Ceará, for supervising the equipment developed to these studies in the Brazilian Equatorial Coast. The symbols used in Figure 9 are a courtesy of the Integration and Application Network, University of Maryland Center for Environmental Science (ian.umces.edu/symbols/).

References