Journal cover Journal topic
Biogeosciences An interactive open-access journal of the European Geosciences Union
Journal topic
Biogeosciences, 15, 6105–6125, 2018
https://doi.org/10.5194/bg-15-6105-2018
Biogeosciences, 15, 6105–6125, 2018
https://doi.org/10.5194/bg-15-6105-2018

Research article 18 Oct 2018

Research article | 18 Oct 2018

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

Impacts of anthropogenic inputs on hypoxia and oxygen dynamics in the Pearl River estuary
Bin Wang1,3, Jiatang Hu1,2, Shiyu Li1,2, Liuqian Yu3, and Jia Huang1 Bin Wang et al.
• 1School of Environmental Science and Engineering, Sun Yat-Sen University, Guangzhou, 510275, China
• 2Guangdong Provincial Key Laboratory of Environmental Pollution Control and Remediation Technology, Guangzhou 510275, China
• 3Department of Oceanography, Dalhousie University, Halifax, Nova Scotia, Canada

Correspondence: Jiatang Hu (hujtang@mail.sysu.edu.cn) and Shiyu Li (eeslsy@mail.sysu.edu.cn)

Abstract

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 < 2mg L−1) 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 km2. 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×1011m3 yr−1 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×105t yr−1 of dissolved inorganic nitrogen (DIN) and 9.9×103t yr−1 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}×{\mathrm{10}}^{\mathrm{6}}$t yr−1 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

### 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 2Conceptual 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, NH4 represents ammonia nitrogen, NO23 represents nitrite and nitrate nitrogen, PON represents particulate organic nitrogen, DON represents dissolved organic nitrogen, DPO4 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 O2 L−1) is given by the following:

$\begin{array}{ll}& \frac{\partial \text{DO}}{\partial t}+u\frac{\partial \text{DO}}{\partial x}+v\frac{\partial \text{DO}}{\partial y}+w\frac{\partial \text{DO}}{\partial z}-\frac{\partial }{\partial x}\left({E}_{x}\frac{\partial \text{DO}}{\partial x}\right)\\ & -\frac{\partial }{\partial y}\left({E}_{y}\frac{\partial \text{DO}}{\partial y}\right)-\frac{\partial }{\partial z}\left({E}_{z}\frac{\partial \text{DO}}{\partial z}\right)\\ \text{(1)}& & =\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−1) represent velocity components in x, y, and z coordinates; and Ex, Ey, and Ez (m s−2) 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}\phantom{\rule{0.125em}{0ex}}{\mathrm{O}}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{L}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\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}\phantom{\rule{0.125em}{0ex}}{\mathrm{O}}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{L}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\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 DOsat represents the DO concentration at saturation (mg O2 L−1), which is dependent on salinity and temperature; ka is the surface mass transfer coefficient (day−1); 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}\phantom{\rule{0.125em}{0ex}}{\mathrm{O}}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{L}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\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}=\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−1), and DOsed represents DO concentrations (mg O2 L−1) 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 (DOBC (mg O2 L−1)), air–sea re-aeration (DOREA (mg O2 L−1)), water column production (DOWCP (mg O2 L−1)), and SOD (DOSOD (mg O2 L−1)) (Table A1). Equations of the four numerical oxygen species are given below:

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

where “tran” represents the physical transport processes, i.e. advection $\left(u\frac{\partial }{\partial x}+v\frac{\partial }{\partial y}+w\frac{\partial }{\partial z}\right)$ and diffusion $\left(-\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)\right)$. REA ($\mathrm{mg}\phantom{\rule{0.125em}{0ex}}{\mathrm{O}}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{L}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{day}}^{-\mathrm{1}}$), WCP ($\mathrm{mg}\phantom{\rule{0.125em}{0ex}}{\mathrm{O}}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{L}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{day}}^{-\mathrm{1}}$), and SOD ($\mathrm{mg}\phantom{\rule{0.125em}{0ex}}{\mathrm{O}}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{L}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\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 DOBC concentrations are only controlled by advection and diffusion. By assigning the initial conditions and lateral boundary conditions of DOBC the same as those for DO, the mixing curve of DOBC will overlap the theory mixing curve shown in Fig. 3b. It follows that the DOBC represents the conservative effects, while the DOREA, DOWCP, and DOSOD 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, RivDO50%, RivDO+50%, RivNtr50%, RivNtr+50%, RivPOC50% and RivPOC+50%, with the basic information of each simulation presented in Table 1.

Table 1Overview 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−1 based on observations nearby. The RivDO50% 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 RivNtr50% 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)}& & \text{Hypoxic area}=\sum p\cdot \mathrm{\Delta }s\text{(11)}& & \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=\frac{{N}_{\mathrm{h}}}{{N}_{\mathrm{s}}}\cdot \mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathit{%},\end{array}$

