Impacts of anthropogenic inputs on hypoxia and oxygen dynamics in the Pearl River estuary

In summer, the Pearl River estuary (PRE) experiences hypoxia, largely driven by the high input of fresh water with low dissolved oxygen (DO), abundant nutrients and particulate organic carbon from the Pearl River network. In this study, we used a well-validated physical–biogeochemical model together with a DO species-tracing method to study the responses of hypoxia and oxygen dynamics to the anthropogenic perturbations of riverine inputs (i.e. DO, nutrients, and particulate organic carbon) in July–August 2006. Model results showed that hypoxia in the PRE was most sensitive to riverine inputs of particulate organic carbon, followed by DO concentrations and nutrients. Specifically, a 50 % decrease (increase) in riverine input of particulate organic carbon led to a 47 % decrease (64 % increase) in hypoxic area, with the sediment oxygen demand and water column production being the two most important processes contributing to changes in DO concentration. Changes in the riverine inputs of DO and nutrients had little impact on the simulated hypoxia because of the buffering effects of re-aeration (DO fluxes across the air–sea interface); i.e. the re-aeration responded to the changes in surface apparent oxygen utilization (AOU) associated with river-induced variations of oxygen source and sink processes. The PRE features shallow waters (with averaged depth of 10 m) in which oxygen provided by the re-aeration could penetrate to bottom waters via vertical diffusion and largely offset the changes in DO contributed by other oxygen source and sink processes. This study highlights the importance of re-aeration in reducing hypoxia variability in shallow estuaries.

1 Introduction

Recent decades have seen a decline in dissolved oxygen (DO) concentrations in most of the coastal oceans because of intensifying anthropogenic disturbances, leading to an increase in the occurrence and intensity of hypoxic conditions (Diaz and Rosenberg, 2008). Relations between the riverine nutrient loading and the hypoxic conditions (DO < 2 mg L⁻¹) have been well documented in many coastal hypoxic systems, such as the Changjiang Estuary (Li et al., 2011; Ning et al., 2011), the Chesapeake Bay (Du and Shen, 2015; Hagy et al., 2004), and the northern Gulf of Mexico (NGOM) (Forrest et al., 2011; Justić et al., 2003). The classic paradigm for explaining the relations is that excessive nutrient inputs to the coastal oceans stimulate the high primary productivity there, and the subsequent decomposition of the organic matter in the bottom water consumes a significant amount of DO that leads to hypoxia. As a result, nutrient reduction has been proposed to alleviate hypoxia in many hypoxic systems (e.g. the Chesapeake Bay – Scavia et al., 2006, and the NGOM – Justić et al., 2003). Recent years have also seen an increasing number of studies showing that climate variation contributes to the spreading hypoxia in coastal oceans. The climate variation can change the ocean circulation or the vertical stratification to alter the balance between the oxygen source and sink processes (Rabalais et al., 2010). A modelling study conducted in the Chesapeake Bay has shown high correlations between the climate variation, stratification, and the observed DO (Du and Shen, 2015). In addition, global warming, as a symptom of climate variation, is another factor that can enhance hypoxia. For example, Laurent et al. (2018) predicted prolonged and more severe hypoxia in the northern Gulf of Mexico under a projected future (2100) climate state in which the global warming leads to a reduction in oxygen solubility and increased stratification.

Figure 1 (a) A bathymetric map showing the Pearl River network and the Pearl River estuary, (b) computational cross sections for the 1-D river network model, and (c) the model grid for 3-D estuary model. Red numbers in (a) represent islands which are not marked on the map: 1 is Qi'ao Island, 2 is Hengqin island, 3 is Gaolan Island, and 4 is Inner Lingding Island.

The Pearl River estuary (PRE) is located on the Pearl River delta (Fig. 1a) and has a drainage area of 452 000 km². Previous studies have reported some summer hypoxic events in the PRE and explored the underlying mechanisms. Yin et al. (2004) suggest that stratification and estuarine circulations are two primary processes controlling hypoxia in the PRE. Rabouille et al. (2008) compare the hypoxic conditions among four hypoxic systems and demonstrate the significance of tidal mixing to break hypoxia in the PRE. Zhang and Li (2010) further suggest that the contributions of biogeochemical processes to hypoxia in the PRE are also important. By conducting the oxygen balance analysis, they show that sediment oxygen demand (SOD) is the dominant sink for oxygen. A more recent study by Wang et al. (2017) further point out that the balance of oxygen in the PRE is mainly controlled by the source and sink processes occurring in local and adjacent waters, among which the re-aeration (DO fluxes across the air–sea interface) and SOD determine the spatial distributions and durations of hypoxia in the PRE.

As a distinct river-dominated estuary, the PRE receives 3.3×10¹¹ m³ yr⁻¹ of fresh water (Ou et al., 2009; Zhang and Li, 2010) along with a large amount of nutrients from the Pearl River network (Fig. 1a), i.e. 5.6×10⁵ t yr⁻¹ of dissolved inorganic nitrogen (DIN) and 9.9×10³ t yr⁻¹ of dissolved inorganic phosphorus (DIP) (Hu and Li, 2009). Both dissolved inorganic nitrogen and phosphorus loadings have increased by about 60 % from 1970 to 2000 and are predicted to increase by 2 times in 2050 due to the fast-growing agriculture and urbanization (Strokal et al., 2015). Understanding the response of hypoxia and oxygen dynamics to changes in nutrient loading in the PRE is hence valuable for hypoxia prediction and management.

In addition to the nutrient loading, the particulate organic carbon (POC) is another important form of anthropogenic input that influences the hypoxia in the estuary ($\sim \mathrm{2.5}\times {\mathrm{10}}^{\mathrm{6}}$ t yr⁻¹ from the Pearl River network – Zhang et al., 2013). The POC can fuel the SOD when deposited and mineralized in the sediment layers, which has been found to dominate the DO depletions within the bottom waters of the PRE (Yin et al., 2004; Zhang and Li, 2010). In coastal systems POC is often derived from the dead phytoplankton (Green et al., 2006), while in the PRE the POC mainly originates from the riverine inputs (Ye et al., 2017; Yu et al., 2010). This suggests the importance of studying the impact of riverine POC on hypoxia in the PRE.

In some cases, hypoxia may also be induced by advection of low-oxygen waters (Grantham et al., 2004; Montes et al., 2014; Wang, 2009; Wang et al., 2012). For example, Wang (2009) demonstrates that the hypoxia development in the Changjiang estuary is largely due to the Taiwan Warm Current bringing low-oxygen waters to the hypoxic zone. As to the PRE, the impact of riverine input of low-oxygen waters on hypoxia is also worth investigating considering the large amount of river discharge entering the estuary and that there has been hypoxia observed in its upper reaches (He et al., 2014).

Collectively the previous studies show that both natural and anthropogenic processes greatly contribute to hypoxia in the PRE. Understanding the respective roles of these two types of processes is important for faithfully predicting future hypoxic events under enhanced human activities and climate variations, which is useful for designing effective management strategies to prevent or remediate the hypoxic conditions in the PRE. Here we focus on the role of human activities, i.e. different anthropogenic inputs, on hypoxia and oxygen dynamics in the PRE, whereas the role of natural processes will be reported in our future work. Specifically, we explore the impact of varying anthropogenic inputs (riverine nutrients, POC, and DO) on hypoxia and oxygen dynamics in the PRE by using a three-dimensional (3-D) coupled physical–biogeochemical model. The DO species-tracing method introduced in Wang et al. (2017) is applied to isolate the effects of each oxygen source and sink process and to elucidate their interactions in this shallow and river-dominated estuarine system.

2 Method

2.1 Model description and validation

2.1.1 Model description

Physical model

Our physical model is a 1-D–3-D coupled model which incorporates the Pearl River network and the PRE (see Fig. 1b, c for locations) into a single framework to resolve the dynamic interactions between these two regions (Hu and Li, 2009). This coupled model was first developed with the biological and sediment models for the Pearl River estuary system to study the water, nutrients, and sediment flux budgets between the river network and estuary (Hu and Li, 2009; Hu et al., 2011). Thereafter, it was extended to study the hypoxia (Wang et al., 2017) and the nutrient fluxes across the water–sediment interface (Liu et al., 2016) in the PRE.

