Destruction and reinstatement of coastal hypoxia in the South China Sea off the Pearl River estuary

We examined the evolution of intermittent hypoxia off the Pearl River estuary based on three cruise legs conducted in July 2018: one during severe hypoxic conditions before the passage of a typhoon and two post-typhoon legs showing destruction of the hypoxia and its reinstatement. The lowest ever recorded regional dissolved oxygen (DO) concentration of 3.5 µmol kg⁻¹ (∼ 0.1 mg L⁻¹) was observed in bottom waters during leg 1, with an ∼ 660 km² area experiencing hypoxic conditions (DO < 63 µmol kg⁻¹). Hypoxia was completely destroyed by the typhoon passage but was quickly restored ∼ 6 d later, resulting primarily from high biochemical oxygen consumption in bottom waters that averaged 14.6 ± 4.8 µmol O₂ kg⁻¹ d⁻¹. The shoreward intrusion of offshore subsurface waters contributed to an additional 8.6 ± 1.7 % of oxygen loss during the reinstatement of hypoxia. Freshwater inputs suppressed wind-driven turbulent mixing, stabilizing the water column and facilitating the hypoxia formation. The rapid reinstatement of summer hypoxia has a shorter timescale than the water residence time, which is however comparable with that of its initial disturbance from frequent tropical cyclones that occur throughout the wet season. This has important implications for better understanding the intermittent nature of hypoxia and predicting coastal hypoxia in a changing climate.

1 Introduction

Coastal hypoxia has been increasingly exacerbated near the mouths of large rivers as a consequence of anthropogenic nutrient inputs (Gilbert et al., 2010; Rabalais et al., 2014; Breitburg et al., 2018). The rise in the size, intensity and frequency of eutrophication-induced hypoxia exposes coastal oceans to a higher risk of elevated N₂O and CH₄ production; enhanced ocean acidification; and associated reductions in biodiversity, shifts in community structures, and negative impacts on food security and livelihoods (Diaz and Rosenberg, 2008; Vaquer-Sunyer and Duarte, 2008; Naqvi et al., 2010).

Coastal hypoxia can be intermittent due to the dynamic nature of estuarine and coastal environments, where winds, tides, river discharge and circulation patterns strongly affect the ventilation of oxygen-deficient waters (Wang and Justić, 2009; Lu et al., 2018; Zhang et al., 2019). Constraints on oxygen supply can be easily eroded by changes in physical forcings, leading to the temporal alleviation of hypoxia (Laurent and Fennel, 2019). Despite the wide application of oxygen budget analysis and modelling to diagnose the dominant processes driving the formation and maintenance of hypoxia (Yu et al., 2015; Li et al., 2016; Lu et al., 2018), the evolution of intermittent hypoxia, such as the destruction and reinstatement of hypoxia from disturbance by tropical cyclones, remains to be better characterized (Testa et al., 2017). Specifically, the identification of key processes and timescale constraints for these hypoxia destruction and recovery processes is of critical importance in order to predict site-specific hypoxia and its cascading effects and to forecast the long-term impact of hypoxia under a changing climate with higher-intensity extreme events (Knutson et al., 2010; Mendelsohn et al., 2012).

Large riverine nutrient loadings and the resulting eutrophication have recently tipped the lower Pearl River estuary (PRE) and adjacent shelf areas into seasonally hypoxic systems (Yin et al., 2004; Rabouille et al., 2008; Su et al., 2017; Qian et al., 2018; Cui et al., 2019; Zhao et al., 2020). Modelling results have shown that summer hypoxia off the PRE is largely intermittent owing to high-frequency variations in wind forcing and tidal fluctuations (Wang et al., 2017b; Huang et al., 2019). Hypoxia is often interrupted by the passage of typhoons but redevelops quickly with a tendency towards rapid oxygen declines (Su et al., 2017; Huang et al., 2019). The prevailing southwest monsoon usually favours the expansion of a quasi-steady-state freshwater bulge outside the entrance of the PRE (Gan et al., 2009b; Lu et al., 2018) that promotes water column stability. However, it remains unclear how the interaction between wind stress, tidal forcing and freshwater buoyancy affects the bottom oxygen conditions when the winds shift in the downwelling-favourable easterly or southeasterly direction, especially in the wake of tropical cyclones.

It is known that aerobic respiration of organic matter is largely responsible for the oxygen depletion here (Su et al., 2017; Qian et al., 2018). Common methods to quantify oxygen sinks under hypoxic conditions include the analysis of oxygen budgets based on the mass balance of oxygen and estimates of community/bacterial respiration or nitrification rates using field incubations (Zhang and Li, 2010; Li et al., 2015; Su et al., 2017; Cui et al., 2019). However, the budget analysis of oxygen typically assumes a steady-state system (Zhang and Li, 2010; Cui et al., 2019) given that the change of oxygen over time is often much smaller than the oxygen depletion and advection or diffusion fluxes (Cui et al., 2019). The respiration or nitrification rates estimated from (enriched) incubation experiments also merely indicate the oxygen consumption rate (OCR) under a specific low-oxygen condition at the time of sampling (He et al., 2014; Su et al., 2017). The actual magnitude of in situ OCR during the hypoxia formation has rarely been reported at large scales. The role of lateral advection (and/or upwelling) also remains to be better quantified (Zhang and Li, 2010; Lu et al., 2018; Qian et al., 2018; Cui et al., 2019).

We investigated the destruction and reinstatement of summer hypoxia off the PRE to examine the effects of freshwater inputs and winds on water column stability and the maintenance, destruction and formation of hypoxia. With the aid of a three-endmember mixing model, we partitioned physical- and biochemical-induced oxygen sinks and calculated the OCR and timescale for hypoxia regeneration after its destruction by a typhoon. The impacts of tropical cyclones on the evolution of seasonal hypoxia in river-dominated ocean margins are further discussed.

2 Materials and methods

2.1 Study area and cruise background

The shelf of the northern South China Sea (NSCS) receives an average annual freshwater discharge of ∼ 10 000 m³ s⁻¹ originating from the Pearl River, the 17th largest river in the world (Cai et al., 2004; Dai et al., 2014). Nearly four-fifths of freshwater discharge occurs during the wet season, typically from April to September (Dai et al., 2014). The riverine freshwater extends offshore to form a widespread plume over the shelf in summer (Gan et al., 2009b; Cao et al., 2011; Chen et al., 2017a) via eight outlets through three sub-estuaries (i.e. Lingdingyang, Modaomen and Huangmaohai; Fig. 1b). On the inner shelf, coastal upwelling interacts with the buoyant plume, which is propelled by the prevailing southwest monsoon and intensified along the eastward widened shelf (Gan et al., 2009a; Chen et al., 2017b). Climatologically, about seven tropical cyclones per year impacted the NSCS from 1949–2019, half of which featured maximum wind speeds greater than 32.7 m s⁻¹.

Figure 1 (a) Map of the study area on the shelf of the northern South China Sea (NSCS), showing the track of typhoon Son-Tinh (circles) across the NSCS during 16–24 July 2018. The colour of the circles represents the magnitude of wind speed. Additionally, the smaller circles denote tropical depression (wind speeds ≤ 17.1 m s⁻¹) and the larger circles denote tropical storm (wind speeds within 17.2–32.6 m s⁻¹). The arrows denote the locations of the typhoon as marked with time and wind speed. The grey lines are the depth contours at 50 and 200 m. (b) Sampling stations on the NSCS shelf off the Pearl River estuary in summer 2018. The pink, green, purple and orange circles denote the stations surveyed in all three legs, only leg 1 and leg 2, only leg 1, and only leg 2, respectively. Time-series observations were conducted at station F303 as marked by the star, and vertically high-resolution samplings were conducted at stations marked with bold circles. (c) The wind speed and (d) wind direction at Waglan Island (triangle in b) from May to August 2018. Bars at the bottom of (d) mark times when tropical cyclones impacted the NSCS. (e) The tidal height at the Dawanshan gauge station near station F303 from May to August 2018. The shaded area indicates the cruise periods for leg 1 (grey), leg 2 (pink) and leg 3 (blue).