where Nh is the number of hours in which hypoxia occurs, and Ns is the total number of hours. In this study, the threshold of hypoxia is defined as 3 mg L−1 (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 5Spatial distributions of bottom DO (a–f) and hypoxic frequency (g–l) during May–October. Hypoxia is defined as DO concentration below 3 mg L−1.

Figure 6Spatial 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 (R2=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 RivDO50% 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 RivDO50% 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 RivDO50%. 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 7The 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 8The spatial distribution of DOBC (a, b, c) and DOREA (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 DOBC and DOREA in the bottom water are shown in Fig. 8. Differences in DOWCP and DOSOD concentrations between simulations are much smaller and hence omitted here.

Halving the riverine DO inputs in the RivDO50% simulation yields lower DOBC concentrations but higher DOREA in the bottom water (Figs. 7b and 8). The decrease in DOBC concentrations is largely balanced by the increase in DOREA concentration in RivDO50% 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 DOBC concentrations (Figs. 7b and 8c) but lower DOREA 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.42mg L−1, while the surface domain-averaged DO concentration in the Base and RivDO50% simulations are 6.81 and 6.57 mg L−1. The surface AOU for the RivDO50% simulation is 39 % higher than that for the Base simulation, which is consistent with the 38 % increase in re-aeration rate for the RivDO50% 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 RivNtr50% 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 DOREA 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 RivNtr50% 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 RivPOC50% simulation, the DO concentration increases by 0.56 mg L−1 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 RivPOC50% 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−1 in the bottom waters of the high-hypoxic-frequency zone (Fig. 7d). However, decreasing the riverine inputs of POC in the RivPOC50% 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 GP (day−1) 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 GPmax represents the maximum growth rate at the optimum conditions (day−1); 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 2Comparisons 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 9Budget of −DOSOD 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}\phantom{\rule{0.125em}{0ex}}{\mathrm{L}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\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 DOSOD 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}\phantom{\rule{0.125em}{0ex}}{\mathrm{L}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{day}}^{-\mathrm{1}}$ and decreases the upward advective DO fluxes, reaching the upper layer by 0.34 $\mathrm{mg}\phantom{\rule{0.125em}{0ex}}{\mathrm{L}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\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}\phantom{\rule{0.125em}{0ex}}{\mathrm{L}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\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}\phantom{\rule{0.125em}{0ex}}{\mathrm{L}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{day}}^{-\mathrm{1}}$. Consequently, the net effect of the SOD on the upper DO is 0.15 $\mathrm{mg}\phantom{\rule{0.125em}{0ex}}{\mathrm{L}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{day}}^{-\mathrm{1}}$, which causes a decline of 2.22 mg L−1 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 Δ(−DOSOD) and Δ(DOWCP) represent the increased DO concentrations due to the decrease in the SOD and increase in the water column production. In the RivPOC50% 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 RivPOC50% simulation is weakened, especially in the west of the lower estuary (Fig. 10a).

Figure 10The changes in air–sea re-aeration rates (a, d), −DOSOD (b, e), and DOWCP (c, f) concentrations in the surface layer with respect to the Base simulation (Base model). Positive values of ΔRea, Δ(−DOSOD), and ΔDOWCP 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 11Conceptual 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 m3 s−1 (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 3A 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 4A 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×106t yr−1 (Wang et al., 2004) and 2.5×106t yr−1 (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}\phantom{\rule{0.125em}{0ex}}\mathrm{C}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{day}}^{-\mathrm{1}}$, which is within the range of 183.9–1213 $\mathrm{mg}\phantom{\rule{0.125em}{0ex}}\mathrm{C}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\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}\phantom{\rule{0.125em}{0ex}}\mathrm{C}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\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 DOREA 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 km2 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−1. If turning off the SOD, the surface AOU would change from 0.61 to −1.61mg L−1, causing a change of re-aeration from 0.55 to −1.45$\mathrm{mg}\phantom{\rule{0.125em}{0ex}}{\mathrm{L}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\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}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\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}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\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.