The coupled model uses an explicit coupling approach to incorporate the 1-D model for the Pearl River network and the 3-D model for the PRE through the eight outlets (including Humen, Yamen, Hongqili, Hengmen, Modaomen, Jitimen, Hutiaomen, and Yamen; see Fig. 1a for their locations). At each time step, the 3-D model is forced by the simulated river discharges from the 1-D model, and as a feedback, sends its simulated water levels through eight outlets to the 1-D model as the downstream boundary conditions for the next time step. More detailed descriptions of the model methodology can be referred to in Hu and Li (2009).

The cross-sectional integrated 1-D model solves the Saint Venant equations of mass and momentum conservation by using a Preissmann implicit scheme and an iterative approach in the well-mixed river network. Figure 1b shows that the Pearl River network is discretized into 1726 computational cross sections, 189 nodes (interactions between the different river branches), 5 upper boundaries (i.e. Shizui, Gaoyao, Shijiao, Laoyagang, and Boluo), and 8 lower boundaries (the eight outlets). The upper boundaries of the 1-D model are specified by the real-time observations of river discharges or water levels. The lower boundaries use simulated water levels from the 3-D model. Initial conditions are set to zero for water levels and velocities, and model time step is 5 s.

The 3-D model is based on the Estuaries and Coastal Ocean Model with Sediment Module (ECOMSED; HydroQual Inc., 2002), which has been extensively used to study the hydrodynamics in estuaries. The model has 183×186 horizontal grid cells with a resolution ranging from 400 m inside the Lingdingyang Bay to 4 km near the open boundaries (Fig. 1c), and has 16 terrain-following sigma layers with refined resolution near the surface and bottom layers. The horizontal mixing is parameterized by a Smagorinsky-type formula (Smagorinsky, 1963) and the vertical mixing is calculated by the Mellor–Yamada level 2.5 turbulent closure model (Mellor and Yamada, 1982). The 3-D model is forced by the 6 hourly winds and 3 hourly surface heat fluxes from ERA-Interim (ECMWF reanalysis, http://www.ecmwf.int/en/research/climate-reanalysis/era-interim, last access: 2015). Three open boundaries are specified by a monthly averaged profile of salinity and temperature (Hu and Li, 2009). Tides are introduced at the open boundaries using the water levels from the Oregon State University Tidal Data Inversion Software (OTIS). Freshwater inputs from the Pearl River network to the estuary use the river discharges simulated by the 1-D model.

Figure 2 Conceptual framework for the RCA model with a sediment flux module (Zhang and Li, 2010). DO represents dissolved oxygen, PHYT represents phytoplankton, POC represents particulate organic carbon, DOC represents dissolved organic carbon, NH₄ represents ammonia nitrogen, NO₂₃ represents nitrite and nitrate nitrogen, PON represents particulate organic nitrogen, DON represents dissolved organic nitrogen, DPO₄ represents dissolved inorganic phosphorus, POP represents particulate organic phosphorus, DOP represents dissolved organic phosphorus, DSi represents dissolved silica, BSi represent biogenic silica, and SOD represents sediment oxygen demand.

The physical model is run from 1 November 2005 to 31 December 2006. More detailed descriptions and configurations can be found in Hu et al. (2011) and Wang et al. (2017).

Biogeochemical model

The biogeochemical model is the row–column AESOP model (RCA; HydroQual Inc., 2004) that solves the mass balance equations for 26 state variables involved in five interactive cycles (i.e. the nitrogen cycle, the phosphorus cycle, the carbon cycle, the silicon cycle, and the oxygen dynamics). Interactions between these state variables with atmosphere and sediment are illustrated in Fig. 2.

The equation of DO (mg O₂ L⁻¹) is given by the following:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial \text{DO}}{\partial t}}+u{\displaystyle \frac{\partial \text{DO}}{\partial x}}+v{\displaystyle \frac{\partial \text{DO}}{\partial y}}+w{\displaystyle \frac{\partial \text{DO}}{\partial z}}-{\displaystyle \frac{\partial }{\partial x}}\left(E_x{\displaystyle \frac{\partial \text{DO}}{\partial x}}\right)\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\partial }{\partial y}}\left(E_y{\displaystyle \frac{\partial \text{DO}}{\partial y}}\right)-{\displaystyle \frac{\partial }{\partial z}}\left(E_z{\displaystyle \frac{\partial \text{DO}}{\partial z}}\right)\\ \text{(1)}& {\displaystyle }& {\displaystyle }=\text{WCP}+\text{REA}-\text{SOD},\end{array}$$

where x and y represent the horizontal coordinates and z the vertical coordinate; u, v, and w (m s⁻¹) represent velocity components in x, y, and z coordinates; and E_x, E_y, and E_z (m s⁻²) are dispersion coefficients. The velocity components and dispersion coefficients are computed by the physical model.

The term WCP represents the gross DO production rates in the water column ($\mathrm{mg}{\mathrm{O}}_{\mathrm{2}}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{day}}^{-\mathrm{1}}$), hereafter the water column production, which is the combination of photosynthesis, respiration, nitrification, and oxidation. Detailed equations for each component of the water column production can be seen in Appendix A. According to the DO budget analysis in Wang et al. (2017), the photosynthesis and respiration are two major oxygen source and sink processes in the water column. Considering that photosynthesis and respiration are both closely and directly correlated to the phytoplankton dynamics, they have similar distributions and responses to the external forcing. We therefore use the water column production to represent the net effects of water column on the DO and hypoxia.

The term REA represents the re-aeration ($\mathrm{mg}{\mathrm{O}}_{\mathrm{2}}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{day}}^{-\mathrm{1}}$) at the air–sea interface, given as follows:

$$\begin{array}{}\text{(2)}& \text{REA}=k_{\mathrm{a}}{\mathit{\theta }}_{\mathrm{a}}^{T-\mathrm{20}}\left({\text{DO}}_{\mathrm{sat}}-\text{DO},\right)\end{array}$$

where DO_sat represents the DO concentration at saturation (mg O₂ L⁻¹), which is dependent on salinity and temperature; k_a is the surface mass transfer coefficient (day⁻¹); and θ_a is a temperature coefficient (dimensionless). Values for these parameters can be seen in Table A2.

The term SOD represents the sediment oxygen demand ($\mathrm{mg}{\mathrm{O}}_{\mathrm{2}}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{day}}^{-\mathrm{1}}$) at the water–sediment interface and Δz represents the thickness of the respective bottom grid cell (m).

$$\begin{array}{}\text{(3)}& \text{SOD}={\displaystyle \frac{s\left(\text{DO}-{\text{DO}}_{\mathrm{sed}}\right)}{\mathrm{\Delta }z}}\end{array}$$

where s represents the transfer coefficient between the sediment and overlying water (m day⁻¹), and DO_sed represents DO concentrations (mg O₂ L⁻¹) in the sediment layers. In the RCA, a sediment flux module is incorporated to simulate the depositional flux of particulate organic matter (i.e. particulate organic carbon, particulate organic nitrogen, and particulate organic phosphate), the diagenesis processes in the sediment, and the transport of nutrients and DO from the sediment to the overlying water (Fig. 2). Detailed descriptions about the sediment flux module can be seen in Appendix B.

The simulation period for our biogeochemical model is the same as the physical model. Initial conditions were obtained from a 2-month spin-up simulation, which was repeated 3 times to reach a steady state. River boundary conditions of biogeochemical variables were derived from the monthly observations in 2006 collected by the State Oceanic Administration (including nutrients and DO) and from a previous study (including different classes of dissolved organic carbon, particulate organic carbon, dissolved organic nitrogen, particulate organic nitrogen, dissolved organic phosphorus, and particulate organic phosphorus) (Liu et al., 2016). The open-boundary conditions of biogeochemical variables were specified following Zhang and Li (2010).