Field observations and sampling were conducted aboard the R/V Haike 68 off the PRE on the inner-shelf of the NSCS in the summer of 2018. The cruise was interrupted by the passage of typhoon Son-Tinh across the NSCS, ∼ 350 km south of the PRE (Fig. 1a). Leg 1 (8–14 July) was the cruise period before the typhoon, and leg 2 (21–25 July) and leg 3 (26–29 July) were conducted after its passage (Fig. 1b). During each leg, we collected samples from west to east and along the cross-shelf transects within isobaths of 10–35 m. Almost all stations in leg 1 (56 stations) were revisited during leg 2 (56 stations, including 4 stations differing from leg 1) and nearly half again during leg 3 (27 stations). Eight stations were additionally revisited on the way back to the port on 31 July. Time-series observations with a sampling interval of 1 h were conducted at station F303 for 26 h before leg 2, beginning at 16:00 UTC+8 on 19 July (Fig. 1b). In contrast to the typical southwesterly winds in the NSCS with average monthly wind speeds of < 6 m s⁻¹ from June to September (Su, 2004), easterly winds prevailed during the cruise period due to the typhoon, with the wind speeds increasing up to ∼ 13 m s⁻¹ (Fig. 1c) at Waglan Island to the east of the study area (Fig. 1b).

2.2 Sampling and analysis

Temperature and salinity were determined using a SBE 917 plus conductivity–temperature–depth (CTD) recorder (Sea-Bird Electronics, Inc.). Discrete samples were collected using 5 L free-flow water samplers mounted onto a rosette sampling assembly. Dissolved oxygen (DO), dissolved inorganic carbon (DIC), total alkalinity (TA) and chlorophyll a (Chl a) concentrations were measured at all stations with depth profiles from 1 m below the surface down to ∼ 4–6 m above the bottom, generally at three depth layers. Additional high-resolution vertical samplings were conducted at seven to eight depth layers (Fig. 1b).

Salinity was calibrated against discrete water samples measured by a Multi 340i salinometer (WTW). The DO concentrations were measured onboard within ∼ 12 h using the spectrophotometric Winkler method (Labasque et al., 2004), with a precision better than ± 2 µmol L⁻¹. DIC was measured on ∼ 0.5 mL acidified water samples using an infrared CO₂ detector (Apollo ASC-3) with a precision of ± 2 µmol L⁻¹ (Cai et al., 2004). TA was determined on 25 mL samples in an open-cell setting based on the Gran titration technique (see details in Cai et al., 2010) with a Kloehn digital syringe pump. The analytical precision was ± 2 µmol L⁻¹. Both DIC and TA concentrations were calibrated against certified reference materials provided by Andrew G. Dickson at the Scripps Institution of Oceanography, University of California, San Diego, USA. Chl a concentrations were determined using a Trilogy laboratory fluorometer (Turner Designs, Inc.) after being extracted with 90 % acetone for 14 h at −20 ^∘C (Welschmeyer, 1994) and calibrated using a Sigma Chl a standard.

2.3 Water column stability

Water column stability regulates the ventilation of subsurface waters and replenishment of DO by suppressing turbulent mixing with stratification (Obenour et al., 2012; Lu et al., 2018; Cui et al., 2019), which can be indicated by the buoyancy frequency (also known as the Brunt–Väisälä frequency),

$$\begin{array}{}\text{(1)}& N^{\mathrm{2}}=-(g/\mathit{\rho })(\partial \mathit{\rho }/\partial z),\end{array}$$

where g is the gravitational acceleration, ρ is potential density and z is the height above the seabed. Generally, a positive N² (i.e. N²>0 s⁻²) indicates a stable regime where stratification may suppress turbulence (Tedford et al., 2009), and a larger N² value indicates a more stable water column. CTD temperature and salinity were used to calculate the buoyancy frequency for sampling stations, to indicate the controls of density stratification on the development of hypoxia (Sect. 4.1).

2.4 Oxygen consumption rate

From the perspective of Euler observations and based on mass balance, the DO changes (ΔDO) in subsurface waters over a specified time interval at a specific site can be decomposed into two components: one driven by physical mixing (ΔDO^mix) and the other induced by biochemical processes (ΔDO^bc). Here, we define the biochemical-induced DO consumption with time as the oxygen consumption rate (OCR). A higher OCR value indicates stronger oxygen consumption, and a negative value indicates oxygen production via biochemical processes (e.g. photosynthesis). For revisited stations, ΔDO is the difference in DO values measured between the two sampling periods. The physical-mixing-induced DO variations were derived using a three-endmember mixing model, which constructed the conservative mixing scheme among different water masses: brackish plume water (PW), offshore surface water (SW) and upwelled subsurface water (SUB) (Su et al., 2017; Cui et al., 2019; Zhao et al., 2020). The model is constrained by salinity (S) and potential temperature (θ) according to the following equations:

$$\begin{array}{}\text{(2)}& {\displaystyle }{f}_{\mathrm{PW}}+f_{\mathrm{SW}}+f_{\mathrm{SUB}}=\mathrm{1},\text{(3)}& {\displaystyle }{S}_{\mathrm{PW}}\times f_{\mathrm{PW}}+S_{\mathrm{SW}}\times f_{\mathrm{SW}}+S_{\mathrm{SUB}}\times f_{\mathrm{SUB}}=S^{\mathrm{meas}},\text{(4)}& {\displaystyle }{\mathit{\theta }}_{\mathrm{PW}}\times f_{\mathrm{PW}}+{\mathit{\theta }}_{\mathrm{SW}}\times f_{\mathrm{SW}}+{\mathit{\theta }}_{\mathrm{SUB}}\times f_{\mathrm{SUB}}={\mathit{\theta }}^{\mathrm{meas}}.\end{array}$$

where the superscript “meas” denotes measured values, and f represents the fraction that each endmember contributes to the in situ samples. Assuming that DO concentrations in surface waters before sinking to the depth were equilibrated with the atmosphere and that the subsurface waters were isolated from the atmosphere due to restriction by stratification, these fractions were applied to predict conservative concentrations of DO (DO^mix) resulting solely from conservative mixing.

$$\begin{array}{}\text{(5)}& {\mathrm{DO}}^{\mathrm{mix}}={\mathrm{DO}}_{\mathrm{PW}}\times f_{\mathrm{PW}}+{\mathrm{DO}}_{\mathrm{SW}}\times f_{\mathrm{SW}}+{\mathrm{DO}}_{\mathrm{SUB}}\times f_{\mathrm{SUB}}\end{array}$$

Similarly, ΔDO^mix is the difference in “conservative” DO values between visits, assuming that the bottom water masses where biochemical oxygen consumption prevailed were constrained by strong convergence and that their outflow from the sampling area is insignificant on the timescale of the water residence time (Lu et al., 2018; Li et al., 2020). As a result,

$$\begin{array}{}\text{(6)}& \mathrm{OCR}=-\mathrm{\Delta }{\mathrm{DO}}^{\mathrm{bc}}/\mathrm{\Delta }t=-\left(\mathrm{\Delta }\mathrm{DO}-\mathrm{\Delta }{\mathrm{DO}}^{\mathrm{mix}}\right)/\mathrm{\Delta }t,\end{array}$$