Data availability
Data availability.

All data used to generate figures can be assessed publicly at http://wemt_number1.pc.goabc.cn/wemt_number1/vip_doc/10941675.html, last access: October 2018.

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}\phantom{\rule{0.125em}{0ex}}{\mathrm{O}}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{L}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{day}}^{-\mathrm{1}}$):

$\begin{array}{ll}\text{Phot}& =\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)}& & \cdot \left(\mathrm{1}-{\mathit{\alpha }}_{{\mathrm{NH}}_{\mathrm{4}}}\right)\cdot {G}_{\mathrm{P}}\cdot {P}_{c}],\end{array}$

where αOC represents the oxygen-to-carbon ratio (mg O2:mg C), ${\mathit{\alpha }}_{{\mathrm{NH}}_{\mathrm{4}}}$ represents the phytoplankton's preference for ammonium uptake (dimensionless), GP represents specific phytoplankton growth rate (day−1) which is dependent on the temperature, light, and nutrients (including NO2+NO3, NH4, PO4, Si; see Eq. 14), Pc represents phytoplankton biomass (mg C L−1), and ${\mathit{\alpha }}_{{\mathrm{NO}}_{\mathrm{23}}\mathrm{c}}$ represents the oxygen-to-carbon ratio for nitrate uptake (mg O2:mg N).

The term “resp” represents the respiration ($\mathrm{mg}\phantom{\rule{0.125em}{0ex}}{\mathrm{O}}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{L}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\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 kPR(T) represents the temperature-dependent respiration rate (day−1).

The term “nitrif” represents the nitrification ($\mathrm{mg}\phantom{\rule{0.125em}{0ex}}{\mathrm{O}}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{L}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\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 \frac{DO}{{K}_{\mathrm{nitri}}+\text{DO}},\end{array}$

where αON represents the oxygen-to-nitrogen ratio (mg O2:mg N), k14,15 represents the nitrification rate at 20 C (day−1), θ14,15 represents the temperature coefficient (dimensionless), and Knitri represent the half saturation constant for oxygen limitation (mg O2 L−1).

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

$\begin{array}{ll}\text{Oxid}& ={\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\\ & \cdot \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}\\ & \cdot \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}\\ & \cdot \frac{\text{ExDOC}}{{K}_{\mathrm{LDOC}}+\text{ExDOC}}]\cdot \frac{{P}_{c}}{{K}_{{P}_{c}}+{P}_{c}}\cdot \frac{\text{DO}}{{K}_{\mathrm{DO}}+\text{DO}}\\ \text{(A5)}& & +{k}_{{\mathrm{O}}_{\mathrm{2}}^{\ast }}{\mathit{\theta }}_{{\mathrm{O}}_{\mathrm{2}}^{\ast }}^{T-\mathrm{20}}\cdot {\mathrm{O}}_{\mathrm{2}}^{\ast }\cdot \frac{{P}_{c}}{{K}_{{P}_{c}}+{P}_{c}}\cdot \frac{\text{DO}}{{K}_{{\mathrm{DO}}_{{\mathrm{O}}_{\mathrm{2}}^{\ast }}}+\text{DO}},\end{array}$

where k20,0, k21,0, k22,0, k23,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−1); θ20,0, θ21,0, θ22,0, θ23,0, and ${\mathit{\theta }}_{{\mathrm{O}}_{\mathrm{2}}^{\ast }}$ represent the temperature coefficient (dimensionless); KLDOC represents the Michaelis constant for LDOC (mg C L−1); ${K}_{{P}_{c}}$ represents the half-saturation constant for phytoplankton limitation (mg C L−1); KDO and ${K}_{{\mathrm{DO}}_{{\mathrm{O}}_{\mathrm{2}}^{\ast }}}$ represent the half-saturation constant for DO limitation (mg O L−1). More detailed information on these variables and parameters can be seen in Tables A1 and A2.

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