Figure 3 (a) Schematic diagram illustrating the mixing process of dissolved oxygen in the estuary and (b) schematic plot for dissolved oxygen versus salinity (the solid black curve line) during mixing in the estuary. C1 represents the concentrations in seawater, while C0 represents the concentrations in river water.

2.1.2 Model validation

The physical–biogeochemical model has been validated against available observations during July 1999 in Hu and Li (2009) and July–August 2006 in Wang et al. (2017). We briefly summarize the validation results in 2006 below.

Being the coupling interface between the 1-D model and the 3-D models, the eight outlets serve as the lower boundaries of the 1-D model and the river boundaries of the 3-D model. It follows that the simulation of eight outlets is of great importance to the robustness of the 1-D–3-D coupled model. Model–data comparisons of water levels were conducted for eight stations, including six outlets (i.e. Jiaomen, Hengmen, Modaomen, Jitimen, Hutiaomen, Yamen) and two other stations (i.e. Zhuhai and Wanshan) in Wang et al. (2017), with the locations of the stations in their Fig. 3. The normalized root-mean-square difference (RMSD) of water levels falls within 0.30 of the standard deviation of the observations and the correlation coefficient between the simulated and observed water levels exceeds 0.95. This indicates that the coupled physical model is able to resolve the interactions between the river network and the estuary well. In addition, the tidal variations and the spring–neap tidal cycles in the PRE are well reproduced.

The PRE is characterized by the large extended river plume in the summer. Therefore, the model-simulated salinity and temperature fields were validated against 146 profiles of salinity and temperature collected by an estuary-wide monitoring cruise. The comparisons show small normalized RMSDs (both <0.60 of standard deviations of observations) and high correlations (>0.90 for salinity and >0.80 for temperature) between the model and observations, indicating that the coupled physical model is robust enough to reproduce the broad-scale features and intra-seasonal patterns of the main hydrodynamic features in the PRE.

For validation of biogeochemical fields, the simulated DO concentrations were validated against 53 oxygen profiles collected at 4 different cruises and distributed estuary-wide. The point-to-point comparisons show that the simulated DO concentrations agree well with observations, with the normalized RMSD below 0.8 standard deviation of the observations and the vast majority (85 %) of the normalized errors falling within 1 standard deviation of the observations. Model–data comparisons of bottom DO concentrations further show that the model is able to reproduce the spatial distribution of the observed bottom DO and hypoxia. We have also assessed model skills in resolving source and sink processes associated with DO concentration. We found that the simulated spatial distributions and magnitudes of the re-aeration, respiration, and the SOD rates are similar to those of previous observational studies (see Table 3 in Wang et al., 2017). The simulated chlorophyll-a, primary productivity and particulate organic carbon, which largely determine the respiration and the SOD rates (Zhang and Li, 2010), are also consistent with historical estimations. This suggests that our model is able to reproduce the oxygen dynamics properly.

In short, the model validation in Wang et al. (2017) indicates that our physical–biogeochemical model is robust enough to simulate the hydrodynamics and biogeochemical cycles in the PRE and is skillful in simulating summer hypoxia in 2006.

2.2 The DO species-tracing method

The DO exhibits non-conservative behaviour during mixing in the estuary because of the oxygen source and sink processes described in Sect. 2.1.1 (Fig. 3a). As shown in Fig. 3b, the DO concentrations are controlled by both the conservative (represented by the theory mixing curve) and the non-conservative effects (represented by the shading areas). The conservative effects are associated with physical advection and diffusion, while the non-conservative effects are due to the oxygen source and sink processes (i.e. re-aeration, the water column production, and the SOD). Quantifying the relative contributions of the respective effect is important for understanding the DO dynamics during the mixing in the estuary. In a 0-D system, the non-conservative effects can be easily estimated as the products of time intervals and rates of corresponding source and sink processes. However, in a river- and tide-dominated estuary such as the PRE, this estimation is not straightforward because of the spatial connections of each source and sink process occurring in different locations. To address this problem, the DO species-tracing method (referred to as the physical modulation method in Wang et al., 2017) was introduced and implemented in our previous study to investigate the mechanisms of hypoxia in the PRE. By dividing the DO into different DO species, the tracing method can track the DO contributed by different source and sink processes. For example, Wang et al. (2017) found that about 28 % of surface DO supplied by the re-aeration penetrated to the bottom waters and hence modulated hypoxia in the PRE. In this study, the DO species-tracing method is used to track contributions of each source and sink process to the DO dynamics and hypoxia under the different riverine input scenarios. Interactions between the oxygen source and sink processes will be investigated as well.

The DO species-tracing method is incorporated into the biogeochemical model by explicitly including four numerical oxygen species as model tracers to track the DO contributed by the lateral boundary conditions (DO_BC (mg O₂ L⁻¹)), air–sea re-aeration (DO_REA (mg O₂ L⁻¹)), water column production (DO_WCP (mg O₂ L⁻¹)), and SOD (DO_SOD (mg O₂ L⁻¹)) (Table A1). Equations of the four numerical oxygen species are given below:

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle \frac{\partial {\text{DO}}_{\mathrm{BC}}}{\partial t}}+\text{tran}\left({\text{DO}}_{\mathrm{BC}}\right)=\mathrm{0}\text{(5)}& {\displaystyle }& {\displaystyle \frac{\partial {\text{DO}}_{\mathrm{REA}}}{\partial t}}+\text{tran}\left({\text{DO}}_{\mathrm{REA}}\right)=\text{REA}\text{(6)}& {\displaystyle }& {\displaystyle \frac{\partial {\text{DO}}_{\mathrm{WCP}}}{\partial t}}+\text{tran}\left({\text{DO}}_{\mathrm{WCP}}\right)=\text{WCP}\text{(7)}& {\displaystyle }& {\displaystyle \frac{\partial {\text{DO}}_{\mathrm{SOD}}}{\partial t}}+\text{tran}\left({\text{DO}}_{\mathrm{SOD}}\right)=\text{SOD},\end{array}$$

and

$$\begin{array}{ll}\text{(8)}& {\displaystyle }& {\displaystyle }\text{DO}={\text{DO}}_{\mathrm{BC}}+{\text{DO}}_{\mathrm{REA}}+{\text{DO}}_{\mathrm{WCP}}-{\text{DO}}_{\mathrm{SOD}}{\displaystyle }& {\displaystyle }\text{tran}\left(\text{DO}\right)=\text{tran}\left({\text{DO}}_{\mathrm{BC}}\right)+\text{tran}\left({\text{DO}}_{\mathrm{REA}}\right)+\text{tran}\left({\text{DO}}_{\mathrm{WCP}}\right)\\ \text{(9)}& {\displaystyle }& {\displaystyle }-\text{tran}\left({\text{DO}}_{\mathrm{SOD}}\right),\end{array}$$

where “tran” represents the physical transport processes, i.e. advection $(u\frac{\partial }{\partial x}+v\frac{\partial }{\partial y}+w\frac{\partial }{\partial z})$ and diffusion $(-\frac{\partial }{\partial x}\left(E_x\frac{\partial }{\partial x}\right)-\frac{\partial }{\partial y}\left(E_y\frac{\partial }{\partial y}\right)-\frac{\partial }{\partial z}\left(E_z\frac{\partial }{\partial z}\right))$. REA ($\mathrm{mg}{\mathrm{O}}_{\mathrm{2}}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{day}}^{-\mathrm{1}}$), WCP ($\mathrm{mg}{\mathrm{O}}_{\mathrm{2}}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{day}}^{-\mathrm{1}}$), and SOD ($\mathrm{mg}{\mathrm{O}}_{\mathrm{2}}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{day}}^{-\mathrm{1}}$) are the re-aeration, water column production, and SOD, which represent the net effects of the air–sea interface, the water column, and the water–sediment interface on the oxygen, respectively. Values of these terms are obtained from the biogeochemical model at each time step.

According to Eq. (4), the DO_BC concentrations are only controlled by advection and diffusion. By assigning the initial conditions and lateral boundary conditions of DO_BC the same as those for DO, the mixing curve of DO_BC will overlap the theory mixing curve shown in Fig. 3b. It follows that the DO_BC represents the conservative effects, while the DO_REA, DO_WCP, and DO_SOD that include the oxygen source or sink term represent the non-conservative effects.