where Δt is the duration between the two observations at times t₁ and t₂ (t₂>t₁), and ΔX = $X^{t_{\mathrm{2}}}-X^{t_{\mathrm{1}}}$ (X = DO, DO^mix, DO^bc, etc.). The uncertainty in the calculation of OCR mainly derives from the estimation of conservative values predicted from the three-endmember mixing model. Sources of the composite uncertainty (ε) in derivation of DO^mix are associated with potential temperature (θ), salinity (S) and the dissolved oxygen (DO) values of endmembers.

$$\begin{array}{}\text{(7)}& {\mathit{\epsilon }}_{{\mathrm{DO}}^{\mathrm{mix}}}=\sqrt{{\sum }_i^n\left[{\left(f_i\cdot {\mathit{\sigma }}_{{\mathrm{DO}}_i}\right)}^{\mathrm{2}}+{\left({\mathrm{DO}}_i\cdot {\mathit{\sigma }}_{f_i}\right)}^{\mathrm{2}}\right]},\end{array}$$

where σ_DO and σ_f are uncertainties in the DO concentration and the fraction of each endmember i (i.e. PW, SW and SUB), respectively, and the latter of which can be calculated as

$$\begin{array}{}\text{(8)}& {\mathit{\sigma }}_{f_i}=\sqrt{{\sum }_j^n\left[{\left(\partial f_i/\partial {\mathit{\theta }}_j\cdot {\mathit{\sigma }}_{{\mathit{\theta }}_j}\right)}^{\mathrm{2}}+{\left(\partial f_i/\partial S_j\cdot {\mathit{\sigma }}_{S_j}\right)}^{\mathrm{2}}\right]},\end{array}$$

where j also denotes each endmember.

2.5 Endmember selection and validations of the three-endmember mixing model

The potential-temperature–salinity diagram is shown in Fig. 2a. The three-endmember mixing scheme for the bottom layer has been elucidated in Zhao et al. (2020) that had a spatial coverage similar to this study. We adopted the endmember values of the offshore surface water and upwelled subsurface water from Zhao et al. (2020) given that our sampling was almost exclusively within the 30 m isobaths (Fig. 1b), and these values were consistent with those found in previous studies (Cao et al., 2011; Guo and Wong, 2015; Su et al., 2017). The brackish plume water was assumed to partly subduct to the bottom layer under downwelling-favourable wind conditions (Huang et al., 2019; Li et al., 2021). The endmember values of the brackish plume water were thus determined from the surface water samples near the mouth of the PRE with a salinity of ∼ 16.9, mainly consisting of a mixture of riverine freshwater and offshore surface water. The DIC endmember of brackish plume water here was consistent with the predicted value using the endmember values of riverine/plume water reported by Su et al. (2017), whereas it was higher than that calculated using the endmember values of riverine freshwater from Zhao et al. (2020) since the riverine DIC concentrations might be diluted by abnormally high river discharge in 2017 (Guo et al., 2008). For simplification, DO concentrations in offshore surface water were assumed to be saturated, in equilibrium with the atmosphere, while the upwelled subsurface water was assumed to be oxygen deficient by ∼ 16 % relative to the saturation level. The DO endmember value of brackish plume water was also assumed to be equilibrated with the atmosphere, which should be in order, because the biological productivity was largely limited by high turbidity in shallow estuarine waters. A summary of the endmember values is listed in Table 1. In estimating the OCR, we excluded the above-pycnocline samples collected at depths < 10 m affected by the upper plume waters that are subject to strong air–sea exchanges and/or photosynthetic production of oxygen.

Figure 2 (a) Potential temperature (^∘C) vs. salinity, (b) predicted TA (TA^pre, µmol kg⁻¹) vs. measured TA (TA^meas, µmol kg⁻¹) and (c) ΔDIC^m−p (µmol kg⁻¹) vs. ΔDO^m−p (µmol kg⁻¹) on the NSCS shelf off the PRE. The black-edged circles represent bottom water samples with depths >10 m. The yellow, green and purple triangles in (a) represent the endmember values of brackish plume water (PW), offshore surface water (SW) and upwelled subsurface water (SUB), respectively. The black line in (c) denotes the slope of ΔDIC^m−p plotted against ΔDO^m−p derived from the Model II regression.

Table 1 Summary of the endmember values adopted in the three-endmember mixing model.

^a Adopted from Zhao et al. (2020). ^b Uncertainties were derived from multiple samples collected at the entrance of the PRE. ^c Uncertainties were calculated by propagating errors associated with the estimation of oxygen solubility using Benson and Krause (1984).

The predicted quasi-conservative TA (TA^pre= TA${}_{\mathrm{PW}}\times f_{\mathrm{RW}}+$ TA${}_{\mathrm{SW}}\times f_{\mathrm{SW}}+$ TA_SUB × f_SUB; same for DIC^pre) is mostly consistent with our measured values (Fig. 2b), with a subtle difference of 8 ± 8 µmol kg⁻¹ likely associated with measurement errors, computational errors in the mixing scheme and/or biological processes. The slope of ΔDIC (ΔDIC^m−p = DIC^meas − DIC^pre) vs. ΔDO^m−p (ΔDO^m−p = DO^meas − DO^mix) in bottom waters was −0.93 ± 0.07 (Fig. 2c), similar to that reported by Zhao et al. (2020).

3 Evolution of intermittent hypoxia off the PRE

3.1 Extensive hypoxia before the typhoon

The water columns showed a prominent two-layer structure on the inner NSCS shelf off the PRE (Figs. 3 and 4). In the surface layer, the freshwater plume mainly attached to the coast when veering to the west as constrained by estuarine topography, the Coriolis force (Wong et al., 2003) and easterly winds during leg 1 (Fig. 1d), despite a freshwater bulge that remained near the mouth of the Lingdingyang sub-estuary due to persistent southwesterly winds before the cruise (Fig. 1d) and the weak shelf current there (Gan et al., 2009b; Lu et al., 2018). A strong bloom occurred in the surface plume waters near the Huangmaohai sub-estuary, characterized by high Chl a concentrations of >20 µg L⁻¹ and oversaturated DO of >300 µmol kg⁻¹ (equivalent to a DO saturation level >150 %) (Fig. 3g, j). The freshwater bulge also featured a relatively weak bloom, with Chl a concentrations of ∼ 10 µg L⁻¹ and DO of ∼ 250 µmol kg⁻¹ (equivalent to a DO saturation level of ∼ 125 %). Similar blooms were previously observed to be associated with high nutrient concentrations, a sufficiently long water residence time and an abundance of photosynthetically active radiation (Lu and Gan, 2015). The surface water temperature was mostly >29 ^∘C. Exceptions occurred to the southwest of Hong Kong (Fig. 3a) where the air temperature was 2–3 ^∘C lower at the time of sampling (Fig. S1 in the Supplement).

