Inter-annual Variability of the Carbonate System in the Hypoxic Upper Pearl River Estuary in Winter

Hypoxia has become a universal environmental and ecological problem in recent decades. The Pearl River estuary (PRE), the largest estuary in Southern China, is hypoxic year-round in the upper estuary. This study reports the inter-annual variation between 2005 and 2019 in the carbonate system of the hypoxic upper PRE in winter. In January 2005, both dissolved inorganic carbon (DIC) and total alkalinity (TA) concentrations were >3000 μmol kg–1 at the upstream-most station and decreased sharply downstream. However, DIC and TA were lower, with concentrations of 2300 and 1950 μmol kg–1, respectively, at the upstream-most in January 2019. At salinities >15, both DIC and TA were conservative and reached steady values at the downstream seawater end-member. The upstream-most station was taken as an example to quantify the influences of biogeochemical processes on DIC and TA, including CO2 degassing, organic carbon oxidation, pelagic nitrification, CaCO3 dissolution and benthic release. Among the biogeochemical process, a decrease in CaCO3 dissolution (from 734.4 μmol kg–1 in 2005 to 168.9 μmol kg–1 in 2019) was the major factor driving the decreases of DIC and TA in 2019. In the context of global change, inter-annual variability in biogeochemical process should receive more attention.


INTRODUCTION
Hypoxia in estuaries and coasts is a global environmental and ecological problem (Caballero-Alfonso et al., 2015;Fennel and Testa, 2019). Well-known hypoxic zones include the northern Gulf of Mexico (Rabalais and Turner, 2006), the Baltic Sea (Conley et al., 2011;Bendtsen and Hansen, 2013;Carstensen et al., 2014), Chesapeake Bay (Testa and Kemp, 2012;Li et al., 2016), and the East China Sea off the Changjiang estuary (Li et al., 2002), as well as others. The Pearl River estuary (PRE) is the largest estuary in southern China, and is surrounded by several rapidly developing cities such as Guangzhou, Dongguan, Shenzhen, and others. Hypoxia in the PRE occurs mainly in two zones; one is the year-round hypoxic zone in the upper estuary from Guangzhou to Humen (Zhai et al., 2014), and the other is the seasonal hypoxic zone in summer at the estuary mouth and in waters near Hong Kong (Qian et al., 2018;Wang et al., 2018;Cui et al., 2019). In the hypoxic upper estuary, highly over-saturated CO 2 was observed all year round, with CO 2 partial pressures (pCO 2 ) up to >7000 µatm and dissolved inorganic carbon (DIC) concentration up to >3000 µmol kg −1 in winter (Xu et al., 2005;Dai et al., 2006Dai et al., , 2008Guo et al., 2008;Zhai et al., 2014).
Viewing the entire PRE, Guo et al. (2008) reported the seasonal variability of carbonate parameters in the main channel of the estuary based on field observations from 2000 to 2005; Liang et al. (2020) simulated the spatial distributions of carbonate parameters in the summer and winter of 2006. In the hypoxic upper PRE, Guo et al. (2008) found that both DIC and total alkalinity (TA) showed large seasonal variations, with much higher values in winter (>3000 µmol kg −1 ) than in summer (<1700 µmol kg −1 ). They suggest that seasonal changes in the location of freshwater end-member, biogeochemical processes and hydrologic conditions are the major causes of the variability . More recently, field observations since 2010 have shown lower DIC and TA values (<2500 µmol kg −1 ) in winter in the upper PRE (Dai et al., unpublished data). However, DIC and TA of the river waters were much lower than 2000 µmol kg −1 (Guo et al., unpublished data). This indicates that the very high DIC and TA concentrations previously found in the upper PRE in winters from 2000 to 2005 might be due to strong local additions from biogeochemical processes.
This study aims to: (1) confirm that DIC and TA in the upper PRE decreased after 2005, and (2) determine the key processes dominating the inter-annual variability in these parameters. In order to solve these two questions, we closely compare data collected from January 2005 and January 2019 in the estuary. The 2005 data were previously published , and the 2019 data were collected during a cruise at nearly the same locations in the upper PRE as the 2005 data. During the cruise, in addition to measuring carbonate system parameters, dissolved oxygen (DO) and nutrient measurements, and oxygen consumption and nitrification incubations were conducted using the same methods as in January 2005.