Eqs. (8) and (9) suggest that the DO concentration and its transport flux equal the sum of the concentrations and transport fluxes of the four DO species, the validity of which has been tested and confirmed in Wang et al. (2017). They show that there is little discrepancy between the DO concentrations calculated by Eqs. (9) and (1), with 97 % of the differences within the range −2 % to 6 % of the averaged DO concentrations. The hourly time series of domain-averaged DO calculated by the DO species-tracing method also agree well with that calculated by the biogeochemical model with the R-square coefficient >0.99 and the regression slope close to 1:1. In addition, the horizontal advective fluxes, vertical advective fluxes, and vertical diffusive fluxes calculated by the tracing method are found to agree well with the respective fluxes calculated by the biogeochemical model, indicating that the tracing method is able to satisfactorily reproduce the physical transport processes of DO.

2.3 Model experiments

We conducted three groups of sensitivity experiments to study the response of hypoxia and oxygen dynamics to different scenarios of riverine inputs. Each group has two simulations, where the concentration of one type of the riverine inputs at eight river outlets is decreased or increased by 50 %. These simulations are named Base, RivDO−50%, RivDO+50%, RivNtr−50%, RivNtr+50%, RivPOC−50% and RivPOC+50%, with the basic information of each simulation presented in Table 1.

Table 1 Overview of model experiments.

The Base simulation uses the realistic riverine inputs as mentioned in Sect. 2.1.1. In the Base simulation, the DO concentration in the Humen outlet, the largest river outlet in the PRE, is set to 4 mg L⁻¹ based on observations nearby. The RivDO−50% simulation in which DO concentration from the eight outlets is decreased by 50 % represents the scenario in which hypoxia has developed in the Humen outlet, which has been reported in previous studies (e.g. He et al., 2014). In contrast, the RivDO+50% simulation, in which the DO concentration from the eight outlets is increased by 50 % to be close to that from the open boundaries, represents the scenario in which the riverine input of DO is free from the anthropogenic impact. As for nutrient simulations, the RivNtr+50% and RivNtr−50% simulations increase and decrease nutrient concentrations from all eight outlets by 50 %, respectively. The resulting riverine inputs in the two simulations will be close to the scenarios in 2050 and 1970 as reported by Strokal et al. (2015). Note that, in the nutrient simulations, the concentrations of all nutrients (including dissolved silica, dissolved inorganic phosphorus, ammonia nitrogen, and nitrite and nitrate nitrogen) are set to vary at the same percentage at which the effects of different combinations of changes in nutrients are not considered. The hydrodynamic conditions are identical in all experiments.

The hypoxic extent in different simulations is quantified by the expected hypoxic area and hypoxic volume:

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }\text{Hypoxic area}=\sum p\cdot \mathrm{\Delta }s\text{(11)}& {\displaystyle }& {\displaystyle }\text{Hypoxic volume}=\sum p\cdot \mathrm{\Delta }v,\end{array}$$

where Δs, Δv, and p are the area, the volume, and the hypoxic frequency of each grid cell. The hypoxic frequency p is calculated by the following:

$$\begin{array}{}\text{(12)}& p={\displaystyle \frac{N_{\mathrm{h}}}{N_{\mathrm{s}}}}\cdot \mathrm{100}\mathit{\%},\end{array}$$

where N_h is the number of hours in which hypoxia occurs, and N_s is the total number of hours. In this study, the threshold of hypoxia is defined as 3 mg L⁻¹ (Luo et al., 2008; Rabalais et al., 2010).

Figure 4 (a) Annual cycles of the model-simulated monthly hypoxic area in 2006 of the PRE, (b) annual cycle of the total river discharges in 2006 (blue bars) and during 1999–2010 (error bars represent a standard deviation around the climatological mean values).

Figure 5 Spatial distributions of bottom DO (a–f) and hypoxic frequency (g–l) during May–October. Hypoxia is defined as DO concentration below 3 mg L⁻¹.

Figure 6 Spatial distributions of DO concentrations (a, b, c) and hypoxic frequency (d, e, f) in the bottom layer for DO concentration simulations. The DO concentration is averaged over July and August 2006.

3 Results

3.1 Hypoxia in the Pearl River estuary

As shown in Fig. 4, hypoxia in the PRE starts to develop in April, peaks in August, and disappears in October, which is highly correlated (R²=0.91) with the annual cycle of total river discharges with a time lag of 1 month. Figure 5 shows the model-simulated DO distributions and hypoxic frequency in the bottom layer during the May–October. In May, the hypoxia is confined to the upstream of the Modaomen sub-estuary. In June, the bottom DO declines along the west coast of the PRE and the hypoxia starts to develop near Gaolan Island (see Fig. 1 for its location). The hypoxia extends eastward to near Hengqin Island in July and August. After August, the hypoxia retreats westward and almost disappears in October. Unlike the large spatial extent of hypoxia observed in the Changjiang Estuary (Wang, 2009; Wang et al., 2012) and the NGOM (Rabouille et al., 2008), hypoxia in the PRE is confined to a small area as a result of the SOD and the re-aeration (Wang et al., 2017).

In 2006, oxygen observations are only available in July and August, which have demonstrated the occurrence of hypoxia. No observations are available for validating the model-simulated hypoxia in other months. We have collected and analysed the oxygen observations from 1993 to 2009. However, the available observations are insufficient to resolve the annual cycle of hypoxia in the PRE. To our knowledge, there are currently few studies on the annual cycle of hypoxia in the PRE due to the scarcity of observations. Discussions in this study therefore focus on hypoxia in July–August when the distinct hypoxia was both observed and simulated by the model. Another motivation of focusing on July and August is that these 2 months are among the typical wet seasons in the PRE (Fig. 4b), which is in line with our study on the effects of riverine inputs.

3.2 Response of hypoxia and oxygen dynamics to riverine DO inputs

Figure 6 shows the comparisons of bottom DO concentrations and hypoxic frequency during July–August for different DO simulations. In the RivDO−50% simulation, the spatial distribution of bottom DO is similar to that in the Base simulation except that hypoxia additionally occurs near the river outlets due to the inputs of low-oxygen waters from the upstream river network (Fig. 6b, e). We have also examined the impact of reducing riverine DO in regions farther away from the river outlets by excluding the hypoxic region near the river outlets. In this case the expected hypoxic area in RivDO−50% is only 2 % higher than that in the Base simulation, while the hypoxic volume is 26 % higher (Fig. 7a), indicating that the thickness of hypoxic water is greatly increased in RivDO−50%. In contrast, the RivDO+50% simulation yields higher bottom DO concentrations, leading to reductions in hypoxic area and volume by 23 % and 30 % (Fig. 7a).

Figure 7 The percentage changes in the hypoxic area and hypoxic volume in each simulation relative to the Base simulation (a). The changes in each DO species averaged over the high-hypoxic-frequency zone (denoted as the white contour in Fig. 6) in DO simulations (b), nutrient simulations (c), and POC simulations (d) relative to the Base simulation.

Figure 8 The spatial distribution of DO_BC (a, b, c) and DO_REA (d, e, f) concentrations at the bottom layer for three DO simulations. The white contour represents the high-frequency zone, and the red box represents the west of the lower estuary.

Figure 7b, c, d further show the changes in each DO species averaged over the bottom layer of the high-frequency zone for different simulations relative to the Base simulation. The high-frequency zone here is defined as the area encompassed by the 10 % isoline of July–August averaged hypoxic frequency and is denoted by the white contour in Fig. 8. To provide more insights into the response of different oxygen species to riverine inputs, the spatial distributions of DO_BC and DO_REA in the bottom water are shown in Fig. 8. Differences in DO_WCP and DO_SOD concentrations between simulations are much smaller and hence omitted here.