Table A2Main 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\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−1), subscript i represents the ith G class (i=1, 2, 3), ${k}_{{G}_{i}}$ represents the corresponding reaction rate (day−1), ${\mathit{\theta }}_{{G}_{i}}$ represents the temperature coefficient (dimensionless), and ${J}_{{G}_{i}}$ represents the depositional fluxes of Gi from the overlying water ($\mathrm{g}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\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}{H}_{\mathrm{1}}\frac{\mathrm{d}{C}_{\mathrm{T}\mathrm{1}}}{\mathrm{d}t}& ={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)}& & +{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}},{H}_{\mathrm{2}}\frac{\mathrm{d}{C}_{\mathrm{T}\mathrm{2}}}{\mathrm{d}t}& =-{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)}& & -{f}_{\mathrm{d}\mathrm{1}}{C}_{\mathrm{T}\mathrm{1}})-{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; H1 and H2 represent the thicknesses of aerobic layer and anaerobic layer (m); CT1 and CT2 represent the total concentrations (mg L−1) in aerobic layer and anaerobic layer; Cd0 represents the dissolved concentrations (mg L−1) in the overlying water; fd1 and fd2 represent the dissolved fractions in the aerobic layer and anaerobic layer (dimensionless), fp1 and fp2 represent the particulate fractions in aerobic layer and anaerobic layer (dimensionless), KL01 represents the transfer coefficient between the overlying water and aerobic layer (m day−1), KL12 represents the transfer coefficient between the aerobic layer and anaerobic layer (m day−1), K1 and K2 represent the first-order decay rate of in the aerobic layer and anaerobic layer (day−1), w12 represents the particle mixing rate between the aerobic layer and anaerobic layer (m day−1), w2 represents the sedimentation rate out of the anaerobic layer (m day−1), and JT1 and JT2 represent the total influxes for each class of particulate organic matter into the aerobic layer and anaerobic layer ($\mathrm{g}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\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−1), and Cwater and Csed represent the concentrations of nutrients and DO in the water and sediment ($\mathrm{mg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{3}}\phantom{\rule{0.125em}{0ex}}{\mathrm{day}}^{-\mathrm{1}}$).

Supplement
Supplement.

Author contributions
Author contributions.

JH developed the 1-D–3-D coupled physical model. BW developed the DO-species tracing method and performed the model simulations. The work was supervised by SL and JH. BW wrote the paper with contributions from all co-authors. In addition, all co-authors contributed to the revision of the manuscript.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

This work was supported by the National Natural Science Foundation of China (grant no. 41306105), the Guangdong Natural Science Foundation (grant no. 2014A030313169), the Science and Technology Planning Project of Guangdong Province, China (grant no. 2014A020217003), and the Fundamental Research Funds for the Central Universities (grant no. 17lgzd20).

Edited by: Silvio Pantoja
Reviewed by: two anonymous referees

References

Chung, J., Gary, C., and Heinke, W.: Pearl River Estuary Pollution Project (PREPP), Cont. Shelf Res., 24, 1739–1744, 2004.

Diaz, R. J. and Rosenberg, R.: Spreading Dead Zones and Consequences for Marine Ecosystems, Science, 321, 926–929, https://doi.org/10.1126/science.1156401, 2008.

Du, J. and Shen, J.: Decoupling the influence of biological and physical processes on the dissolved oxygen in the Chesapeake Bay, J. Geophys. Res.-Ocean., 120, 78–93, https://doi.org/10.1002/2014JC010422, 2015.

Fontugne, M. R. and Jouanneau, J.-M.: Modulation of the particulate organic carbon flux to the ocean by a macrotidal estuary: Evidence from measurements of carbon isotopes in organic matter from the Gironde system, Estuar. Coast. Shelf Sci., 24, 377–387, https://doi.org/10.1016/0272-7714(87)90057-6, 1987.

Forrest, D. R., Hetland, R. D., and Dimarco, S. F.: Multivariable statistical regression models of the areal extent of hypoxia over the Texas-Louisiana continental shelf, Environ. Res. Lett., 6, 1–10, https://doi.org/10.1088/1748-9326/6/4/045002, 2011.

Grantham, B. A., Chan, F., Nielsen, K. J., Fox, D. S., Barth, J. A., Huyer, A., Lubchenco, J., and Menge, B. A.: Upwelling-driven nearshore hypoxia signals ecosystem and oceanographic changes in the northeast Pacific, Nature, 429, 749–754, https://doi.org/10.1038/nature02605, 2004.

Green, R. E., Bianchi, T. S., Dagg, M. J., Walker, N. D., and Breed, G. A.: An organic carbon budget for the Mississippi River turbidity plume and plume contributions to air-sea CO2 fluxes and bottom water hypoxia, Estuar. Coast., 29, 579–597, https://doi.org/10.1007/BF02784284, 2006.

Guo, W., Ye, F., Xu, S., and Jia, G.: Seasonal variation in sources and processing of particulate organic carbon in the Pearl River estuary, South China, Estuar. Coast. Shelf Sci., 167, 540–548, https://doi.org/10.1016/j.ecss.2015.11.004, 2015.

Hagy, J. D., Boynton, W. R., Keefe, C. W., and Wood, K. V.: Hypoxia in Chesapeake Bay, 1950–2001: Long-term change in relation to nutrient loading and river flow, Estuaries, 27, 634–658, https://doi.org/10.1007/BF02907650, 2004.

He, B., Dai, M., Zhai, W., Guo, X., and Wang, L.: Hypoxia in the upper reaches of the Pearl River Estuary and its maintenance mechanisms: A synthesis based on multiple year observations during 2000–2008, Mar. Chem., 167, 13–24, https://doi.org/10.1016/j.marchem.2014.07.003, 2014.

Hong, B. and Shen, J.: Linking dynamics of transport timescale and variations of hypoxia in the Chesapeake Bay, J. Geophys. Res.-Ocean., 118, 6017–6029, https://doi.org/10.1002/2013JC008859, 2013.

Hu, J. and Li, S.: Modeling the mass fluxes and transformations of nutrients in the Pearl River Delta, China, J. Mar. Syst., 78, 146–167, https://doi.org/10.1016/j.jmarsys.2009.05.001, 2009.

Hu, J., Peng, P., Jia, G., Mai, B., and Zhang, G.: Distribution and sources of organic carbon, nitrogen and their isotopes in sediments of the subtropical Pearl River estuary and adjacent shelf, Southern China, Mar. Chem., 98, 274–285, https://doi.org/10.1016/j.marchem.2005.03.008, 2006.

Hu, J., Li, S., and Geng, B.: Modeling the mass flux budgets of water and suspended sediments for the river network and estuary in the Pearl River Delta, China, J. Mar. Syst., 88, 252–266, https://doi.org/10.1016/j.jmarsys.2011.05.002, 2011.

HydroQual, Inc.: A Primer for ECOMSED Version 1.3, HydroQual, Inc., Mahwah, NJ, 2002.

HydroQual, Inc.: User's Guide for RCA (Release 3.0), HydroQual, Inc., Mahwah, NJ, 2004.

Jia, G.-D. and Peng, P.-A.: Temporal and spatial variations in signatures of sedimented organic matter in Lingding Bay (Pearl estuary), southern China, Mar. Chem., 82, 47–54, https://doi.org/10.1016/S0304-4203(03)00050-1, 2003.

Justić, D., Rabalais, N. N., and Turner, R. E.: Simulated responses of the Gulf of Mexico hypoxia to variations in climate and anthropogenic nutrient loading, J. Mar. Syst., 42, 115–126, https://doi.org/10.1016/S0924-7963(03)00070-8, 2003.

Laurent, A., Fennel, K., Ko, D., and Lehrter, J.: Climate Change Projected to Exacerbate Impacts of Coastal Eutrophication in the Northern Gulf of Mexico, J. Geophys. Res.-Ocean., 123, 3408–3426, https://doi.org/10.1002/2017JC013583, 2018.

Li, X., Bianchi, T. S., Yang, Z., Osterman, L. E., Allison, M. A., DiMarco, S. F., and Yang, G.: Historical trends of hypoxia in Changjiang River estuary: Applications of chemical biomarkers and microfossils, J. Mar. Syst., 86, 57–68, https://doi.org/10.1016/j.jmarsys.2011.02.003, 2011.

Lin, H. Y., Liu, S., and Han, W. Y.: Potential Trigger CTB, from Seasonal Bottom Water Hypoxia in the Pearl River Estuary, J. Zhanjiang Ocean Univ., 21, 25–29, 2001.

Liu, D., Hu, J., Li, S., and Huang, J.: Validation and application of a three-dimensional coupled water quality and sediment model of the Pearl River Estuary, Huanjing Kexue Xuebao/Acta Sci. Circumstantiae, 36, 4025–4036, https://doi.org/10.13671/j.hjkxxb.2016.0145, 2016.

Luo, L., Li, S. Y., and Wang, D. X.: Modelling of hypoxia in the Pearl River estuary in summer, Adv. Water Sci., 19, 729–735, 2008.

Mellor, G. L. and Yamada, T.: Development of a turbulence closure model for geophysical fluid problems, Rev. Geophys., 20, 851–875, https://doi.org/10.1029/RG020i004p00851, 1982.

Middelburg, J. J. and Herman, P. M. J.: Organic matter processing in tidal estuaries, Mar. Chem., 106, 127–147, https://doi.org/10.1016/J.MARCHEM.2006.02.007, 2007.

Montes, I., Dewitte, B., Gutknecht, E., Paulmier, A., Dadou, I., Oschlies, A., and Garçon, V.: High-resolution modeling of the Eastern Tropical Pacific oxygen minimum zone: Sensitivity to the tropical oceanic circulation, J. Geophys. Res.-Ocean., 119, 5515–5532, https://doi.org/10.1002/2014JC009858, 2014.

Murrell, M. C. and Lehrter, J. C.: Sediment and Lower Water Column Oxygen Consumption in the Seasonally Hypoxic Region of the Louisiana Continental Shelf, Estuar. Coast., 34, 912–924, https://doi.org/10.1007/s12237-010-9351-9, 2011.

Ning, X., Lin, C., Su, J., Liu, C., Hao, Q., and Le, F.: Long-term changes of dissolved oxygen, hypoxia, and the responses of the ecosystems in the East China Sea from 1975 to 1995, J. Oceanogr., 67, 59–75, https://doi.org/10.1007/s10872-011-0006-7, 2011.

Ou, S., Zhang, H., and Wang, D. X.: Dynamics of the buoyant plume off the Pearl River Estuary in summer, Environ. Fluid Mech., 9, 471–492, https://doi.org/10.1007/s10652-009-9146-3, 2009.

Quigg, A., Sylvan, J. B., Gustafson, A. B., Fisher, T. R., Oliver, R. L., Tozzi, S., and Ammerman, J. W.: Going West: Nutrient Limitation of Primary Production in the Northern Gulf of Mexico and the Importance of the Atchafalaya River, Aquat. Geochemistry, 17, 519–544, https://doi.org/10.1007/s10498-011-9134-3, 2011.

Rabalais, N. N., Díaz, R. J., Levin, L. A., Turner, R. E., Gilbert, D., and Zhang, J.: Dynamics and distribution of natural and human-caused hypoxia, Biogeosciences, 7, 585–619, https://doi.org/10.5194/bg-7-585-2010, 2010.

Rabouille, C., Conley, D. J., Dai, M. H., Cai, W. J., Chen, C. T. A., Lansard, B., Green, R., Yin, K., Harrison, P. J., Dagg, M., and McKee, B.: Comparison of hypoxia among four river-dominated ocean margins: The Changjiang (Yangtze), Mississippi, Pearl, and Rhône rivers, Cont. Shelf Res., 28, 1527–1537, https://doi.org/10.1016/j.csr.2008.01.020, 2008.

Scavia, D., Kelly, E. L. A., and Hagy, J. D.: A simple model for forecasting the effects of nitrogen loads on Chesapeake Bay hypoxia, Estuar. Coast., 29, 674–684, https://doi.org/10.1007/BF02784292, 2006.

Smagorinsky, J.: General circulation experiments with the primitive equations: I. The basic equations, Mon. Weather Rev., 91, 99–164, 1963.

Strokal, M., Kroeze, C., Li, L., Luan, S., Wang, H., Yang, S., and Zhang, Y.: Increasing dissolved nitrogen and phosphorus export by the Pearl River (Zhujiang): a modeling approach at the sub-basin scale to assess effective nutrient management, Biogeochemistry, 125, 221–242, https://doi.org/10.1007/s10533-015-0124-1, 2015.

Su, J., Dai, M., He, B., Wang, L., Gan, J., Guo, X., Zhao, H., and Yu, F.: Tracing the origin of the oxygen-consuming organic matter in the hypoxic zone in a large eutrophic estuary: the lower reach of the Pearl River Estuary, China, Biogeosciences, 14, 4085–4099, https://doi.org/10.5194/bg-14-4085-2017, 2017.

Wang, B.: Hydromorphological mechanisms leading to hypoxia off the Changjiang estuary, Mar. Environ. Res., 67, 53–58, https://doi.org/10.1016/j.marenvres.2008.11.001, 2009.

Wang, B., Wei, Q., Chen, J., and Xie, L.: Annual cycle of hypoxia off the Changjiang (Yangtze River) Estuary, Mar. Environ. Res., 77, 1–5, https://doi.org/10.1016/j.marenvres.2011.12.007, 2012.

Wang, B., Hu, J., Li, S., and Liu, D.: A numerical analysis of biogeochemical controls with physical modulation on hypoxia during summer in the Pearl River estuary, Biogeosciences, 14, 2979–2999, https://doi.org/10.5194/bg-14-2979-2017, 2017.

Wang, X.-C., Chen, R. F., and Gardner, G. B.: Sources and transport of dissolved and particulate organic carbon in the Mississippi River estuary and adjacent coastal waters of the northern Gulf of Mexico, Mar. Chem., 89, 241–256, 2004.

Ye, F., Guo, W., Shi, Z., Jia, G., and Wei, G.: Seasonal dynamics of particulate organic matter and its response to flooding in the Pearl River Estuary, China, revealed by stable isotope (δ13C and δ15N) analyses, J. Geophys. Res.-Ocean., 122, 1–22, https://doi.org/10.1002/2017JC012931, 2017.

Ye, H., Chen, C., Sun, Z., Tang, S., Song, X., Yang, C., Tian, L., and Liu, F.: Estimation of the Primary Productivity in Pearl River Estuary Using MODIS Data, Estuar. Coast., 38, 506–518, https://doi.org/10.1007/s12237-014-9830-5, 2014.

Yin, K., Lin, Z., and Ke, Z.: Temporal and spatial distribution of dissolved oxygen in the Pearl River Estuary and adjacent coastal waters, Cont. Shelf Res., 24, 1935–1948, https://doi.org/10.1016/j.csr.2004.06.017, 2004.

Yu, F., Zong, Y., Lloyd, J. M., Huang, G., Leng, M. J., Kendrick, C., Lamb, A. L., and Yim, W. W. S.: Bulk organic δ13C and C/N as indicators for sediment sources in the Pearl River delta and estuary, southern China, Estuar. Coast. Shelf Sci., 87, 618–630, https://doi.org/10.1016/j.ecss.2010.02.018, 2010.

Yu, L., Fennel, K., and Laurent, A.: A modeling study of physical controls on hypoxia generation in the northern Gulf of Mexico, J. Geophys. Res.-Ocean., 120, 5019–5039, https://doi.org/10.1002/2014JC010634, 2015a.

Yu, L., Fennel, K., Laurent, A., Murrell, M. C., and Lehrter, J. C.: Numerical analysis of the primary processes controlling oxygen dynamics on the Louisiana shelf, Biogeosciences, 12, 2063–2076, https://doi.org/10.5194/bg-12-2063-2015, 2015b.

Zhang, H. and Li, S.: Effects of physical and biochemical processes on the dissolved oxygen budget for the Pearl River Estuary during summer, J. Mar. Syst., 79, 65–88, https://doi.org/10.1016/j.jmarsys.2009.07.002, 2010.

Zhang, L., Qin, X., Yang, H., Huang, Q., and Liu, P.: Transported Fluxes of the Riverine Carbon and Seasonal Variation in Pearl River Basin, Environ. Sci., 34, 3025–3034, https://doi.org/10.13227/j.hjkx.2013.08.043, 2013.