Figure 3 Distributions of temperature (^∘C), salinity, DO (µmol kg⁻¹) and Chl a concentrations (µg L⁻¹) at the surface water layer off the PRE during pre-typhoon leg 1 and during post-typhoon legs 2 and 3. The white and magenta contours in (g) and (w) show the hypoxic (DO < 63 µmol kg⁻¹) and oxygen-deficit (DO < 94 µmol kg⁻¹) zones, respectively. Figures were produced using Ocean Data View v. 5.3.0 (http://odv.awi.de, last access: 8 June 2020, Schlitzer, 2018).

At the bottom layer, low-temperature (< 26 ^∘C), high-salinity (>33) deep shelf waters intruded upslope to the 10–20 m isobaths beneath the surface plume (Fig. 4a, d). An extensive hypoxic zone (DO < 63 µmol kg⁻¹) developed beneath the freshwater bulge and extended westwards along the 20–30 m isobaths to the region off the Modaomen sub-estuary (Fig. 4g). To the east, a relatively weak hypoxic centre occurred adjoining the Hong Kong waters. Additionally, a smaller-scale hypoxic zone appeared beneath the surface bloom near the Huangmaohai sub-estuary, a region where hypoxia was also reported but not fully covered by survey measurements (Su et al., 2017; Cui et al., 2019). The general spatial pattern of hypoxic centres was similar to that found in summers of 2014 and 2017 (Su et al., 2017; Zhao et al., 2020) yet with a slight offshore shift. The minimum oxygen level was 3.5 µmol kg⁻¹ (∼ 0.1 mg L⁻¹) in bottom waters at station F303, which was lower than the previously reported minimum (∼ 7 µmol kg⁻¹ at station F304; Su et al., 2017). Within the surveyed region where DO concentrations were interpolated into ${\mathrm{0.5}}^{\prime }\times {\mathrm{0.5}}^{\prime }$ grids using the Kriging interpolation method, the total area of the hypoxic zone reached ∼ 660 km², and the oxygen-deficit zone (DO < 94 µmol kg⁻¹) occupied ∼ 1470 km², which was larger than those in summer 2014 (>280 km² for the hypoxic zone and >800 km² for the oxygen-deficit zone) when we surveyed in a smaller area (Su et al., 2017). Cui et al. (2019) reported that spring–neap tidal oscillations lead to variations in the DO concentration outside of the Lingdingyang sub-estuary with a maximum neighbouring oxygen range of ∼ 15 µmol kg⁻¹. As our observations were conducted during the transition from the neap tide to the spring tide (Fig. 1e), the total areas of the hypoxic zone and the oxygen-deficit zone would be at most ∼ 990 and ∼ 1930 km², respectively, which are 34 %–50 % larger than our observed areas.

Figure 4 Distributions of temperature (^∘C), salinity, DO (µmol kg⁻¹) and Chl a concentrations (µg L⁻¹) at the bottom water layer off the PRE during pre-typhoon leg 1 and during post-typhoon legs 2 and 3. The white and magenta contours in (g) and (w) show the hypoxic (DO < 63 µmol kg⁻¹) and oxygen-deficit (DO < 94 µmol kg⁻¹) zones, respectively.

3.2 Destruction of hypoxia by the typhoon

The spatial patterns of temperature, salinity, DO and Chl a concentrations all changed drastically from disturbance by intensified easterly winds during the typhoon period. Strong winds drove the warm, low-salinity and oxygen-saturated surface waters to mix downward, increasing temperature and DO concentrations and decreasing salinity in bottom waters (Fig. 4). Time-series observations at station F303 before leg 2 showed a vertically well-mixed water column in the first half, as reflected by homogeneous distributions of intermediate temperature (∼ 28 ^∘C), salinity (∼ 32) and DO levels (∼ 180 µmol kg⁻¹) (Fig. 5a–c). The less well-mixed water column in the second half of the time-series observations likely resulted from the weakened winds, showing an upward intrusion of bottom waters which were slightly warmer than surface waters which lost heat to the low-temperature atmosphere (Fig. S1). With subdued winds that shifted to southwesterly in the following 2 d (Fig. 1c, d), the offshore spreading of the river plume suppressed the upward intrusion of slightly warm bottom waters and facilitated the restoration of a two-layer water column, as observed in the second half of the time-series observations (Fig. 5) and during leg 2 (Figs. 3 and 4). Large phytoplankton blooms were identified in the surface plume, widely spreading from the mouth of the Lingdingyang sub-estuary to near the Huangmaohai sub-estuary (Fig. 3k), potentially fuelled by nutrients mixed upward from the depth in addition to riverine inputs (Wang et al., 2017b; Qiu et al., 2019). The maximum Chl a concentration was >40 µg L⁻¹ off the Modaomen sub-estuary, accompanied by an extraordinarily high DO concentration of >350 µmol kg⁻¹ (Fig. 3h, k). The surface water temperature was ∼ 29.0 ± 0.5 ^∘C, which was higher than that in bottom waters by 0.8 ^∘C during leg 2 (Figs. 3b and 4b).

Figure 5 Time-series observations of (a) temperature (^∘C), (b) salinity, (c) DO (µmol kg⁻¹) and (d) Chl a concentrations (µg L⁻¹) at station F303 (see Fig. 1b) from 19–20 July 2018 after the typhoon passage, showing the complete destruction and the subsequent rapid development of stratification.

In bottom waters, hypoxia had been completely destroyed due to strong reaeration in the wake of the typhoon travelling across the NSCS, replaced by a homogenous spatial distribution of relatively high DO concentrations (∼ 171 ± 16 µmol kg⁻¹) during leg 2 (Fig. 4h). The cross-shore gradients in temperature and salinity were also largely relaxed, with isotropically elevated temperatures up to ∼ 28 ^∘C and decreased salinity (Fig. 4b, e). The mid-depth distributions of temperature, salinity and DO concentrations showed similar patterns as at the bottom layer (Fig. S2). Although the water column remained relatively well-mixed in the subsurface layer, freshwater buoyancy and weakened winds facilitated the revitalization of density stratification and subsequent oxygen decline below the pycnocline. Indeed, the bottom water DO concentration at station F303 decreased by ∼ 18 µmol kg⁻¹ compared to that in the time-series observations and was lower than that at adjacent stations by ∼ 9–22 µmol kg⁻¹ when revisited during leg 2 on 22 July (Fig. S3).

3.3 Reinstatement of hypoxia after the typhoon

With the dying out of the typhoon after its landfall to the west of the study area on 23 July (Fig. 1a), the wind speed decreased to < 5 m s⁻¹ on 25 July, while the wind direction remained from the southeast before it shifted to southwesterly on 29 July (Fig. 1c, d). The DO concentrations in bottom waters were noticeably lower (139–164 µmol kg⁻¹ shallower than 20 m isobaths), starting around July 23 in the second half of leg 2 (Figs. 1b and 4h). During leg 3 from 26–29 July, hypoxia was present again in bottom waters, due to favourable conditions for its formation (Figs. 3 and 4). The surface layer was warmed up to over 30 ^∘C (Fig. 3c), increasing vertical thermal gradients relative to leg 2 (Figs. 3 and 4) and strengthening the stratification (Allahdadi and Li, 2017). The freshwater bulge of lower salinity (< 15) advected offshore around the Modaomen sub-estuary, with the offshore migration of the bloom (Fig. 4f, l) likely driven by the interaction between the seaward buoyant current and northeastward shelf current (Pan et al., 2014; Li et al., 2020). The Chl a concentrations near the entrances of the three sub-estuaries remained relatively high (>10 µg L⁻¹), and the DO concentrations remained at high levels of >250 µmol kg⁻¹, which were ∼ 20 % over the saturation levels (Fig. 3i, l).

Similar to leg 1 and leg 2, the surface waters penetrated into the subsurface layer along the coast (Fig. 4f), likely forced by the downwelling-favourable winds (Fig. 1d) (Huang et al., 2019; Li et al., 2021), augmenting temperature and DO concentrations but bringing down salinity, particularly in the mid-depth layer from 3–5 m at stations nearest the shore to ∼ 15 m at stations furthest from the shore (Fig. S2). The downward penetration of surface waters, nonetheless, seemed to be restricted to the ∼ 10 m isobath and thus offset only a limited amount of the oxygen reduction caused by biochemical consumption (Koweek et al., 2020). DO concentrations in bottom waters were reduced to ∼ 46 µmol kg⁻¹ off the Modaomen sub-estuary along the 20 m isobaths and to < 94 µmol kg⁻¹ to the southwest of Hong Kong (Fig. 4i), indicating the reinstatement of hypoxia. When sites were revisited on 31 July, the re-emerging hypoxia was found to have been strengthened, with oxygen levels down to ∼ 37 µmol kg⁻¹. We found that hypoxia formed in bottom waters with a temperature of ∼ 28 ^∘C off the Modaomen sub-estuary during leg 3, while the oxygen-deficit zone to the southwest of Hong Kong showed a relatively low temperature of < 27 ^∘C (Fig. 4c), likely due to the cross-isobath transport of deep shelf waters, arising from local topographic effects (Dai et al., 2014; Wang et al., 2014) as during leg 1. In this sense, the shoreward intrusion of deep shelf waters is not a prerequisite for the initiation of hypoxia formation off the PRE, but it contributes to the reinstatement of hypoxia southwest off Hong Kong.

4 Physical and biogeochemical controls on the evolution of intermittent hypoxia

4.1 Water column stability

A stable water column is a key prerequisite for the formation and maintenance of hypoxia in coastal oceans (Wang and Justić, 2009; Obenour et al., 2012; Testa and Kemp, 2014; Lu et al., 2018; Zhang et al., 2019), which restricts the oxygen supply by suppressing advective and diffusive mixing with oxygen-rich waters (Murphy et al., 2011; Cui et al., 2019). Many studies have demonstrated that density stratification becomes enhanced and stabilizes subsurface waters when freshwater flows over seawater (Gan et al., 2009b; MacCready et al., 2009; Bianchi et al., 2010), allowing oxygen depletion over a longer timescale (Fennel and Testa, 2019). In the second half of the time-series observations, the two-layer structure re-emerged with the spreading of the river plume, which suppressed the vertical mixing in the subsurface layer, as reflected by the receded upward intrusion of slightly warm, saline bottom waters (Fig. 5a, b). The Chl a concentrations in the surface layer increased, followed by an elevated DO level in surface waters and a lowered DO level in subsurface waters (Fig. 5c, d). The surface-to-bottom salinity difference showed large values within the surface plume area, which almost covered the bottom hypoxic zones (Fig. 4g, i). Exceptions only occurred to the hypoxic zone off the Modaomen sub-estuary, where the surface-to-bottom salinity differences were relatively small but the temperature differences were large (i.e. ΔT_b−s < −4 ^∘C in Fig. 6a–c) due to the shoreward intrusion of cold deep shelf waters (Fig. 4a–c). The regions occupied by the surface plume and the shoreward-intruded deep shelf waters therefore overlapped, resulting in a more stable water column where a patchy hypoxic zone could persist for more than 5 d (Cui et al., 2019).

Figure 6 The surface-to-bottom temperature (a–c) and salinity (d–e) distributions off the PRE during pre-typhoon leg 1 and post-typhoon during legs 2 and 3. ΔT_b−s and ΔS_b−s represent the difference in temperature and salinity, respectively, between the bottom and surface layer.

Off the PRE, when not influenced by freshwater inputs, the surface layer showed a relatively small N² close to 0 s⁻² and vertically well-mixed temperature, salinity and DO concentrations (e.g. in the top ∼ 10 m at station A14, Fig. 7c). However, in the presence of the freshwater plume, the surface layer became more stable, with a larger N² of >1 × 10⁻³ s⁻², or even >5 × 10⁻³ s⁻² (e.g. stations A8 and A11, Fig. 7a, b). The DO concentrations decreased sharply at the base of the surface plume (salinity ∼ 30). Using the eddy diffusivity for density (K_z) of < 5 × 10⁻⁶ m² s⁻¹ for N² larger than 1 × 10⁻³ s⁻² (Cui et al., 2019), we estimated the vertical diffusion for DO concentrations (VDIF = K_z × (∂DO $/$ ∂z)) of ∼ 0.25 g m⁻² d⁻¹ with a maximum of 0.54 g m⁻² d⁻¹ in the top 10 m at stations A8 and A11, which was comparable to the results from Cui et al. (2019). It therefore acted as a barrier layer, with weak dissipation of oxygen into the subsurface waters. The inherent pycnocline between the offshore surface water and deep shelf waters mainly driven by steep temperature gradients (Qu et al., 2007), such as at station A14, acted as a second pycnocline in the plume region (e.g. stations A8 and A11, Fig. 7a) yet with weaker stratification likely from increased shear stresses in shallower waters (Pan and Gu, 2016; Li et al., 2020). At the edge of the freshwater plume (Figs. 1b and 3a), such as station F304, the second pycnocline was more stable than that induced by freshwater inputs (Fig. 7d). This three-layer structure, separated by two pycnoclines, therefore effectively decreased oxygen influx from the surface and facilitated oxygen depletion in bottom waters.

Figure 7 Profiles of temperature (^∘C) (blue dashed lines), salinity (pink solid lines), dissolved oxygen (DO, µmol kg⁻¹) (purple dots) and buoyancy frequency N² (s⁻²) (bold black solid lines) at stations A8, A11, A14 and F304 (see Fig. 1b), with visits for both pre-typhoon (leg 1) and post-typhoon periods (legs 2 and 3). The vertical distributions of N² have been smoothed by the Gaussian method.

Water column stability also largely depends on wind stress and/or direction in coastal waters (Wilson et al., 2008; Wang and Justić, 2009). Higher wind stress usually de-stratifies the water column, leading to stronger turbulent mixing, air–sea gas exchange and reaeration (Chen et al., 2015; Lu et al., 2018; Huang et al., 2019), relieving hypoxic conditions (Ni et al., 2016; Wei et al., 2016). During the typhoon period, the wind speed rose to as high as 13 m s⁻¹ (Fig. 1c), which was large enough to break the stratification (Geng et al., 2019) driven by freshwater inputs and the inherent thermocline and resulted in a vertically well-mixed water column (Fig. 5). The surface waters became undersaturated (∼ 90 % of the oxygen saturation level) due to the upward mixing of low-oxygen waters, which in turn favoured the ventilation of bottom waters and the breakdown of hypoxic conditions (Hu et al., 2017). Despite wind speeds still as high as 10 m s⁻¹ during leg 2 (Fig. 1c), stratification was regenerated in the top ∼ 10 m of the plume region, which had relatively large surface-to-bottom salinity differences (Fig. 6e) and high N² (Fig. 7e–h) under easterly winds (Fig. 1d). This suggests that freshwater input-induced stratification suppressed turbulent mixing driven by wind stress, favouring the initiation of hypoxia development even under downwelling-favourable conditions. Most modelling work simulated the formation of hypoxia off the PRE under the typical summer southwest monsoon (Wei et al., 2016; Lu et al., 2018; Li et al., 2020; Yu et al., 2020); when southwest winds blow in a more southward direction, a larger hypoxic zone develops (Wei et al., 2016). In contrast to southwesterly winds that facilitate the wide eastward spreading of the surface plume (Gan et al., 2009b), the downwelling-favourable easterly winds tend to constrain the surface plume to flow westward near the coast (Li et al., 2021) as shown during leg 2 and even drive the surface oxygen-saturated waters to penetrate into the deep along the coast (Figs. S2 and 4), resulting in an offshore or westward shift of the hypoxic zones that only partly overlapped with the surface plume in terms of their localities. If the easterly winds last for a longer time than the hypoxia formation timescale, stronger blooms in the surface plume (Fig. 3k) would enhance the bottom hypoxia with abundant supply of fresh, labile organic matter, but the downwelling-favourable winds would also destroy the bottom hypoxia if the wind stresses become strong enough (Li et al., 2021).

4.2 Oxygen sinks and hypoxia formation timescale

4.2.1 Mixing-induced oxygen sinks

With the restoration of density stratification from leg 2 to leg 3, DO concentrations in bottom waters were generally reduced by >25 µmol kg⁻¹, with two hotspots showing reductions up to 75 µmol kg⁻¹: one located offshore of the Lingdingyang sub-estuary, to the southwest of Hong Kong, and the other between the Modaomen and Huangmaohai sub-estuaries (Fig. 8a). The mixing-induced DO changes were positive along the coast within 10–20 m isobaths, as the oxygen-saturated surface waters penetrated downward to re-aerate the bottom waters driven by the downwelling-favourable easterly winds (Huang et al., 2019). The mixing-induced oxygen sinks mainly occurred in bottom waters southwest off Hong Kong, with an average of −5.7 ± 0.8 µmol kg⁻¹, which was higher than in other regions west of the PRE (e.g. beyond the 20 m isobath; −1.4 ± 0.8 µmol kg⁻¹) (Fig. 8b). The mixing-induced oxygen sinks can be attributed to the shoreward intrusion of oceanic, cold oxygen-undersaturated subsurface waters – as reflected by lowered temperature (< 27 ^∘C; Fig. 4c), which usually act as a non-local driver of coastal hypoxia by lowering the initial DO concentration (Wang, 2009; Qian et al., 2017). The cold, saline oceanic subsurface waters also completely occupied the bottom layer beyond the 20 m isobath during leg 1, where extensive hypoxia developed (Fig. 4a, g). This upwelling-induced reduction in the initial DO level amounted to 8.6 ± 1.7 % of the oxygen decline to the southwest of Hong Kong, suggesting coastal upwelling played a minor role in the hypoxia formation.

Figure 8 Distributions of (a) total DO changes (ΔDO, µmol kg⁻¹), (b) mixing-induced DO changes (ΔDO^mix, µmol kg⁻¹) and (c) the biochemical-induced oxygen consumption rate (OCR, µmol O₂ kg⁻¹ d⁻¹) between leg 3 and leg 2 on the inner NSCS shelf off the PRE. (d) The biochemical-induced changes in DIC (ΔDIC^bc, µmol kg⁻¹) vs. DO (ΔDO^bc, µmol kg⁻¹) in bottom waters with depths >10 m from leg 2 to leg 3. The black line denotes the slope of ΔDIC^bc plotted against ΔDO^bc derived from the Model II regression.

The contribution of oxygen-deficient coastal upwelling to the hypoxia formation varies in different ocean marginal systems, which is largely dependent on the source of the subsurface water masses and biogeochemical reactions along the pathway of the intrusion. Qian et al. (2017) showed the apparent oxygen utilization (AOU) values of ∼ 50 µmol kg⁻¹ in the source waters of the nearshore Kuroshio branch at the shelf break northeast of Taiwan and an AOU increment of ∼ 40 µmol kg⁻¹ throughout its travel time of ∼ 60 d to the vicinity off the Changjiang river (Yang et al., 2013). In eastern boundary upwelling systems along the northeast Pacific Ocean, the source waters typically have a low DO concentration of ∼ 80–160 µmol kg⁻¹, even near or below hypoxic levels from a depth of ∼ 100–200 m (Grantham et al., 2004). The DO deficits were further exacerbated by respiration when the source waters transited shoreward over the shelf. Comparing with these systems, the source water originating from the low-latitude high-temperature oligotrophic NSCS (Wong et al., 2007) has a higher DO level (AOU ∼ 35 µmol kg⁻¹). As the NSCS shelf is narrower than the East China Sea shelf, the source water spends a shorter time intruding towards the coastal zone after diverting its direction from the continental slope (Gan et al., 2009a; Wang et al., 2014). These might explain the relatively low contribution of coastal upwelling to the oxygen depletion on the inner NSCS shelf off the PRE.

4.2.2 Biochemical-induced oxygen sinks

The biochemical-induced DO and DIC changes from leg 2 to leg 3 also showed a good relationship with a slope of −0.92 ± 0.17 (Fig. 8d), consistent with the slope of biochemical-induced variations in DO against DIC concentrations throughout the sampling legs (Fig. 2c), implying that aerobic respiration of organic matter indeed dominated the oxygen consumption. The distribution of OCR estimates in bottom waters almost mirrored variations in the total DO pattern from leg 2 to leg 3 (Fig. 8a, c). This biochemically mediated OCR ranged from 0.9 µmol O₂ kg⁻¹ d⁻¹ at offshore non-hypoxic stations to 19.5 ± 0.4 µmol O₂ kg⁻¹ d⁻¹ in hypoxic waters, with an average of 14.6 ± 4.8 µmol O₂ kg⁻¹ d⁻¹ in the oxygen-deficit zone. The uncertainty introduced by the mixing scheme was 0.63–0.98 µmol O₂ kg⁻¹ d⁻¹, accounting for a deviation of 4 %–27 %. The geographical locations of high-OCR zones at the bottom layer (OCR >10 µmol O₂ kg⁻¹ d⁻¹) were not fully overlapped with those of surface blooms (Figs. 3l and 8c). Cui et al. (2019) suggested that the patchy distribution of bottom hypoxia was closely associated with the river plume front which traps organic particles (Hetland and DiMarco, 2008), accelerating their settlement and deposition in the overlapping zones between the river plume and shelf salt wedge (Zhang and Li, 2010) to fuel the high OCR there. In fact, only the west hypoxic centre was located beneath the surface bloom from leg 2 to leg 3 (Figs. 3j–l and 4g–i). Lu et al. (2018) proposed another explanation that small detritus of external sources could be accumulated and remineralized in the bottom flow convergence zone as regulated by highly variable coastline and bottom topography. Our observed hypoxic zones generally coincided with the modelled strong convergence zones (Lu et al., 2018; Li et al., 2020), implying sufficient organic matter supply fuelled the high OCR to restore hypoxia in bottom waters.

Our estimated OCR is comparable in magnitude to the community/bacterial respiration rate from previous studies in this study area (9.6 µmol O₂ kg⁻¹ d⁻¹; Su et al., 2017; 7.9 to 19.0 µmol O₂ kg⁻¹ d⁻¹; Cui et al., 2019; 16.8 ± 8.9 µmol O₂ kg⁻¹ d⁻¹; Li et al., 2019) and within the range in other estuaries and coastal systems (Dortch et al., 1994; Robinson, 2008), despite a lower limit due to a potential overestimation of the actual time if significant oxygen consumption started later than our observations during the first half of leg 2 and an assumed negligible oxygen flux supplied from the surface by diffusion. From leg 2 to leg 3, the vertical oxygen diffusion flux was ∼ 0.18 g m⁻² d⁻¹, amounting to a supply of ∼ 0.6–1.1 µmol O₂ kg⁻¹ d⁻¹ to the bottom waters with a thickness of 5–10 m. The exclusion of this vertical oxygen diffusion flux in our estimates therefore might underestimate the OCR by ∼ 6 % on average. We also assumed that the OCR was uniform in the subsurface waters several metres above the seabed, because we collected samples at three depth layers during leg 2 and leg 3, with usually only one layer below the pycnocline for most stations. For stations with the middle layer also below the pycnocline, such as station F304 during leg 3, the profile of DO concentrations showed almost constant values below the pycnocline (Fig. 7k). Sediment oxygen demand might be significant near the seabed or in its overlying water column in shallow waters (Kemp et al., 1992; Zhang and Li, 2010). However, our sampling at the bottom layer was 4–6 m above the seabed, and the water depth of the hypoxic zones was nearly 20 m (Fig. 4). In such sampling areas, Cui et al. (2019) found that oxygen losses by sediment oxygen demand (i.e. benthic respiration) were much smaller than the bacterial respiration in the water column based on both incubation experiments and oxygen budget analysis. We thus assumed the sediment oxygen demand was negligible and the microbial respiration in the water column dominated the estimated OCR.

However, the estimated OCR here differs in implications from that derived from the steady-state budget analysis in which the oxygen consumption was completely offset by the physical transport of oxygen (Zhang and Li, 2010; Cui et al., 2019). In our estimates, the OCR was derived during the hypoxia formation in which the in situ DO concentrations decreased over time (i.e. in a non-steady-state system) as oxygen sinks exceed sources. To shift from such a non-steady-state system towards a balance between oxygen consumption and replenishment for the maintenance of hypoxia, the OCR might decrease or the physical-induced oxygen supply increases. Kalvelage et al. (2015) found that aerobic respiration rates decreased with decreasing in situ DO concentrations, i.e. from the upper boundary of oxygen minimum zone (OMZ) towards the OMZ core off Peru. The oxygen enriched incubation of unfiltered water samples revealed that the OCR changed little when the in situ DO concentration was higher than ∼ 90 µmol kg⁻¹ but also decreased when the in situ DO concentration decreased to ∼ 30 µmol kg⁻¹ (He et al., 2014). Despite the bottom DO concentrations of ∼ 180 µmol kg⁻¹ at station F303 during the time-series observations, DO declined at a rate of ∼ 9 µmol O₂ kg⁻¹ d⁻¹ from 20–22 July, when the winds remained strong, and the OCR decreased to ∼ 5.5 µmol O₂ kg⁻¹ d⁻¹ from leg 2 (DO ∼ 164 µmol kg⁻¹) to leg 3 (DO ∼ 127 µmol kg⁻¹) (from 22–29 July) (Fig. S3). The OCR declined likely as inhibited by the reduced supply of labile organic matter (Yuan et al., 2010) – the strong blooms shifted offshore with the river plume, and the Chl a concentrations over the hypoxic zone also decreased from leg 2 to leg 3 (Fig. 3k, l). If we considered the shoreward intrusion of deep shelf waters following up the hypoxia development instead of contributing to lowering the initial DO level, their intrusion onto the inner shelf supplied oxygen to the oxygen-depleted waters, such as the cold, saline subsurface waters occupied the bottom layer covering the hypoxic zone on the inner shelf during leg 1 (Fig. 4a, g). As density stratification limited the vertical diffusivity of oxygen from the upper layer, bottom hypoxia was therefore progressively formed and maintained until the decreasing OCR almost achieved equilibrium with the oxygen supply from the lateral advection.

4.2.3 Hypoxia formation timescale

The above-estimated OCR and the contribution of shoreward-intruded deep shelf waters made it possible to roughly estimate the hypoxia formation timescale, i.e. the time that the DO level in a known volume of water takes to decrease below the hypoxia threshold from an assumed initial DO concentration (Fennel and Testa, 2019). In this study, for the reinstatement of hypoxia after the typhoon, the initial DO level in subsurface waters could be taken as ∼ 180 µmol kg⁻¹, which varied within ± 5 µmol kg⁻¹ throughout the time-series observations (Fig. 5c). Considering the OCR for the hypoxic zone was at most ∼ 20 µmol O₂ kg⁻¹ d⁻¹ and the negligible contribution of lateral advection to oxygen loss off the Modaomen sub-estuary (Fig. 4), it took nearly 6 d for the drawdown of DO to reach concentrations of ∼ 63 µmol kg⁻¹ within a limited area. Scaling to a larger area, it would instead take 8–12 d if we choose the average OCR of the oxygen-deficit zone, i.e. ∼ 15 ± 5 µmol O₂ kg⁻¹ d⁻¹ (Fig. 8c). These estimates of hypoxia formation timescales are shorter than the bottom water residence time (>15 d, Li et al., 2020), which favoured the hypoxia formation to the west off the PRE (Su et al., 2017; Zhao et al., 2020).

For a more common scenario of hypoxia formation starting from late spring, we assumed that the shoreward-intruded oxygen-deficient offshore subsurface waters almost occupied the bottom layer on the inner shelf instead of the well-oxygenated offshore surface waters, as during leg 1 (Fig. 4a, d). Based on the above estimates that the shoreward-intruded subsurface waters totally reduced the initial DO level by 8.6 ± 1.7 % of the oxygen decline for hypoxia formation, the initial DO level was lowered by ∼ 11 µmol kg⁻¹ to be ∼ 183 µmol kg⁻¹ before biochemical oxygen consumption. The hypoxia hotspot will then first occur ∼ 6 d after its initiation on the inner NSCS shelf off the PRE and extend to a larger hypoxic zone within the water residence time of >15 d (Lu et al., 2018; Li et al., 2020). This result is larger than that estimated by Fennel and Testa (2019) – 4 d – using the modelled OCR of ∼ 34 µmol O₂ kg⁻¹ d⁻¹ in the water column with a sediment oxygen demand of ∼ 2.1 g m⁻² d⁻¹, which were only applicable to the hypoxia formation in the Lingdingyang sub-estuary with shallower waters (∼ 5 m) (Zhang and Li, 2010). This result is still at the lower end of the hypoxia formation timescale in large river-dominated shelves globally (e.g. East China Sea off the Changjiang estuary, Northern Gulf of Mexico and northwestern Black Sea), which varies from 8 to 89 d for hypoxia to develop once initiated (Fennel and Testa, 2019). This short hypoxia formation timescale likely owes to a high OCR in relatively warm subsurface waters fuelled by abundant labile organic matter (Su et al., 2017; Zhao et al., 2020).

4.3 Imprint of tropical cyclones on the evolution of coastal hypoxia

Tropical cyclones dramatically alter the physical stability of the water column and attenuate or even disrupt hypoxic conditions at low and mid-latitudes, such as off the Changjiang Estuary (Ni et al., 2016; Zhang et al., 2020) and the PRE (Su et al., 2017; Huang et al., 2019), in Chesapeake Bay (Testa and Kemp, 2014; Testa et al., 2017), and in the northern Gulf of Mexico (Wang and Justić, 2009; Feng et al., 2012). These intense, episodic storms thus strongly impact the duration and intensity of oxygen depletion in coastal bottom waters (Rabalais et al., 2009; Wang et al., 2017a; Zhang et al., 2020), driving the seasonal hypoxia to be intermittent. As shown in Fig. 1, at least four named tropical cyclones impacted the study area from May to August in 2018, most of which shifted the wind direction from prevailing southwesterly to easterly or southeasterly and increased the wind speed up to over 9 m s⁻¹, being able to destroy water column stratification and interrupt hypoxia formation (Geng et al., 2019). This was not an exception, as annually there were approximately six tropical cyclones travelling across the NSCS from May to September, when seasonal hypoxia develops (Qian et al., 2018; Wang et al., 2018), during the period from 1975–2019 (Fig. 9a). On average, five out of the six tropical cyclones on an annual basis had the maximum wind speeds exceeding 9 m s⁻¹, which could easily destroy the stability of the water column and replenish oxygen into the bottom waters. Indeed, when tropical cyclones impacted the NSCS over the last 4 decades, the local maximum wind speeds typically reached over 9 m s⁻¹ and often were larger than 15 m s⁻¹ (Fig. 9b); the local wind direction was inclined to be from the east, accounting for 56.2 % in frequency, followed by the north (19.8 %) and the west (16.0 %) (Fig. 9d). The long-lasting easterly winds (Fig. 9e) were more likely to confine the river plume to the coast (Xu et al., 2019) and even force the riverine freshwater to subduct down to the depth (Fig. S2), strengthening the reaeration of the oxygen-poor bottom waters (Wang et al., 2017b, 2018; Huang et al., 2019). We also found that the time interval between two successive tropical cyclones was mostly less than 15 d, especially from July to September (Fig. 9c), close to the timescale for hypoxia formation. The frequency of the time interval < 15 d showed a decreasing trend over 1975–2019 (Fig. 9c). Specifically, it averaged 4.1 yr⁻¹ before the year of 2000 in contrast to 2.8 yr⁻¹ after 2000. This might partly explain why less hypoxia was formed and even observed before the year of 2000 (Yin et al., 2004; Rabouille et al., 2008). Therefore, frequent disturbance by tropical cyclones is one of the vital controls on the intermittent hypoxia in low-latitude river-dominated ocean margins.

Figure 9 Statistics of tropical cyclones passing the northern South China Sea (NSCS) from May to September over 1975–2019. (a) Numbers of tropical cyclones. TD, TS, STS, TY and STY represent tropical depressions (the maximum wind speed near the centre is between 10.8–17.1 m s⁻¹ over its lifetime), tropical storms (17.2–24.4 m s⁻¹), strong tropical storms (24.5–32.6 m s⁻¹), typhoons (32.7–41.4 m s⁻¹) and strong typhoons (41.5–50.9 m s⁻¹), respectively. (b) The maximum wind speed of each tropical cyclone. The black line and grey shadow denote the annual average and range of the maximum wind speeds, respectively. (c) The time interval between two successive tropical cyclones. The black line and grey shadow denote the annual average and range of the time intervals, respectively. The negative time interval indicates an overlap in time of two successive tropical cyclones. (d) The wind rose of the intensity of tropical cyclones. (e) The wind rose of the duration of tropical cyclones. The wind speed in (b) and wind direction in (d) and (e) were recorded at Waglan Island (Fig. 1b).

Although strong winds constrain the development of hypoxia, tropical cyclones potentially promote intensive oxygen depletion after the transient dissipation of hypoxia (Rabalais et al., 2009). Heavy precipitation delivered by storms usually increases riverine freshwater loading to the coastal ocean (Zhou et al., 2012), resulting in intensified stratification when winds weaken (Wilson et al., 2008; Su et al., 2017). Indeed, heavy precipitation happened in the Pearl River Delta region after the landing of typhoon Son-Tinh, leading to a wide spreading of the river plume with lower salinity (< 15) than that during leg 1 and the intensification of stratification (Fig. 6), even though the river discharge of the PRE did not increase significantly (Li et al., 2021). Enhanced vertical mixing and/or freshwater discharge supplied large amounts of nutrients to the surface layer to fuel phytoplankton blooms following large storms (Zhao et al., 2009; Ni et al., 2016; Wang et al., 2017a), as shown in Fig. 3, such that strong blooms occurred in the surface plume along the coast with much higher Chl a concentrations during leg 2 than during leg 1. The fresh autochthonous organic matter, together with the resuspended sedimentary organic carbon, provides sufficient substrates for microbial respiration in a re-stratified water column, leading to renewed or even exacerbated bottom water oxygen depletion (Zhou et al., 2012; Song et al., 2020). Along the east coast of North America, lowered DO concentrations were observed after storms (Paerl et al., 2000; Tomasko et al., 2006). Hypoxia was re-established across a larger area when Hurricane Katrina crossed the southeast Louisiana coast (Rabalais et al., 2009). Off the PRE, we also found that hypoxia reoccurred in the wake of a more extensive freshwater plume and enhanced eutrophication after the passage of typhoon Son-Tinh (e.g. DO ∼ 45 µmol kg⁻¹ at stations F202 and F302). Whether it can develop into more severe hypoxia compared to that found initially during leg 1 depends on the net OCR and water column stability, up until the passage of the next storm, typhoon Bebinca (Fig. 1d).

In a changing climate, tropical cyclone features are believed to have a trend with higher intensity but lower frequency (Knutson et al., 2010; Chu et al., 2020). In the NSCS, the annual mean number of tropical cyclones in each decade decreased from 1950–2019, but the number of strong typhoons – the maximum wind speed near the centre is between 41.5–50.9 m s⁻¹ over its lifetime – increased in the last decade (Table 2). The time interval between two successive tropical cyclones therefore increased by 2–3 d in the last 2 decades than prior to 2000 (Fig. 9c). The less-frequent disturbance in the stability of the water column by tropical cyclones and the elongated time interval between two successive tropical cyclones likely favour more persistent hypoxia. However, intensified tropical cyclones would destroy the hypoxia more completely. Despite stronger blooms after tropical cyclones as during leg 2 (Fig. 3k), the winds shifting back to prevailing southwest monsoon winds might drive the blooms offshore (Fig. 3l), reducing the downward transport of fresh, labile organic matter to fuel intensive oxygen consumption in subsurface waters. In this sense, the development of coastal hypoxia depends on a trade-off between the intensity and frequency of tropical cyclones in a warmer ocean.

Table 2 Summary of average frequency of tropical cyclones in each decade from 1950–2019. TD, TS, STS, TY and STY represent tropical depressions (the maximum wind speed near the centre is between 10.8–17.1 m s⁻¹ over its lifetime), tropical storms (17.2–24.4 m s⁻¹), strong tropical storms (24.5–32.6 m s⁻¹), typhoons (32.7–41.4 m s⁻¹) and strong typhoons (41.5–50.9 m s⁻¹), respectively.

5 Conclusions

We have demonstrated the evolution of intermittent hypoxia in summertime as disturbed by the typhoon passage on the inner NSCS shelf off the PRE and examined the controls on maintenance, destruction and reinstatement of hypoxia in this dynamic river-dominated marginal system. Eutrophication-induced hypoxia off the PRE was exacerbated with an enlarged area of ∼ 660 km² and the lowest ever recorded regional DO concentration of 3.5 µmol kg⁻¹ (∼ 0.1 mg L⁻¹). This area of hypoxia was likely underestimated due to spring–neap tidal oscillations, which has recently been demonstrated by a numerical model (Li et al., 2020). Freshwater inputs suppressed wind-driven turbulent mixing and stabilized the water column, restricting the ventilation of the subsurface water and facilitating the formation and maintenance of hypoxia. We estimated for the first time the upwelling-induced oxygen decline and in situ OCR over the destruction and reinstatement of hypoxia, which took place on a timescale of 6–12 d. This in situ OCR could help constrain the source and sink terms for modelling oxygen dynamics in subtropical coastal waters. As the hypoxia timescale is shorter than water residence time but comparable with the disturbance of hypoxia from frequent tropical cyclones or high-wind events throughout the summer season, it could largely explain the intermittent nature of hypoxia off the PRE. Comparing with the hypoxia off large river systems at higher latitudes, tropical cyclones greatly influence the evolution of coastal hypoxia in tropical/subtropical regions. Nevertheless, how the changes in intensity and frequency of tropical cyclones would impact the development of coastal hypoxia in a warmer ocean clearly warrants further studies with carefully designed long-term observations and coupled physical–biogeochemical modelling.