Halving the riverine DO inputs in the RivDO−50% simulation yields lower DO_BC concentrations but higher DO_REA in the bottom water (Figs. 7b and 8). The decrease in DO_BC concentrations is largely balanced by the increase in DO_REA concentration in RivDO−50% simulation, which ultimately reduces the magnitude of changes in hypoxic extent responding to the reduced riverine DO input. On the contrary, the RivDO+50% simulation leads to higher DO_BC concentrations (Figs. 7b and 8c) but lower DO_REA in the bottom water (Figs. 7b and 8f), which together reduces the net increase in bottom DO (Fig. 7b).

The re-aeration buffering effects can be explained by the surface apparent oxygen utilization (AOU, the difference between the actual DO concentration and its saturation at a known temperature and salinity). As shown in Eq. (2), the re-aeration is a function of surface AOU. Halving the riverine DO inputs decreases the DO concentrations in entire water column and therefore increases the surface AOU, which ultimately results in an increase in re-aeration rate. In our model simulations, the surface domain-averaged saturated DO concentration is ∼7.42 mg L⁻¹, while the surface domain-averaged DO concentration in the Base and RivDO−50% simulations are 6.81 and 6.57 mg L⁻¹. The surface AOU for the RivDO−50% simulation is 39 % higher than that for the Base simulation, which is consistent with the 38 % increase in re-aeration rate for the RivDO−50% simulation.

3.3 Response of hypoxia and oxygen dynamics to riverine nutrient inputs

As shown in Fig. 7a, perturbing riverine nutrient inputs by 50 % has a relatively weak impact on hypoxic extent (changes are within 10 %). Among all the oxygen sink and source processes, the water column production and re-aeration are the two that are most sensitive to variations in nutrient inputs. Halving the nutrient inputs by 50 % in the RivNtr−50% simulation remarkably reduces the primary productivity and water column production rates, which in turn increases the surface AOU that facilitates the re-aeration. The increase in DO_REA in the bottom water via vertical diffusion offsets ∼60 % of the total DO loss associated with the reduced nutrient inputs in the high-hypoxic-frequency zone (Fig. 7c). As a result, the hypoxic area and hypoxia volume only increase by about 10 % in the RivNtr−50% simulation relative to the Base simulation. In contrast, the RivNtr+50% simulation yields higher water column production and a lower re-aeration rate, while the changes in the two balance each other out, and hence the simulation only leads to 4 % and 6 % decreases in hypoxic area and hypoxic volume relative to the Base simulation.

3.4 Response of hypoxia and oxygen dynamics to riverine POC inputs

As shown in Fig. 7, perturbing the riverine inputs of POC by 50 % leads to significant changes in DO concentrations and hypoxic extent. In the RivPOC−50% simulation, the DO concentration increases by 0.56 mg L⁻¹ in the high-hypoxic-frequency zone and the hypoxic area and hypoxic volume decrease by 50 % and 64 %. On the contrary, RivPOC+50% simulation leads to a significant decrease in the DO concentration, causing an extension of hypoxic area by 64 % and a doubling of hypoxic volume (Fig. 7a).

As to oxygen dynamics, the RivPOC−50% simulation leads to a significant decline in the SOD rate (Fig. 7d) and increase in the water column production rate (Fig. 7d) as a result of the lower inputs of POC weakening the light attenuation in PRE. The combination of lower SOD and higher water column production rates increases oxygen concentration by 0.81 mg L⁻¹ in the bottom waters of the high-hypoxic-frequency zone (Fig. 7d). However, decreasing the riverine inputs of POC in the RivPOC−50% simulation simultaneously weakens the re-aeration due to the decreased surface AOU. As a result, nearly 27 % of the increased DO concentrations is offset by the decreased re-aeration in the high-hypoxic-frequency zone. In contrast, increasing the riverine inputs of POC in the RivPOC+50% simulation increases the SOD rates but weakens the water column production rates, which consequently reduces bottom water oxygen; nevertheless, 26 % of the oxygen loss is offset by the enhanced re-aeration in this simulation.

To understand the impact of changing the riverine inputs of POC on the water column production rates, we further examine how phytoplankton growth responds to varying riverine inputs of POC. The equation for the phytoplankton growth rate G_P (day⁻¹) can be written as follows:

$$\begin{array}{}\text{(13)}& G_{\mathrm{P}}=G_{\mathrm{Pmax}}\cdot G\left(T\right)\cdot G\left(I\right)\cdot G\left(N\right),\end{array}$$

where G_Pmax represents the maximum growth rate at the optimum conditions (day⁻¹); G(T), G(I), and G(N) represent the limitations by temperature, light, and nutrients. These limitation coefficients are non-dimensional scale values ranging from 0 to 1, with 0 representing no growth and 1 no limitation. The two POC simulations and the Base simulation have identical physical processes and hence the same temperature limitation. Table 2 shows that changing the riverine inputs of POC has little impact on nutrient limitation but leads to large variations in light limitation, suggesting that the riverine inputs of POC can significantly affect the phytoplankton growth through light shading effects.

Table 2 Comparisons of nutrient limitation and light limitation on the growth of phytoplankton for Base and two POC simulations. Values are averaged over the bottom layer of the PRE. The lower values represent the stronger limitation.

Figure 9 Budget of −DO_SOD for the upper layer, middle layer, and bottom layer in the PRE for the Base simulation. Blue arrows represent sediment oxygen demand, red arrows represent the vertical diffusion, orange arrows represent vertical advection, and green arrows represent horizontal advection. Positive values mean the source effects, while the negative values mean the sink effects of the sediment oxygen demand on DO concentrations (unit: $\mathrm{mg}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{day}}^{-\mathrm{1}}$).

Considering the important role of re-aeration in POC simulations, we further quantify how re-aeration responds to the SOD by conducting a diagnostic analysis of DO_SOD in July and August (Fig. 9). Three vertical layers are defined: the upper layer (top 20 % of the water column), middle layer (middle 60 % of the water column), and bottom layer (20 % of the water column above the sediment). Note that horizontal diffusion is omitted in the diagnostic analysis because its magnitude is much smaller than other terms. Diagnostic analysis of other DO species can be seen in Figs. 11 and 12 of Wang et al. (2017). As shown in Fig. 8, the SOD can affect the DO concentrations in the upper layer indirectly through the interactions with the vertical advection, the vertical diffusion, and the horizontal advection as explained below. First, the SOD consumes bottom DO by 0.53 $\mathrm{mg}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{day}}^{-\mathrm{1}}$ and decreases the upward advective DO fluxes, reaching the upper layer by 0.34 $\mathrm{mg}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{day}}^{-\mathrm{1}}$. Second, the deoxygenation induced by SOD can increase the vertical DO gradient and facilitate the downward vertical diffusion of oxygen by 0.02 $\mathrm{mg}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{day}}^{-\mathrm{1}}$ from the upper layer. Finally, the decreased upper DO concentrations affect the horizontal outfluxes of DO and ultimately result in a higher net horizontal advective flux of 0.21 $\mathrm{mg}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{day}}^{-\mathrm{1}}$. Consequently, the net effect of the SOD on the upper DO is 0.15 $\mathrm{mg}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{day}}^{-\mathrm{1}}$, which causes a decline of 2.22 mg L⁻¹ in DO concentrations in the surface layer. Figure 9 shows contributions of the SOD and the water column production rates to the changes in surface DO. The positive values of Δ(−DO_SOD) and Δ(DO_WCP) represent the increased DO concentrations due to the decrease in the SOD and increase in the water column production. In the RivPOC−50% simulation, decreasing the POC input decreases the SOD rate but increases the water column production rate, which in combination increase the DO concentrations in the surface layer. As a result, the re-aeration in the RivPOC−50% simulation is weakened, especially in the west of the lower estuary (Fig. 10a).

Figure 10 The changes in air–sea re-aeration rates (a, d), −DO_SOD (b, e), and DO_WCP (c, f) concentrations in the surface layer with respect to the Base simulation (Base model). Positive values of ΔRea, Δ(−DO_SOD), and ΔDO_WCP concentrations represent higher re-aeration rates and higher DO concentrations caused by the changes in the sediment oxygen demand rate and the water column production rate, respectively.

Figure 11 Conceptual schematic of the oxygen dynamics in response to riverine inputs in the PRE. The white boxes represent the state variables in the water column, the orange boxes represent the source and sink processes associated with the oxygen dynamics. The positive signs represent the sources, while the negative signs represent the sinks for DO concentrations.