Study Area
The Pearl River is the 2nd largest river in China and the 13th largest river in the world in terms of freshwater discharge, with an annual freshwater discharge of 3.3 × 10 11 m 3 year −1 (Zhao, 1990;Dai et al., 2014). The Pearl River has three major tributaries: the West River, the North River and the East River, and many local rivers and reticulated channels in the delta (Figure 1). The West River originates from the Yunnan-Guizhou Plateau and flows though Yunnan Province, Guizhou Province, the Guangxi Zhuang Autonomous Region and Guangdong Province. The North River and the East River originate from Jiangxi Province and flow through Guangdong Province. The three major tributaries contribute ∼95% of the freshwater discharge of the Pearl River system (Zhao, 1990). The drainage basins of the West River and North River have an abundance of carbonate minerals, while the drainage basin of the East River is rich in silicate minerals. Therefore, the bicarbonate and calcium concentrations in the West and North Rivers are higher than in the East River (Chen and He, 1999). The drainage basin of the Pearl River system is located in the south Asia tropical and subtropical areas, and ∼80% of the freshwater discharge occurs during the wet season from April to October, driven by the Asian monsoon (Zhao, 1990;Dai et al., 2014). The highest discharges of the three main tributaries occur from June to July, and lowest discharges are in winter (December-February, Figure 2).
All runoff from the Pearl River system discharges into the northern South China Sea (NSCS) via eight major outlets (from east to west: the Humen, Jiaomen, Hongqimen, Hengmen, Modaomen, Jitimen, Hutiaomen, and Yamen outlets) and through three sub-estuaries: the Lingdingyang, Modaomen, and Huangmaohai (Figure 1). Lingdingyang is traditionally called the PRE, which collects the freshwater from the eastern four outlets, including the East River and branches of the North and West Rivers.
The lower Pearl River drainage basin has several rapidly developing cities, such as Guangzhou, Dongguan, and others. Rapid economic and population growth have caused serious environmental problems, such as water pollution and hypoxia. In the early 2000s, whole water column hypoxia and very high ammonia concentrations (>800 µmol kg −1 in winter and >300 µmol kg −1 in summer) were observed in the upper estuary from Guangzhou to Humen . Strong aerobic respiration and nitrification induced very high pCO 2 and CO 2 degassing rates (Dai et al., 2006Guo et al., 2008). Due to dilution by seawater and relatively weaker biogeochemical processes downstream, nutrient concentrations and pCO 2 decreased sharply downstream of the Humen Outlet, i.e., in the Lingdingyang (Dai et al., 2006;Guo et al., 2009).
Over the past few decades, sewage discharge and treatment rates have both significantly increased. For example, in the city Guangzhou, urban domestic sewage discharge increased from 7.4 billion tons year −1 in 1995 to 14.3 billion tons year −1 in 2015 (Figure 3A). At the same time, the sewage treatment rate also increased, from <10% in 1995 to 48% in 2005, and then to >85% after 2010 ( Figure 3B). Subsequently, direct sewage discharge without treatment decreased sharply after 2005 ( Figure 3C) (Bulletin of Statistics of Guangzhou's Environmental Status). 1 As a result, the water quality of the channels flowing through Guangzhou and Dongguan has improved considerably. During the January 2005 cruise, the water was black, had a noticeably bad smell, and wildlife such as fish and birds were not observed. However, in 2019, we observed greener water, flock of egrets, and the absence of previous odors. Therefore, 2005 is a  representative year of the huge impacts of direct sewage discharge and severe pollution, while 2019 is a representative year of water quality restoration.

Water Sample Collection
Cruises were conducted in the PRE from January 19-23, 2005, and from January 12-18, 2019, onboard the Yue-Dongguan-Yu 00589 and Yue-Zhu-Yu 31008. At each station, depth profiles of salinity and temperature were recorded with a YSI 6600 multi-parameter monitoring system (2005) or an AML BASE.X (sensor exchangeable Instrument) Conductivity-Temperature-Depth/pressure (CTD) sensor package (2019). Water samples were collected with 2 L or 5 L Niskin bottles.
Sub-samples for DO, pH, DIC/TA were taken with Tygon R tubing free of air bubbles, with ample sample overflow in order to minimize any contamination from atmospheric O 2 or CO 2 . Samples for DO were sampled with 60 mL BOD (biological oxygen demand) bottles and fixed with Winkler reagents (Carpenter, 1965). Samples for pH were taken into FIGURE 3 | Urban sewage discharges, sewage treatment rate, and sewage discharges without treatment from the city Guangzhou. Each column represents one year. Data were from the Bulletin of Statistics of Guangzhou's Environmental Status (1995, http://sthjj.gz.gov.cn/gkmlpt/content. 120 mL amber glass bottles and poisoned with 100 µL of a saturated HgCl 2 solution. During the January 2005 cruise, DIC samples were sampled into 40 mL amber glass vials with 50 µL saturated HgCl 2 added; TA samples were sampled into 120 mL high-density Polyethylene (HDPE) bottles with 100 µL saturated HgCl 2 added. During the 2019 cruise, samples for DIC and TA measurements were taken into 250 mL borosilicate bottles with grounded stoppers and poisoned with 200 µL of the saturated HgCl 2 solution. All DIC and TA samples were stored in the dark at room temperature until analysis. Samples for nutrient and calcium (Ca 2+ ) measurements were filtered with 0.45 µm polyacetate filters. Samples for nitrite (NO 2 − ), nitrate and nitrite (N + N), and soluble reactive phosphorus (SRP) were frozen and stored at −20 • C until analysis. Samples for ammonium (NH 4 + ) and silicate (Si) measurements were stored at −4 • C until analysis. Samples for Ca 2+ measurements were stored in the dark at room temperature.
CO 2 partial pressures of surface water and the atmosphere was measured with an underway system integrating an air-water equilibrator and a Li-Cor 7000 or Picarro G2101-i CO 2 analyzer. Six air-based CO 2 standard gases ranging 200 to 10000 parts per million provided by the National Reference Material Research Center of China were used to calibrate the measurements. Data processing was the same as previously described in Guo et al. (2009). Wind speed was measured with an R.M. Young model 05106 marine wind monitor at ∼10 m above the sea surface. The accuracy of wind speed measurements was ±0.3 m s −1 .
In February 18-19 of 2019, surface water samples were collected at gauge stations in the lower North and East Rivers (Shijiao in the North River and Boluo in the East River) to measured DIC, TA, nutrients, and calcium. The samples were taken with a 2 L plexiglass sampler. The sub-sampling, treatment and storage methods are the same as those during the January 2019 cruise.

Rate Incubation Experiments
Bulk oxygen consumption incubations were conducted at the upstream-most station (Figure 1) following the method of Dai et al. (2006). Briefly, surface water was pumped into a clean 400 L plastic box. After oxygenation, the water was transferred into two 20 L HDPE cubitainers (Thermo Fisher Scientific) immediately. 10 mL of a saturated HgCl 2 solution was added to one of the cubitainers as a control. The incubation was conducted in the dark and the temperature was controlled to match that of the surface water by using running water circulated outside of the cubitainers. Samples were taken every 3-6 h to measure DO. Bulk oxygen consumption rates were estimated from the change in the DO concentration over time.
Incubation experiments for estimating nitrification rates were conducted at the same stations as the bulk oxygen consumption incubations following the method of Dai et al. (2008). Briefly, surface water was sampled with 5 L Niskin bottles. Triplicate samples were filled into 4 L narrow-neck amber glass bottles. Allylthiourea (ATU, final concentration of 100 mg L −1 ) and NaClO 3 (final concentration of 10 mg L −1 ) were added to the bottles to inhibit the oxidation of NH 4 + and NO 2 − , respectively. The third bottle was used as a control without reagent addition. Water samples were taken every 4-8 h to measure NO 2 − . NH 4 + and NO 2 − oxidation rates were estimated from the change in the concentration of NO 2 − over time.

Sample Analysis
Dissolved oxygen samples were measured onboard with the classic Winkler method, and calibrated using a potassium iodate standard solution provided by the Reference Material Center of the Marine Ecology and Environmental Laboratory of the 2nd Institute of Oceanography (Ministry of Natural Resources of the People's Republic of China). The precision of the measurements was better than ±2 µmol kg −1 . pH samples were measured onboard within 2 h of sampling with a Corning 350 or Orion 3 Star pH meter and an Orion 8102BN Ross combination electrode calibrated against three NBS buffers (4.01, 7.00, and 10.01 at 25 • C, provided by Thermo Fisher Scientific). Samples for pH measurements and buffers were placed in a constant temperature bath at 25 ± 0.01 • C for about 1 h before their pH values were measured. Therefore, the measured pH values were at 25 • C (pH 25 ) and given on the NBS scale. DIC and TA samples were measured onboard within 1 day of sampling during the 2005 cruise and within 1 week after the 2019 cruise. Analyses of DIC and TA followed the methods previously described by Cai et al. (2004). DIC was measured by collecting and quantifying the CO 2 released from the sample upon acidification with a non-dispersive infrared detector  using an Apollo SciTech model AS-C3 DIC Analyzer with a precision of ±2 µmol kg −1 . TA was determined by Gran titration with hydrochloric acid using an automated alkalinity titrator (Apollo SciTech model AS-ALK1+) with a precision of better than ±2 µmol kg −1 . Both DIC and TA measurements were calibrated with certified reference material (CRM) provided by Dr. A. G. Dickson at the Scripps Institution of Oceanography, to obtain an accuracy of better than ±2 µmol kg −1 . For the January 2019 cruise, Ca 2+ was measured using the ethylene glycol tetraacetic acid (EGTA) titration method developed by Lebel and Poisson (1976) and modified by Cao and Dai (2011), with a precision of ±0.06%. The measurements were calibrated with artificial seawater based CaCl 2 solution. Salinity standard seawater of the International Association for the Physical Sciences of the Oceans (IAPSO) (Batch P158) was adopted as a reference. The measured Ca/salinity ratio of the IAPSO standard seawater was 291.6, which was slightly higher than the reported value of Batch P79 and P86 (290.5 and 290.9) (Olson and Chen, 1982), but slightly lower than the reported values of 292.0 for Batch P78 (Olson and Chen, 1982), 293.0 for Batch P67 (Kanamori and Ikegami, 1980), and 292.3 for Batch P147 (Cao and Dai, 2011). The deviations of the Ca/salinity ratios might be attributed to differences in the ratios of the different batches of the IAPSO standard seawater, and/or the system error among the different laboratories or operators. Nevertheless, our Ca/salinity ratio was very close to the average of the above literature values (291.7). The precision of our measurements (±0.06%) were ±5.8 µmol kg −1 for seawater with salinity of 33. For the January 2005 cruise, Ca 2+ was measured with the ethylenediaminetetraacetic acid (EDTA) titrimetric method (China-EPA, 1987), with a precision of ±0.41%, (±39.5 µmol kg −1 for seawater with salinity 33).
Ammonium samples were measured onboard within 3 h of sampling using the indophenol blue method, with a detection limit of 0.1 µmol L −1 (Pai et al., 2001). NO 2 − was determined using the pink azo dye method; N + N was determined by coppercadmium column reduction and the pink azo dye photometric method; SRP was determined with the phosphomolybdenum blue photometric method; silicate samples were determined using the silicomolybdic blue photometric method (Hansen and Koroleff, 1999). NO 2 − , N + N, SRP and silicate samples were all measured with an AA3 Auto-Analyzer (Bran-Lube, GmbH). The detection limits for NO 2 − , N + N, SRP and silicate were 0.04, 0.10, 0.08, and 0.16 µmol L −1 , respectively, and the analytical precision was better than 1% for N + N and silicate and 2% for SRP (Han et al., 2012).

Data Processing
The saturation state index ( ) was defined as the product of the concentrations of Ca 2+ ([Ca 2+ ]) and carbonate ([CO 3 2− ]) divided by the apparent solubility product of CaCO 3 (K sp * ) (Eq. 1).
[Ca 2+ ] was the measured Ca 2+ concentration. [CO 3 2− ] was calculated from DIC and TA with the program CO2SYS (Version 14) (Pierrot et al., 2006). The dissolution constants of carbonic acid are from Millero (2010). The CO 2 solubility coefficient is from Weiss (1974) and the sulfate dissociation constant from Dickson (1990). The PO 4 3− and SiO 2 data are the measured cruise data. K sp * at 1 standard atmosphere pressure was calculated from temperature and salinity according to the formula of Mucci (1983). The pressure correction of the K sp * of calcite was from Ingle (1975); the pressure correction of the K sp * of aragonite was from Millero (1979). The calculated was 0.02-0.37 higher than the values calculated with CO2SYS, which might suggest the contribution of excess Ca 2+ produced by the dissolution of CaCO 3 in the upper estuary.
DO saturation degree (DO%) was defined as the ratio of the measured DO to the DO at saturation (DO sat ) (Eq. 2). DO sat was calculated according to the empirical formula of Benson and Krause (1984).
( 2) The net CO 2 flux (F CO2 ) between surface water and the atmosphere (or air-sea CO 2 flux) was calculated using the following formula: where s is the solubility of CO 2 (Weiss, 1974), pCO 2 is the pCO 2 difference between the surface water and the atmosphere, and k is the CO 2 transfer velocity. k was parameterized using the empirical function of Sweeney et al. (2007). U 10 is the wind speed at 10 m above sea level (in m s −1 ). Here, the daily average wind speeds were adopted in the FCO 2 calculations.

Salinity, Temperature and DO Saturation Degree
The freshwater sources of the upper PRE are mainly the North River and the East River. The long-term (2000-2018) monthly average freshwater discharge rates of the North and East Rivers are ∼500 and ∼420 m 3 s −1 in winter, equivalent to 25-40% of the freshwater discharge in summer (Figure 2). The low freshwater discharge in winter may have resulted in the higher salinities observed in the upper PRE in winter than in summer (salinity was 0 from Guangzhou to Humen in summer, Guo et al., 2008). During both the January 2005 and 2019 cruises, the salinity of the upstream-most station at Guangzhou was <0.5, and increased downstream (Figures 4A,B). Salinity at the Humen Outlet ranged from 10 to 20, and increased to >25 in Outer Lingdingyang and 30-33 at the estuary mouth.
Surface water temperatures ranged from 14.6 to 18.8 • C in January 2005, with the higher temperatures found in the low latitude zone (<22 • N). In January 2019, temperatures ranged from 15.9 to 20.5 • C, with higher temperatures in the Lingdingyang and outside of the estuary mouth than in the upper estuary (Figures 4C,D).
The DO saturation degree in the Guangzhou region ranged from ∼4-31% (12-93 µmol kg −1 ) in January 2005, and most of this region was hypoxic. Downstream from Guangzhou, the DO saturation degree increased to 60-80% at the Humen Outlet and 80-100% in the Lingdingyang. Outside of the Lingdingyang, DO was slightly oversaturated (Figure 4E). In January 2019, DO saturation in the Guangzhou region ranged from 15 to 65% (45-200 µmol kg −1 ), and hypoxic zones were observed at the two east-west ends of this region. Downstream of Guangzhou, the degree and pattern of DO saturation was similar to that observed in January 2005 ( Figure 4F). Comparing the data of the two January cruises reveals that DO was slightly higher in 2019 than in 2005 in the upper estuary (Figures 4E,F).

Spatial Distribution of Carbonate System Parameters, DO and Nitrogen-Containing Nutrients
In January 2005, DIC and TA values were as high as ∼3300 and 3100 µmol kg −1 , respectively, at the upstream-most station, and decreased sharply downstream. At salinities >15, both DIC and TA showed either a linear decrease or increase with salinity. However, the distributions of DIC and TA in January 2019 were different. At the upstream-most station, DIC and TA concentrations were 2300 and 1950 µmol kg −1 , respectively, ∼1000 µmol kg −1 lower than those measured at this station in 2005. They both decreased with salinity and reached their minima at salinity of ∼5. At salinities >5, both DIC and TA increased with salinity. Similar to the 2005 survey, at salinities >15, both DIC and TA were conservative (Figures 5A-1,B-1). As DIC and TA concentrations in the low salinity zone were much lower in 2019 than in 2005, their concentrations in almost the entire estuarine mixing zone were lower in 2019 than in  Guo et al. (2008), and nutrient data are from Dai et al. (2008). Legends are the same for all panels. Figures 5A-2,B-2). At salinities >30, TA values in January of 2005 and 2019 were consistent, which indicated the stable characteristics of the seawater compared to the estuarine water ( Figure 5B-1).

(
In January 2005, the pH at the upstream-most station was ∼7.2. pH decreased downstream and reached a minimum (6.9) at a salinity of ∼5 (Figure 5C-1), which was due mainly to the high rates of nitrification and aerobic respiration . At salinities >5, pH increased downstream and reached ∼8.11 at a salinity of ∼33 (Figure 5C-1). In January 2019, the distribution of pH in the upper estuary was slightly different from that in 2005. Although the pH at the upstream-most station was as low as 7.06, the pH in the middle of the Guangzhou region was as high as ∼7.4 (Figure 5C-2). pH decreased again to 7.1 further downstream (60 km upstream Humen Outlet), and then increased with salinity to 8.05 at a salinity of 32.5. Compared with 2005, DO in the Guangzhou region in 2019 was generally higher. However, at the same site downstream of Guangzhou (from the distance of −50 km downstream), the DO saturation degree in 2019 resembled January 2005 (Figure 5D-2). The relatively higher pH at stations in the Guangzhou region in January 2019 was accompanied by higher DO concentrations (Figures 5C-2,D-2).
In January 2005, NH 4 + concentrations in the most upper estuary were as high as 800 µmol kg −1 and decreased sharply downstream to 3 µmol kg −1 at salinities >30 (Figures 5E-1,E-2). N + N concentrations increased from 177 µmol kg −1 at the upstream-most station to a maximum of 372 µmol kg −1 at a salinity of ∼5, and then decreased to 10 µmol kg −1 at salinities >30 (Figures 5F-1,F-2). The sharp decrease in NH 4 + concentrations and maximum N + N concentrations in the low salinity zone were due to the presence of strong nitrification Guo et al., 2008). In January 2019, the NH 4 + concentration at the upstream-most station was 150 µmol kg −1 , much lower than in 2005. Consequently, NH 4 + concentrations across the entire estuarine mixing zone were much lower in 2019 than in 2005. However, the N + N concentration at upstreammost stations (up to 463 µmol kg −1 ) was higher than in 2005 (Figures 5E-1,E-2,F-1,F-2). This might be due mainly to the nitrification process during sewage treatment (EPA, 1993). The total concentration of NH 4 + and N + N at the upstream-most station was also much lower in 2019 than in 2005 (1009.6 vs. 622.7 µmol kg −1 ), which might be due to the denitrification process during sewage treatment (EPA, 1993).

Distribution of Calcium, and the Saturation State Index of Calcium Carbonate
Calcium concentrations ranged from 1798 to 9410 µmol kg −1 in January 2005 and from 989 to 9529 µmol kg −1 in January 2019. Although generally the Ca 2+ concentration increased with salinity, the overall distribution of Ca 2+ was slightly different between the two winter cruises. In 2005, Ca 2+ concentrations in the low salinity zone were higher than in 2019. If the conservative mixing line between the North River and seawater was taken as a reference, Ca 2+ concentrations showed additions in the low salinity zone in January of both 2005 and 2019, and the addition in 2005 was higher than in 2019.
The saturation state indices of both aragonite and calcite were <1 in the upper PRE (upstream of Humen Outlet, salinity <15), suggesting under-saturation for both aragonite and calcite minerals. In the area downstream of the Humen Outlet, Ca increased to ∼5.2 in 2005 and 4.0 in 2019, and Ar increased to 3.2 in 2005 and 2.4 in 2019, at the seawater end-member (Figure 6). The higher saturation state indices of calcite and aragonite in the high salinity zone in 2005 were due mainly to a local phytoplankton bloom, which was consistent with the relatively higher pH and DO, but lower DIC, in 2005.

Bulk Oxygen Consumption Rate
In January of 2005 and 2019, DO decreased almost linearly during incubations without HgCl 2 addition, while DO concentrations remained almost constant for the samples with HgCl 2 added (control samples, Figure 7). The bulk oxygen consumption rates were 55.88 and 49.44 µmol L −1 d −1 in 2005 and 2019, respectively ( Table 1).
In addition, we may use in situ DO concentrations to estimate the bulk oxygen consumption rate. In January 2005, DO saturation at the upstream-most station was 311.6 µmol kg −1 , while the in situ DO concentration was 14.9 µmol kg −1 . Therefore, the DO consumption was 296.7 µmol kg −1 . Assuming a residence time of 5 days , the DO consumption rate would be 59.3 µmol L −1 d −1 . This value is consistent with the value of 55.9 µmol L −1 d −1 estimated based on the incubations (Table 1). In January 2019, the in situ DO concentration at Guangzhou was 45.6 µmol kg −1 and the DO saturation was 302.7 µmol kg −1 . Therefore, the in situ DO consumption was 257.1 µmol kg −1 . If a residence time of 5 days is adopted , the bulk oxygen consumption rate is estimated as 51.4 µmol L −1 d −1 , which is reasonably consistent with values based on the incubations (49.4 µmol L −1 d −1 , Table 1).
Compared to other estuaries, the bulk oxygen consumption rate measured in the upper PRE is similar to that estimated in the polluted Seine River estuary in the spring and fall of 1996 [20-68 µmol L −1 d −1 assuming an average depth of 5 m, Garnier et al. (2001)]. However, it is higher than that found in the inner Scheldt estuary in the winter of 2003 [18-21 µmol L −1 d −1 , Gazeau et al. (2005)], or the polluted Huangpu River flowing through the city Shanghai (a branch of the lower Changjiang) in the winter and fall of 2005, which was 4-11 µmol L −1 d −1 (Zhai et al., 2007). However, it is lower than in salt marsh waters of the southeastern United States measured during the fall of 1995 and summer of 1996 [80 µmol L −1 d −1 , Cai et al. (1999)].

Nitrification Rate
In January 2005 and 2019, NO 2 − concentrations decreased in the samples with ATU added (inhibiting the oxidation of NH 4 + to NO 2 − ) during the incubations, suggesting the conversion of NO 2 − to NO 3 − . However, the NO 2 − concentration increased in the incubations with NaNO 3 added (inhibiting the oxidation  (Table 1). Nitrification rates in January of 2019 were higher than in January of 2005, but lower than in spring and summer [up to 31.5 µmol L −1 d −1 , Dai et al. (2008)]. Although NH 4 + concentrations were lower in 2019 compared to 2005, the nitrification rate increased, which suggests that NH 4 + concentrations are controlled by multiple factors. As the influencing factors of nitrification are beyond the scope of this study, we will not discuss it further here.

Processes Dominating DIC and TA Addition/Removal in the Upper PRE
Compared to January 2005, DIC and TA concentrations at the upstream-most station in January 2019 were both ∼1000 µmol kg −1 lower, but relatively consistent in the seawater end-member between both years. Influenced by the water characterized by very high DIC and TA at the upstream-most station, DIC and TA in the estuarine mixing zone at salinities <20 in January 2019 were much lower than those in January 2005 (Figure 5). It should be noted that DIC was slightly lower and pH and DO were slightly higher at the seawater end-member in 2005 (Figures 5A-1,C-1,D-1), which might be due to a weak local phytoplankton bloom.
The DIC and TA minima at salinity of ∼5 in 2019 were observed where the East River converges into the main channel of the estuary. Therefore, they might be due to the influence of East River water which was characterized by lower DIC and TA (the blue triangles in Figures 5A-1,B-1). DIC and TA in the East River were 605.3 and 566.6 µmol kg −1 in February 2019. The characteristically low DIC and TA concentrations of East River water were also reported in Guo et al. (2008). However, DIC and TA did not show this minima in January 2005, which might be due to large DIC and TA additions from local biogeochemical processes, similar to the case in the Guangzhou region.
Although both DIC and TA decreased downstream from the Guangzhou region, the very high DIC and TA values at the upstream-most station (3329.2 and 3093.5 µmol kg −1 in January 2005 and 2301.3 and 1947.2 µmol kg −1 in January   Biogeochemical processes in the Guangzhou region of the upper PRE were very strong (Dai et al., 2006. Processes influencing DIC and TA concentrations include aerobic respiration, nitrification, air-water exchange of CO 2 , and others . Benthic release is also an important DIC and TA source for the water column (Cai et al., 2015). Additionally, calcium carbonate dissolution or precipitation might also impact the carbonate system in estuarine environments (Abril et al., 2004;Macreadie et al., 2017;Su et al., 2020). This complicated mixing scheme in the upper PRE makes it difficult to quantify the DIC and TA budget at all stations; thus, we take the upstream-most station as an example to quantify the influences of biogeochemical processes. At this station, the freshwater source was mainly the North River.

Air-Water CO 2 Exchange
The upper PRE was a strong CO 2 source. CO 2 evasion decreases DIC in the water, but has no influence on TA. The average air-water CO 2 exchange rate in the Guangzhou region was 175.5 ± 14.5 and 174.5 ± 1.5 mmol m −2 d −1 in January of 2005 and 2019, respectively ( Table 1). If a residence time of 5 days and average water depth of 5 m are adopted (Zhao, 1990;Guo et al., 2008), the influence of CO 2 evasion would have decreased DIC by 175.5 and 174.5 µmol kg −1 in January 2005 and 2019, respectively. This shows the influence of CO 2 evasion on DIC was nearly identical in 2019 and 2005.

Pelagic Nitrification
The stoichiometric relationships of NH 4 + and NO 2 − to DO and H + are expressed in Eqs. 5 and 6 (Dai et al., 2006). Only the oxidation of NH 4 + to NO 2 − influences TA, and the ratio of change in TA to NH 4 + concentrations is 2. Nitrification decreases TA in the water, but has no influence on DIC. At the upstreammost station, NH 4 + oxidation rates were 3.86 µmol L −1 d −1 in January 2005 and 6.53 µmol L −1 d −1 in January 2019 (Table 1). Therefore, the rates of change in TA were 7.7 and 13.1 µmol L −1 d −1 in 2005 and 2019, respectively. Similarly, taking a residence time of 5 days , the TA change would be −38.

Organic Carbon Oxidation
The total oxygen consumption rate includes the oxygen consumption due to nitrification and organic carbon oxidation. During nitrification, the stoichiometric ratio of O 2 to NH 4 + is 1.5 during the NH 4 + oxidation process, and the ratio of O 2 to NO 2 − is 0.5 during the NO 2 − oxidation process (Eqs. 5 and 6). The NH 4 + oxidation rate was 3.86 µmol L −1 d −1 in 2005 and 6.53 µmol L −1 d −1 in 2019, and the NO 2 − oxidation rate was 4.62 µmol L −1 d −1 in 2005 and 7.75 µmol L −1 d −1 in 2019 (Table 1). Therefore, the oxygen consumption rate induced by nitrification was 8.10 µmol O 2 L −1 d −1 in January 2005 and 13.69 µmol O 2 L −1 d −1 in January 2019, respectively.
As the total oxygen consumption rate at this station was 55.88 µmol L −1 O 2 d −1 in January 2005 and 49.44 µmol L −1 O 2 d −1 in January 2019, the DO consumption rate induced by the oxidation of organic carbon (excluding oxidation of nitrogen) was 47.78 µmol L −1 O 2 d −1 in January 2005 and 35.75 µmol L −1 O 2 d −1 in January 2019. According to the stoichiometric ratio of organic carbon oxidation (excluding nitrification of NH 4 + , Eq. 7), the DIC production rate due to organic carbon oxidation would be 47.78 µmol C L −1 d −1 in January 2005 and 35.75 µmol C L −1 d −1 in January 2019. If a residence time of 5 days is assumed , the DIC addition due to organic carbon oxidation would have been 238.9 µmol kg −1 in January 2005 and 178.8 µmol kg −1 in January 2019.
(CH 2 O) 106 + 106 O 2 → 106 CO 2 + 106 H 2 O (7) CaCO 3 Dissolution As the saturation state of CaCO 3 in the upper PRE ( < 0.5) was well below the saturation level ( = 1), CaCO 3 should dissolve. CaCO 3 dissolution influences both DIC and TA. We estimated the Ca 2+ addition resulting from CaCO 3 dissolution according to the difference between the observed and estimated conservative mixing Ca 2+ concentrations. To estimate the conservative Ca 2+ concentration, we assume two end-member mixing between North River water and seawater. The relationship of the conservative mixing line of Ca 2+ with salinity is expressed by Eq. 8.
The relatively lower CaCO 3 dissolution in January 2019 might be due to the relatively higher pH (Figures 5C-1,C-2). Low pH is favorable for CaCO 3 dissolution. Therefore, the relatively lower pH in 2005 might enhance CaCO 3 dissolution. However, CaCO 3 dissolution increases pH (Eq. 9), which adds the complexity of the relationship between pH and CaCO 3 dissolution.

Benthic Release
Denitrification can also influence DIC and TA (Chen, 2002), but as the water column was generally oxygenated (although DO < 60 µmol kg −1 ) we assumed denitrification was negligible. However, both nitrification and denitrification might occur in the sediments of the PRE, which also influences the carbonate system in the water (Xu et al., 2005). Other anaerobic reactions (reduction of Fe/Mn and/or sulfate) in the sediment might also affect the carbonate system in the water. Anaerobic reactions in sediment consume protons and increase TA in the sediment porewater. Combined with the organic carbon degradation, DIC in porewater in the PRE was much higher than in the water column (Cai et al., 2015). The interactions between sediment and water at the water-sediment interface cause the sediment to release DIC and TA to the water column. A similar phenomenon of sedimentary reactions altering the carbonate system in the water column was also observed in the northern Gulf of Mexico (Hu et al., 2017;Berelson et al., 2019).
As we didn't measure benthic fluxes during the cruises, we used the benthic DIC flux measured in November 2013 and the benthic TA/DIC ratio reported in the literature to estimate the influence of benthic release on the water column DIC and TA. In November 2013, the benthic DIC release rate at the upstream station was 1200 ± 300 mmol m 2 d −1 (Cai et al., 2015).
To account for the large influence of temperature on the rates of biogeochemical processes (Rabus et al., 2002;Govorushko, 2012), we calibrated the biogeochemical process rates in the sediment according to Eq. 10.
Q 10 refers to the multiple of the increase in the rate of biochemical reactions per temperature increase of 10 • C; k refers to the slope of the linear regression between temperature and the logarithm of rates.
The reported Q 10 of anaerobic reactions in sediment ranges from 2 to 3 (Pomeroy and Wiebe, 2001;Kirchman et al., 2009), so a Q 10 of 2.5 was taken in the estimation. The temperature in November 2013 when the sampling of Cai et al. (2015) was conducted was 22 • C, while temperatures during our January 2005 and 2019 samplings were 16.0 and 16.4 • C. Lowering the temperature by 6 and 5.6 • C decreases the reaction rates by 1.73 and 1.67, respectively. Assuming that other conditions during our cruises were the same as November 2013, the benthic DIC release rate was 692.5 µmol m 2 d −1 in January 2005 and 718.3 µmol m 2 d −1 in January 2019.
The ratio of the benthic TA/DIC flux is 0.57-0.95 in the hypoxic northern Gulf of Mexico (Berelson et al., 2019). If this ratio range is similar in the PRE, the benthic TA flux would be 394.7-657.9 mmol m −2 d −1 in January 2005 and 409.5-682.4 mmol m −2 d −1 in January 2019.
The Guangzhou region of the upper PRE was well-mixed vertically, so the released benthic DIC and TA were likely homogenized throughout the whole water column. Taking an average water depth of 5 m and a residence time of 5 days (Zhao, 1990;Guo et al., 2008), the DIC and TA additions from the benthic flux would be 695.2 and 394.7-657.9 µmol kg −1 in January 2005, and 718.3 and 409.5-682.4 µmol kg −1 in January 2019. Therefore, considering the uncertainties in the benthic DIC and TA flux estimates, no conspicuous changes were observed between January 2005 and January 2019.

Inter-annual Variations in DIC and TA in the Upper PRE
The quantification of the above biogeochemical processes at the upstream-most station is shown in Table 2  The only DIC sink for upper PRE water in January of both 2005 and 2019 was CO 2 evasion to the atmosphere, while the DIC sources were organic carbon oxidation, CaCO 3 dissolution and benthic release. The DIC removal due to CO 2 evasion to the atmosphere was almost the same during the two winter cruises (175.5 and 174.5 µmol kg −1 in 2005 and 2019, respectively). Although the DIC addition from organic carbon oxidation in January 2019 was slightly lower than in January 2005 (238.9 µmol kg −1 in 2005 vs. 178.8 µmol kg −1 in 2019), DIC removal due to CO 2 evasion was almost offset by the DIC addition resulting from organic carbon oxidation during both cruises relative to the size of the overall DIC budget. The most important sources of DIC were CaCO 3 dissolution and benthic release. As the benthic release rates were estimated based on previously reported benthic DIC fluxes and only the temperature effect was considered, estimates of benthic DIC release were similar during the two winter cruises as the temperatures were similar. A comparison between the two winter cruises indicates the largest differences were in the DIC addition due to CaCO 3 dissolution; DIC addition due to this process was 565.5 µmol kg −1 lower in January 2019 than in January 2005, controlling the much lower overall DIC addition observed in January 2019.
The only major sink for TA of the water was pelagic nitrification. As the nitrification rate in January 2019 was slightly higher than in January 2005, the TA removal in 2019 (65.3 µmol kg −1 ) was slightly greater than in 2005 (38.6 µmol kg −1 ). As was the case with DIC, the benthic TA release during both the 2005 and 2019 cruises was similar. However, the TA addition due to CaCO 3 dissolution in January 2019 was 1131.0 µmol kg −1 lower than in 2005, and dominated the decrease in the local TA addition in 2019 in the upper PRE.
Although the above estimates generally explain the interannual variations in DIC and TA in the upper PRE between 2005 and 2019, there are differences between the sum of the estimated and the observed addition values. These differences might be due to the different dynamic conditions during the different cruises as well as the uncertainties in estimating benthic DIC and/or TA fluxes. The uncertainties in water residence time, average depth, and other parameters might also add to the uncertainties of the estimates. In addition, the influence of complicated mixing on the biogeochemical processes across the entire estuarine mixing zone is yet to be quantified. The relationship between decreasing CaCO 3 dissolution and the improvement of water quality also needs further examination. Ecological environmental improvement is a universal goal in China and across the world, and this may in turn influence the cycles of biogenic elements. Similar inter-annual variability in biogeochemical parameters or mass fluxes may manifest in other estuaries and coastal zones, which could also influence ocean element cycling and even local climate change. This phenomenon deserves more attention.

CONCLUSION
Dissolved inorganic carbon and TA distributions in the PRE were determined in January of 2005 and 2019. Although they show little difference in the high-salinity lower estuary, both DIC and TA in the heavily perturbed upper estuary showed a ∼1000 µmol kg −1 decrease in 2019 compared to 2005. The upstream-most station was used as an example to quantify the influences of biogeochemical processes on DIC and TA. For DIC, the CO 2 degassing was almost compensated by organic carbon oxidation, while benthic release and CaCO 3 dissolution contributed to the major additions. For TA, benthic release and CaCO 3 dissolution were the major additions, and pelagic nitrification contributed to slight removal. Although the above biogeochemical processes differed between the 2 years, the decrease in CaCO 3 dissolution dominated the DIC and TA decreases in the upper PRE in 2019 compared to 2005. In the context of global change, inter-annual variability in the biogeochemical parameters of estuaries and coasts might be universal and deserves more attention.

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

AUTHOR CONTRIBUTIONS
XG collected the oxygen consumption rate data from the January 2005 cruise. All authors collected the data from the January 2019 cruise except LW. XG designed the cruise, took the DIC, TA and pH samples. XS and YG conducted the nitrification and oxygen consumption incubation experiments. YL collected the nutrient and calcium samples. YX measured DO and pCO 2 . TH measured the ammonium samples. YG and YL measured the N + N, phosphate and silicate samples. XS and YL measured the DIC, TA and calcium samples. YG and YL collected the North and East Rivers data in February 2019. XG drafted the manuscript and all authors participated the discussion on data interpretation.