4 Discussion

4.1 Comparability of hypoxia in 2006

In this study, we performed a series of numerical experiments together with the application of the DO species-tracing method to study the effects of different anthropogenic inputs on hypoxia and oxygen dynamics in the PRE. This study is the first attempt to quantitatively estimate the interactions between each DO source and sink processes (e.g. DO buffering effects) under the anthropogenic perturbations in the PRE. The year 2006 was selected because distinct hypoxia was observed, and the available observations are relatively more abundant than in other years. In addition, it is a wet year with the annual averaged total river discharge over 10 000 m³ s⁻¹ (interannual variations of total discharges during 1999–2010 in the PRE can be seen in Fig. S1 in the Supplement). Discussions only focus on hypoxia in July and August 2006 when oxygen observations are available. However, conclusions drawn here should be applicable to other years because previous studies have reported similar locations and spatial extents of hypoxia in other years (Lin et al., 2001; Zhang and Li, 2010). The mechanisms underlying hypoxia in summer 2006 found here are also consistent with previous studies on hypoxia in this region, such as the strong re-aeration (Zhang and Li, 2010), the dominance of the SOD (Yin et al., 2004; Zhang and Li, 2010), and the important contributions of the allochthonous POC (Hu et al., 2006; Yu et al., 2010).

Table 3 A summary of characteristics of hypoxia among three systems (i.e. Chesapeake Bay, northern Gulf of Mexico, and PRE). WCR represents the water column respiration, which is the sum of respiration, nitrification, and oxidation.

4.2 Relative contributions of different anthropogenic inputs

Numerical experiments show that hypoxia in the PRE is more sensitive to the riverine inputs of POC rather than the nutrient loading (Fig. 6a). This is distinct from other hypoxic systems such as the Chesapeake Bay (Hagy et al., 2004) and the NGOM (Justić et al., 2003) that have observed close relation between nutrient loading and hypoxia. We attribute this to the different characteristics of hypoxia in these systems (Table 3). In the Chesapeake Bay, the dominant oxygen sink leading to hypoxia is the water column respiration, which is associated with high primary productivity stimulated by the excessive nutrient loading (Hong and Shen, 2013). In contrast, the bottom water DO depletions are dominated by the SOD in the NGOM (Murrell and Lehrter, 2011; Yu et al., 2015b) and the PRE (Yin et al., 2004; Zhang and Li, 2010). Hypoxia in the NGOM can be well simulated with appropriate parameterization of SOD while neglecting the water column processes (Yu et al., 2015a).

However, the relative contributions of autochthonous POC (i.e. the POC generated by settling of phytoplankton after death) versus allochthonous POC to the SOD are different in the NGOM and the PRE. In the NGOM, the autochthonous POC serves as the major source of POC (Green et al., 2006), which means that increasing the nutrient loading can facilitate the SOD by increasing the depositional fluxes of dead phytoplankton and ultimately promote the formation of hypoxia. In the PRE, the relative contributions of autochthonous versus allochthonous POC inputs to the SOD and hypoxia have long been a topic of debate. Some studies suggest that allochthonous POC dominates in wet seasons due to the high river discharges (Ye et al., 2017; Yu et al., 2010), while others argue that autochthonous inputs can also play an important role (Guo et al., 2015; Su et al., 2017). Previous studies (Guo et al., 2015; Hu et al., 2006; Ye et al., 2017) show that the ratios of allochthonous POC to autochthonous POC have distinct spatial and seasonal variabilities in the PRE. Generally, the allochthonous contributions dominate inside the estuary and gradually decrease seaward as the impact of the river discharges weakens (Hu et al., 2006; Jia and Peng, 2003). In our study, the high-hypoxic-frequency zone is near the Modaomen sub-estuary, which receives high depositional fluxes of allochthonous POC. Therefore, the allochthonous inputs have dominant contributions to the SOD and summer hypoxia in the high-hypoxic-frequency zone.

Table 4 A summary of the differences in physical and biogeochemical processes associated with the relative contributions of autochthonous versus allochthonous POC between the PRE and NGOM.

^a Zhang et al. (2013). ^b Wang et al. (2004). ^c Ye et al. (2014). ^d Quigg et al. (2011). ^e Hu and Li (2009). ^f Rabouille et al. (2008).

The different POC sources in the NGOM and the PRE might be explained by their distinct physical and biogeochemical processes (Table 4). Firstly, the relative magnitudes of autochthonous versus allochthonous POC are different in the two hypoxic systems. The allochthonous inputs of POC in the NGOM and PRE are at the same magnitude: 3.8×10⁶ t yr⁻¹ (Wang et al., 2004) and 2.5×10⁶ t yr⁻¹ (Zhang et al., 2013). However, the autochthonous inputs in the two systems are different. According to our model results, the primary productivity in the PRE is 310.8±427.5 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{day}}^{-\mathrm{1}}$, which is within the range of 183.9–1213 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{day}}^{-\mathrm{1}}$ reported by Ye et al. (2014). However, the observed primary productivity in the NGOM ranges from 330 to 7010 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{day}}^{-\mathrm{1}}$ (Quigg et al., 2011), the upper range of which is much higher than that in the PRE. The relatively lower primary productivity in the PRE is a result of the stronger phosphorus limitation (DIN:DIP ratio of 126 in the PRE versus 33 in the NGOM) and the light shading effects of high suspended sediment concentrations. The dominant role of the allochthonous POC in highly turbid estuaries has been reported in previous studies (Fontugne and Jouanneau, 1987; Middelburg and Herman, 2007). Secondly, the fates of the allochthonous POC in the two systems are different due to the difference in the residence time between the systems. In the PRE, the residence time is 3–5 days during the wet season, which is much shorter than in the NGOM (∼95 days). It follows that the allochthonous POC cannot be degraded completely and hence can significantly fuel the SOD in the PRE. The difference in surface salinity distribution can also be used to explain the different relative roles of allochthonous POC in the two hypoxic systems. Previous studies have suggested a good correlation between the relative contributions of allochthonous POC and the salinity, namely that the contributions of allochthonous POC generally decrease as salinity increases seaward (Fontugne and Jouanneau, 1987; Middelburg and Herman, 2007). Similar correlations have also been reported in the PRE (Yu et al., 2010) and NGOM (Wang et al., 2004). The surface salinity in the high-hypoxia-frequency zone varies between 0 and 10 psu during the wet season based on our model results, while the surface salinity in the hypoxic zone of the NGOM is saltier than 24 psu, even in the wet season, according to the results from a well-validated physical model in Yu et al. (2015a). This implies a more important role of allochthonous POC in the PRE than in the NGOM. Finally, compositions of the allochthonous POC are different in the two hypoxic systems. Zhang and Li (2010) mentioned that contributions of labile POC to the allochthonous POC are higher in the PRE than in the NGOM.

4.3 The importance of re-aeration in PRE

Model results also highlight the importance of re-aeration in regulating DO dynamics and hypoxia migration in the PRE. On the one hand, based on our previous study by applying the same physical–biogeochemical model and tracing method as here (Wang et al., 2017), the re-aeration together with the SOD are the most important process controlling DO dynamics. Nearly 28 % of the surface DO_REA can reach the bottom layer, exerting a strong constraint on the spatial extent and duration of hypoxia in the PRE. When turning off the re-aeration, the high SOD will lead to persistent hypoxia covering an area of over 3000 km² in the PRE. On the other hand, the re-aeration responds rapidly to the perturbations of riverine inputs, which moderates the DO changes impacted by the perturbations. A conceptual diagram of these processes is illustrated in Fig. 11. Compared with other hypoxic systems, the re-aeration in the PRE is of great importance because of the shallow topography and the strong re-aeration, which enable the surface oxygen supplied by re-aeration to penetrate to the bottom water. Re-aeration can thus greatly influence spatial migration of hypoxia under the perturbations of riverine inputs in the PRE. Furthermore, the shallow topography in the PRE allows the bottom SOD to indirectly affect the surface DO by decreasing the upward DO advective fluxes, which also facilitates strong re-aeration in the PRE. As we have described in Sect. 3.3, the bottom SOD can lead to a decrease in surface DO concentrations by 2.22 mg L⁻¹. If turning off the SOD, the surface AOU would change from 0.61 to −1.61 mg L⁻¹, causing a change of re-aeration from 0.55 to −1.45 $\mathrm{mg}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{day}}^{-\mathrm{1}}$. This indicates that the SOD could shift the role of re-aeration from a strong oxygen sink to a strong source.

One counter-example to the shallow PRE is the NGOM, in which the hypoxic zone is deeper such that the surface water and bottom hypoxic water are detached. Also, the observed SOD varies from 0.06 to 0.70 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{day}}^{-\mathrm{1}}$ in the summer season in the NGOM (Murrell and Lehrter, 2011), which is much lower than those in the PRE (0.72–3.89 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{day}}^{-\mathrm{1}}$; Chung et al., 2004). These characteristics together with the supersaturated DO concentrations in the surface water due to the high primary productivity make the re-aeration primarily an outgassing process in the NGOM (Yu et al., 2015b).

In the other hypoxic system, the Chesapeake Bay as described earlier, extended the discussion on the importance of re-aeration is limited by a lack of observations and relevant studies of re-aeration. Nevertheless, according to our results, we can speculate that the re-aeration might be quite important in the Chesapeake Bay because the strong water column respiration can draw down the surface DO concentrations and enhance the re-aeration. However, the penetration of the oxygen supplied by re-aeration to the bottom layer is hard to estimate without applying the DO species-tracing method like in our study or method, similarly to the Chesapeake Bay. In general, more relevant studies are required to examine the role of the re-aeration on hypoxia in the Chesapeake Bay.

5 Conclusion

This study uses a physical–biogeochemical model to simulate the DO dynamics and hypoxia in the PRE and investigate their responses to anthropogenic perturbations in riverine inputs. Model results based on a simulation in 2006 show that hypoxia in the PRE starts in April, peaks in August, and disappears in October. Perturbing riverine inputs has strong impacts on DO dynamics and hypoxia. The hypoxic extent in the PRE is most sensitive to riverine input of particulate organic carbon, followed by oxygen and nutrients. This is different from other hypoxic systems (i.e. NGOM and Chesapeake Bay) because of the distinct physical and biogeochemical features in the PRE, i.e. the shallow topography, high water exchange rates and dominance of the SOD for DO depletion within bottom waters.

Model results also highlight the importance of re-aeration on hypoxia, which has strong buffering effects on the oxygen dynamics in the PRE. River-induced changes in source and sink processes can trigger an opposite shift in re-aerations by altering the surface AOU. In turn, the re-aeration can moderate the DO changes and hypoxia shifts responding to the changes in the oxygen source and sink processes. The important role of re-aeration in the PRE is due to the shallow waters and strong SOD in the estuary. Firstly, because of the shallow topography, the SOD can affect the surface DO indirectly by decreasing the surface AOU and consequently shifting re-aeration from an oxygen sink to a strong source process. Secondly, the shallow waters enable the oxygen supplied by the re-aeration to diffuse to bottom waters and compensate the DO loss by the SOD.

Appendix A: Each component of water column production

The water column production (WCP) used in this study represents the net effects of water column on DO, which is a combination of the photosynthesis, respiration, nitrification, and oxidation:

$$\begin{array}{}\text{(A1)}& \text{WCP}=\text{phot}-\text{resp}-\text{nitrif}-\text{oxid}.\end{array}$$

The first term “phot” represents the photosynthesis ($\mathrm{mg}{\mathrm{O}}_{\mathrm{2}}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{day}}^{-\mathrm{1}}$):

$$\begin{array}{ll}{\displaystyle }\text{Phot}& {\displaystyle }=\left[{\mathit{\alpha }}_{\mathrm{OC}}\cdot {\mathit{\alpha }}_{{\mathrm{NH}}_{\mathrm{4}}}\cdot G_{\mathrm{P}}\cdot P_{\mathrm{c}}+\left({\mathit{\alpha }}_{{\mathrm{NO}}_{\mathrm{23}}\mathrm{c}}\right)\right.\\ \text{(A2)}& {\displaystyle }& {\displaystyle }\left.\cdot \left(\mathrm{1}-{\mathit{\alpha }}_{{\mathrm{NH}}_{\mathrm{4}}}\right)\cdot G_{\mathrm{P}}\cdot P_c\right],\end{array}$$

where α_OC represents the oxygen-to-carbon ratio (mg O₂:mg C), ${\mathit{\alpha }}_{{\mathrm{NH}}_{\mathrm{4}}}$ represents the phytoplankton's preference for ammonium uptake (dimensionless), G_P represents specific phytoplankton growth rate (day⁻¹) which is dependent on the temperature, light, and nutrients (including NO₂+NO₃, NH₄, PO₄, Si; see Eq. 14), P_c represents phytoplankton biomass (mg C L⁻¹), and ${\mathit{\alpha }}_{{\mathrm{NO}}_{\mathrm{23}}\mathrm{c}}$ represents the oxygen-to-carbon ratio for nitrate uptake (mg O₂:mg N).

The term “resp” represents the respiration ($\mathrm{mg}{\mathrm{O}}_{\mathrm{2}}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{day}}^{-\mathrm{1}}$):

$$\begin{array}{}\text{(A3)}& \text{Resp}={\mathit{\alpha }}_{\mathrm{OC}}\cdot k_{\mathrm{PR}}\left(T\right)\cdot P_c,\end{array}$$

where k_PR(T) represents the temperature-dependent respiration rate (day⁻¹).

The term “nitrif” represents the nitrification ($\mathrm{mg}{\mathrm{O}}_{\mathrm{2}}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{day}}^{-\mathrm{1}}$):

$$\begin{array}{}\text{(A4)}& \text{Nitrif}=\mathrm{2}\cdot {\mathit{\alpha }}_{\mathrm{ON}}\cdot k_{\mathrm{14},\mathrm{15}}{\mathit{\theta }}_{\mathrm{14},\mathrm{15}}^{T-\mathrm{20}}\cdot {\mathrm{NH}}_{\mathrm{4}}\cdot {\displaystyle \frac{DO}{K_{\mathrm{nitri}}+\text{DO}}},\end{array}$$

where α_ON represents the oxygen-to-nitrogen ratio (mg O₂:mg N), k_14,15 represents the nitrification rate at 20 ^∘C (day⁻¹), θ_14,15 represents the temperature coefficient (dimensionless), and K_nitri represent the half saturation constant for oxygen limitation (mg O₂ L⁻¹).

The term “oxid” represents the oxidation of dissolved organic carbon and dissolved sulfide ($\mathrm{mg}{\mathrm{O}}_{\mathrm{2}}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{day}}^{-\mathrm{1}}$):

$$\begin{array}{ll}{\displaystyle }\text{Oxid}& {\displaystyle }={\mathit{\alpha }}_{\mathrm{OC}}\cdot \left[k_{\mathrm{20},\mathrm{0}}{\mathit{\theta }}_{\mathrm{20},\mathrm{0}}^{T-\mathrm{20}}\cdot \text{RDOC}+k_{\mathrm{21},\mathrm{0}}{\mathit{\theta }}_{\mathrm{21},\mathrm{0}}^{T-\mathrm{20}}\cdot \text{LDOC}\right.\\ {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{\text{LDOC}}{K_{\mathrm{LDOC}}+\text{LDOC}}}+k_{\mathrm{22},\mathrm{0}}{\mathit{\theta }}_{\mathrm{22},\mathrm{0}}^{T-\mathrm{20}}\cdot \text{ReDOC}\\ {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{\text{ReDOC}}{K_{\mathrm{LDOC}}+\text{ReDOC}}}+k_{\mathrm{23},\mathrm{0}}{\mathit{\theta }}_{\mathrm{23},\mathrm{0}}^{T-\mathrm{20}}\cdot \text{ExDOC}\\ {\displaystyle }& {\displaystyle }\left.\cdot {\displaystyle \frac{\text{ExDOC}}{K_{\mathrm{LDOC}}+\text{ExDOC}}}\right]\cdot {\displaystyle \frac{P_c}{K_{P_c}+P_c}}\cdot {\displaystyle \frac{\text{DO}}{K_{\mathrm{DO}}+\text{DO}}}\\ \text{(A5)}& {\displaystyle }& {\displaystyle }+k_{{\mathrm{O}}_{\mathrm{2}}^{\ast }}{\mathit{\theta }}_{{\mathrm{O}}_{\mathrm{2}}^{\ast }}^{T-\mathrm{20}}\cdot {\mathrm{O}}_{\mathrm{2}}^{\ast }\cdot {\displaystyle \frac{P_c}{K_{P_c}+P_c}}\cdot {\displaystyle \frac{\text{DO}}{K_{{\mathrm{DO}}_{{\mathrm{O}}_{\mathrm{2}}^{\ast }}}+\text{DO}}},\end{array}$$

where k_20,0, k_21,0, k_22,0, k_23,0, and $k_{{\mathrm{O}}_{\mathrm{2}}^{\ast }}$ represent the oxidation rates of refractory dissolved organic carbon (RDOC), labile dissolved organic carbon (LDOC), reactive dissolved organic carbon (ReDOC), algal exudate dissolved organic carbon (ExDOC), and dissolved sulfide at 20 ^∘C (day⁻¹); θ_20,0, θ_21,0, θ_22,0, θ_23,0, and ${\mathit{\theta }}_{{\mathrm{O}}_{\mathrm{2}}^{\ast }}$ represent the temperature coefficient (dimensionless); K_LDOC represents the Michaelis constant for LDOC (mg C L⁻¹); $K_{P_c}$ represents the half-saturation constant for phytoplankton limitation (mg C L⁻¹); K_DO and $K_{{\mathrm{DO}}_{{\mathrm{O}}_{\mathrm{2}}^{\ast }}}$ represent the half-saturation constant for DO limitation (mg O L⁻¹). More detailed information on these variables and parameters can be seen in Tables A1 and A2.

Table A1 List of state variables in the water quality model (RCA) and the DO species-tracing method.

Table A2 Main parameters and constants for the water quality model (RCA).

Appendix B: Sediment flux module

In this study, a sediment flux module is used to receive the depositional fluxes of particulate organic carbon, particulate organic nitrogen, and particulate organic phosphorus, which are collectively referred to as particulate organic matter, from the overlying water. After that, the diagenesis of particulate organic matter will occur in the sediment and produce soluble end products. The fluxes of nutrients and SOD across the water–sediment interface will be determined by differences in the dissolved concentrations between the resulting sediment and overlying water combined with the transfer coefficient.

In the sediment flux module, particulate organic matter is classified into three G classes (G1 is reactive, G2 is refractory, and G3 is inert) with the different reaction rates. The kinetic equation for diagenesis is as follows:

$$\begin{array}{}\text{(B1)}& H{\displaystyle \frac{\mathrm{d}{G}_i}{\mathrm{d}t}}=-k_{G_i}{\mathit{\theta }}_{G_i}^{T-\mathrm{20}}{G}_{i}H+J_{G_i},\end{array}$$

where H is the depth of sediment (m), G represents the particulate organic carbon, the particulate organic nitrogen, or the particulate organic phosphorus (mg L⁻¹), subscript i represents the ith G class (i=1, 2, 3), $k_{G_i}$ represents the corresponding reaction rate (day⁻¹), ${\mathit{\theta }}_{G_i}$ represents the temperature coefficient (dimensionless), and $J_{G_i}$ represents the depositional fluxes of G_i from the overlying water ($\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{day}}^{-\mathrm{1}}$).

After the deposition and diagenesis, further reactions of organic matter (including particulate organic carbon, dissolved organic carbon, particulate organic nitrogen, dissolved organic nitrogen, particulate organic phosphorus, and dissolved organic phosphorus) will occur in both the aerobic layer (denoted layer 1) and anaerobic layer (denoted layer 2). The mass balance equations can be expressed as a general form:

$$\begin{array}{ll}{\displaystyle }{H}_{\mathrm{1}}{\displaystyle \frac{\mathrm{d}{C}_{\mathrm{T}\mathrm{1}}}{\mathrm{d}t}}& {\displaystyle }=K_{\mathrm{L}\mathrm{01}}\left(f_{\mathrm{d}\mathrm{1}}{C}_{\mathrm{T}\mathrm{1}}-C_{\mathrm{d}\mathrm{0}}\right)+w_{\mathrm{12}}\left(f_{\mathrm{p}\mathrm{2}}{C}_{\mathrm{T}\mathrm{2}}-f_{\mathrm{p}\mathrm{1}}{C}_{\mathrm{T}\mathrm{1}}\right)\\ \text{(B2)}& {\displaystyle }& {\displaystyle }+K_{\mathrm{L}\mathrm{12}}\left(f_{\mathrm{d}\mathrm{2}}{C}_{\mathrm{T}\mathrm{2}}-f_{\mathrm{d}\mathrm{1}}{C}_{\mathrm{T}\mathrm{1}}\right)-K_{\mathrm{1}}{H}_{\mathrm{1}}{C}_{\mathrm{T}\mathrm{1}}+J_{\mathrm{T}\mathrm{1}}{\displaystyle },H_{\mathrm{2}}{\displaystyle \frac{\mathrm{d}{C}_{\mathrm{T}\mathrm{2}}}{\mathrm{d}t}}& {\displaystyle }=-w_{\mathrm{12}}\left(f_{\mathrm{p}\mathrm{2}}{C}_{\mathrm{T}\mathrm{2}}-f_{\mathrm{p}\mathrm{1}}{C}_{\mathrm{T}\mathrm{1}}\right)-K_{\mathrm{L}\mathrm{12}}\left(f_{\mathrm{d}\mathrm{2}}{C}_{\mathrm{T}\mathrm{2}}\right.\\ \text{(B3)}& {\displaystyle }& {\displaystyle }\left.-f_{\mathrm{d}\mathrm{1}}{C}_{\mathrm{T}\mathrm{1}}\right)-K_{\mathrm{2}}{H}_{\mathrm{2}}{C}_{\mathrm{T}\mathrm{2}}-w_{\mathrm{2}}{C}_{\mathrm{T}\mathrm{2}}+J_{\mathrm{T}\mathrm{2}},\end{array}$$

where subscript 0, 1, and 2 represent the overlying water; the aerobic layer, and the anaerobic layer; H₁ and H₂ represent the thicknesses of aerobic layer and anaerobic layer (m); C_T1 and C_T2 represent the total concentrations (mg L⁻¹) in aerobic layer and anaerobic layer; C_d0 represents the dissolved concentrations (mg L⁻¹) in the overlying water; f_d1 and f_d2 represent the dissolved fractions in the aerobic layer and anaerobic layer (dimensionless), f_p1 and f_p2 represent the particulate fractions in aerobic layer and anaerobic layer (dimensionless), K_L01 represents the transfer coefficient between the overlying water and aerobic layer (m day⁻¹), K_L12 represents the transfer coefficient between the aerobic layer and anaerobic layer (m day⁻¹), K₁ and K₂ represent the first-order decay rate of in the aerobic layer and anaerobic layer (day⁻¹), w₁₂ represents the particle mixing rate between the aerobic layer and anaerobic layer (m day⁻¹), w₂ represents the sedimentation rate out of the anaerobic layer (m day⁻¹), and J_T1 and J_T2 represent the total influxes for each class of particulate organic matter into the aerobic layer and anaerobic layer ($\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{day}}^{-\mathrm{1}}$).

Fluxes of nutrients and DO across the water–sediment interface can be represented as follows:

$$\begin{array}{}\text{(B4)}& J=s\left(C_{\mathrm{water}}-C_{\mathrm{sed}}\right),\end{array}$$

where s represents the transfer coefficient (m day⁻¹), and C_water and C_sed represent the concentrations of nutrients and DO in the water and sediment ($\mathrm{mg}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{day}}^{-\mathrm{1}}$).