Interannual variabilities, long-term trends, and regulating factors of low-oxygen conditions in the coastal waters off Hong Kong

. The summertime low-oxygen conditions in the Pearl River Estuary (PRE) have experienced a significant expansion in spatial extent associated with notable deoxygenation in recent decades. Nevertheless, there is still a lack of quantitative understanding of the long-term trends and interannual variabilities in oxygen conditions in the PRE as well as the driving factors. Therefore, the long-term deoxygenation in the a subregion of the PRE (the coastal waters off Hong Kong) was 15 comprehensively investigated in this study using monthly observations during 1994-2018. To evaluate the changes in scope and intensity of oxygen conditions, an indicator (defined as the Low-oxygen Index, LOI) that integrates several metrics related to low-oxygen conditions was introduced as the result of a principal component analysis (PCA). Moreover, primary physical and biogeochemical factors controlling the interannual variabilities and long-term trends in oxygen conditions were discerned, and their relative contributions were quantified by the multiple regression analysis. Results showed that the regression models 20 explained over 60% of the interannual variations in LOI. Both the wind speeds and concentrations of dissolved inorganic nitrogen (DIN) played a significant role in determining the interannual variations (by 39% and 49%, respectively) and long-term trends (by 39% and 56%, respectively) in LOI. Due to the increasing nutrient loads and alterations in physical conditions (e.g. the long-term decreasing trend in wind speeds), coastal eutrophication was exaggerated and massive marine-sourced organic matter was subsequently produced, thereby resulting in an expansion of intensified low-oxygen conditions. The 25 deteriorating eutrophication has also driven a shift in the dominant source of organic matter from terrestrial inputs to in situ primary production, which has probably led to an earlier onset of hypoxia in summer. In summary, the Hong Kong waters has have undergone considerable deterioration of low-oxygen conditions driven by substantial changes in anthropogenic eutrophication and external physical factors.


Introduction
Dissolved oxygen (DO) plays a vital role in maintaining the good functioning of aquatic ecosystems. Hypoxia (DO < 2 mg/L) could lead to a marked reduction in habitat for aquatic organisms (Ludsin et al., 2009) and imposes detrimental effects on ecosystem community structure and energy flow (Diaz and Rosenberg, 2008). In recent decades, long-term exacerbation on hypoxia in terms of its spatial extent and intensity has been documented in estuaries and coastal waters worldwide, including 35 the Baltic Sea (Conley et al., 2011;Meier et al., 2019), the northern Gulf of Mexico (Obenour et al., 2013;Laurent and Fennel, 2019), Chesapeake Bay (Li et al., 2016;Ni et al., 2020), the Yangtze River Estuary (Zhu et al., 2011;Zhang et al., 2021), and the Pearl River Estuary Hu et al., 2021). In addition, changes in the phenology of hypoxia were also reported.
For example, in Chesapeake Bay, hypoxic volume has shown a significant increase in early summer but a slight decrease in late summer since 1985(Murphy et al., 2011Testa et al., 2018) . Zhou (2014) also found that the timing of maximum hypoxic 40 volume in Chesapeake Bay was advanced from late July to early July during 1985-2010. A great number of studies have indicated that the exacerbation of hypoxia in coastal systems was closely related to human activities, such as urbanization and industrialization (Breitburg et al., 2018). Due to the anthropogenic influence, massive organic matter and nutrients were discharged into estuaries and coastal waters. Terrestrial organic matter could lead to intense microbial respiration (Rabalais et al., 2010) and excessive nutrient inputs could further stimulate the growth of phytoplankton 45 and exacerbate eutrophication, with a dramatic increase in oxygen demand from marine-sourced organic matter (Fennel and Testa, 2018). Meanwhile, physical processes such as stratification (Rabalais et al., 1991), convergence and migration of water masses (Li et al., 2021), and upwelling (Feng et al., 2014) could regulate the spatial extent and intensity of hypoxia as well.
These processes are closely linked to wind forcing and freshwater discharge (Feng et al., 2012;Yu et al., 2015). In general, the physical and biogeochemical processes exert joint impacts on the generation and development of hypoxia, but different 50 mechanisms may predominate in different systems due to their distinctive natural conditions (e.g. topography) and pressure from anthropogenic pollution. Ni (2020) quantified the contributions of estuary warming, sea level rise, and nutrient load reduction to the long-term changes in hypoxia in Chesapeake Bay through numerical simulation experiments, suggesting that ocean warming was the dominant factor. Forrest (2011) investigated the effects of various processes on the interannual variations of hypoxia in the northern Gulf of Mexico by statistical methods and pointed out that the east-west winds and 55 nutrient loads each accounted for a considerable contribution. While in the Yangtze River Estuary, studies showed that vertical density stratification, which was heavily influenced by a combination of freshwater inputs, various water masses, and winds, was the key factor controlling the interannual changes in hypoxia (Chi et al., 2020).
With the rapid socioeconomic development, the Pearl River Estuary (PRE) has received a large amount of pollutants and nutrients, resulting in a series of environmental problems, including eutrophication, red tide, and hypoxia Li 60 et al., 2020). Since the 1980s, low-oxygen (DO < 4 mg/L) and hypoxic conditions have been reported in the upper reach of Lingdingyang Bay Cui et al., 2018;Hu et al., 2021), Modaomen Bay, Huangmaohai Bay (Su et al., 2017;Shi et al., 2019;Wang et al., 2017;Zhang and Li, 2010) and coastal waters adjacent to Hong Kong (Yin et al., 2004;Su et al., 2017;Shi et al., 2019). Previous studies have shown that hypoxia in the PRE typically occurred in the bottom waters during summer (Yin et al., 2004), driven by strong stratification and sediment oxygen consumption (Zhang and Li, 2010;Wang et al., 65 2017). Due to relatively shallow topography, short water residence time (Rabouille et al., 2008) and short maintenance of stratification (Luo et al., 2009;Lu et al., 2018) hypoxia in the PRE appeared to be episodic and localized (Rabouille et al., 2008). However, this long-standing point of view has been challenged by recent observations showing the emergence of large low-oxygen and hypoxic extents. The area affected by low oxygen in the bottom waters of the PRE was estimated to be around 1,000 km 2 in 2010 (Wen et al., 2020) and ~1,500 km 2 in 2015 . With the increasing availability of observations, 70 an apparent expansion of hypoxia with large interannual variations has been revealed from the data during 1976-2017 (Hu et al., 2021). Nevertheless, due to the scarcity of observations in both time and space and significant differences in sampling periods and locations (sometimes the water quality measurement methods as well) between available datasets, a clear understanding of the long-term trend and interannual changes in hypoxia in the PRE as well as the associated drivers is still lacking, especially from a quantitative perspective. 75 In this study, we utilize observational oxygen and related data collected by the Hong Kong Environmental Protection Department (HKEPD) at certain coastal sites off Hong Kong (see details in section 2.1 below), which have significant merits in terms of temporal coverage (~30 years) and consistency of sampling locations, to perform a quantitative analysis on the long-term oxygen changes (trend and interannual variability) in the region. Moreover, we also aim to discern the key factors controlling the interannual variability and long-term trends in the low-oxygen conditions and to quantify the relative 80 contribution of each primary factor using multiple regression models (Murphy et al., 2011;Forrest et al., 2011;Wang et al., 2021). It is important to note that the HKEPD data with good spatiotemporal continuity allowed us to better estimate the longterm deoxygenation in the coastal waters off Hong Kong, which was close to a hotspot area of low-oxygen conditions in the eastern PRE (Hu et al., 2021) and subject to frequent occurrences of low-oxygen and hypoxic events as well (Yin et al., 2004;4 1994 as well as total organic carbon (TOC) and total nitrogen (TN) measured in the sediments during 1998-2018. The survey stations can be divided into three subregions: (1) the northwestern subregion, including stations NM5 (with water depth of 20 m), NM6 (5 m), and NM8 (8 m); (2) the southern subregion, including stations SM20 (7 m), SM17 (12 m), SM18 (21 m), and SM19 (24m); and (3) the eastern subregion, including stations MM8 (31 m), MM13 (28 m), and MM14 (25 m). Water samples were collected from the surface (1 m below the sea surface), middle (half of the depth at each station), and bottom (1 100 m above the sediments) layers, respectively. Details on the sampling procedures and measurements were described in Xu et al (2010).
In addition, the monthly data of wind speeds and directions used for analysis were estimated using the daily wind observations during 1994-2018 provided by the Waglan Island automatic weather station (Figure 1) of the Hong Kong Observatory. It should be noted that the duration of southwestern winds was defined as the number of its occurrence in days 105 during summer. As for the freshwater inputs from the Pearl River, the monthly data during 1994-2018 were calculated using the discharge data obtained from three major hydrological stations (i.e. Gaoyao, Shijiao, and Boluo) of the Pearl River Water Resources Commission of the Ministry of Water Resources.

Statistical methods
Several metrics, including the cross-sectional area and the layer thickness of low oxygen (DO < 4 mg/L), oxygen 110 deficiency (DO < 3 mg/L) and hypoxia (DO < 2 mg/L) as well as the mean and minimum DO concentrations in the bottom waters, were used to depict the oxygen conditions in the region. Firstly, the observed DO profiles were interpolated by "Natural-Neighbor" method through MATLAB along the three subregions with a grid resolution of 600 m (distance) × 0.3 m (depth). The total areas of DO below 4 mg/L, 3 mg/L, and 2 mg/L were then calculated as the cross-sectional areas of low oxygen, oxygen deficiency, and hypoxia, respectively. The associated layer thickness was defined as the averaged thickness 115 of the grids with DO below the corresponding levels (i.e. 4 mg/L, 3 mg/L, and 2 mg/L). Regarding the island between stations NM8 and SM20, the spatial interpolations were performed directly with all the observed data and then the areas covered by the island were masked out roughly based on its size (Figures 2, 3, A14), as the topographic data of the island was not available.
Such a treatment has little influence on the estimation of vertical low-oxygen areas because low-oxygen conditions were seldom found in stations NM8 and SM20. Moreover, the same treatment procedure was applied to the data in each month of 120 25 years to generate an interpolation set for every month, making it make it consistent when investigating the interannual variations in low-oxygen conditions.
In order to investigate the main variation (interannual changes) of oxygen conditions, we have introduced an indicator integrating the above metrics (except the hypoxic area and thickness, Table A1) through the PCA (principal component analysis) technique, which can reduce the dimensionality of a dataset to make it more interpretable with minimum information 125 loss (Cadima et al., 2016). The two metrics related to hypoxia were excluded from PCA because the occurrence of hypoxia was relatively rare and its interannual variation was not as significant as that of low oxygen and oxygen deficiency. The results of PCA analysis (Table A2) showed that the first component explained most of variance (86.40%) for the six input variables, while the remaining components explained less variance (13.60%). The first component was highly correlated with the interannual variations of the cross-sectional areas (with a correlation coefficient r of 0.96, p < 0.01) and the thickness (r = 0.96, 130 p < 0.01) of low oxygen as well as the bottom DO concentrations (r = -0.90, p < 0.01), and it was thereafter referred as Lowoxygen Index (LOI, Equ. 1) to describe the interannual severity of low-oxygen conditions comprehensively. where DOmean and DOmin represent the mean and the minimum DO concentrations in the bottom waters, respectively; Area4 (Area3) and Thickness4 (Thickness3) represent the cross-sectional area and the thickness of low oxygen (oxygen deficiency), 135 respectively.
As the low-oxygen conditions within Hong Kong waters were jointly affected by physical and biogeochemical processes, we attempted to quantify the relative contributions of multiple relevant factors including wind, freshwater, water temperature, and nutrients to interannual variability and long-term trends of the oxygen conditions through multiple regression. As for the selection of the wind variable in use, the daily wind data were processed into monthly average wind speed (WS), southwestern 140 wind duration (SWWD), southwestern wind cumulative stress (SWCS), and southeastern wind cumulative stress (SECS) in summer (June-August) to examine the effect of wind speed and direction ( Figure A2). Then, a suite of multiple regressions was carried out to fit the LOI for each wind-related variable. As shown in Table A3, the fitting effect of LOI was better when using WS, which also has the highest correlation with LOI among the wind-related variables, revealing that WS explained the most interannual variation of LOI among the wind-related factors. Therefore, WS was eventually adopted to be the wind-145 related input variable in the multiple regression with freshwater discharge (flow), the monthly spatial-average of bottom temperature (T), and surface DIN in the summer. The resulting regression coefficients were then standardized by multiplying the ratio between the standard deviation of each input variable (e.g. WS) and the standard deviation of LOI to evaluate their interannual contributions (Equ. 2).

= ×
(2) 150 Where Csti and Ci represent the regression coefficients of WS, flow, T, and DIN after and before standardization, respectively; SDi represent the standard deviation of WS, flow, T, and DIN; SDLOI represents the standard deviation of LOI.
In addition, the dataset was randomly split into a training dataset (70%) and a testing dataset (30%) in order to provide a more robust data fitting with estimates on the uncertainties arising from different data selections. Consequently, over 480,700 combinations of training and testing datasets were generated randomly from this splitting process and were used to build up a 155 variety of regression models. Coefficient of determination (R 2 ) was used to measure the fitting effect in training and testing datasets. Of all the established models, the fitting effect of training datasets (e.g. R 2 train) and coefficients of the four variables were similar, but the predictive skills in testing dataset (e.g. R 2 test) varied in a large range ( Figure A35, Table A4). Besides, larger standard deviation occurred in coefficients in cases with worse testing effects. To provide a more robust estimation for the fitting, only those with R 2 over or equal to 0.6 both for the training and testing datasets were selected to quantify the impact 160 of each input variable according to their regression coefficients on average ( Figure A42). Furthermore, based on the selected models (with R 2 ≥ 0.6 for both datasets), we also set up four sensitive experiments in which the long-term trend of each input variable was removed and only interannual fluctuations were retained. The LOI was then re-calculated in each scenario and its change relative to the original LOI was used to assess the impact of each variable to the long-term oxygen trend. In terms of spatial distributions, distinct differences were observed for the hydrological and eutrophication parameters among the three subregions investigated. Due to the profound influence of river discharges, temperature/salinity in the northwestern subregion (NM5-NM8, closer to the river outlets) was noticeably higher/lower when compared to the other two (Figure 2 a-d), varying by 28.61±1.09 °C/14.63±6.23 PSU at the surface in summer. Meanwhile, the DIN concentration in the 190 northwestern subregion was the highest (Figure 3 a-b), reaching up to 1.17±0.40 mg/L at the surface. On the contrary, the eastern subregion (MM8-MM14), which was farthest to the river outlets and more heavily affected by the shelf water, had the lowest temperature (27.85±1.27 °C), highest salinity (29.22±3.10 PSU) and lowest DIN concentration (0.14±0.13 mg/L) in the surface waters. As for Chl a (Figure 3c-d), the highest level appeared in the southern subregion (SM17-SM20, with 10.19±8.86 μg/L at the surface in summer), while the lowest one was found at the northwestern subregion (with 6. 82±10.67 195 μg/L at the surface).  Figure A53).

Dissolved oxygen and low-oxygen conditions
In addition to the DO levels, we also investigated the interannual changes in the summertime low-oxygen conditions in terms of areal extents (vertical profiles), thickness, and the LOI as defined in section 2.2 ( Figure 4). Significant interannual 210 fluctuations were found for all these metrics; for example, the area and thickness affected by low oxygen fluctuated between  (Figure 4). More specifically, before 2000 the spatially-averaged DO concentrations in the bottom waters exceeded 4 mg/L and low-oxygen conditions were seldom observed (Figure 4a), while the DO minimums were all above 2 mg/L (i.e. no hypoxic events occurred). However, since 2000 the occurrences of low oxygen and hypoxia have become more 225 frequent, with a significant growth in the LOI and its related metrics, confirming the exacerbation of low-oxygen conditions in the Hong Kong waters.
To further quantify the intensity of long-term deoxygenation in summer, linear regressions were performed for the DO concentrations in different layers and in different subregions and also for the areal extents of low-oxygen conditions during 1994-2018 ( Figure 5). As shown, apparent declining trends were found for the DO series both at the surface (although not 230 significant, Figure 5a) and the bottom (Figure 5b). For the bottom waters, the averaged DO concentrations displayed a decreasing pattern with a rate of 0.03 mg/L per year (equivalent to approximately 0.7% of the climatological DO mean at the bottom), while the DO minimums showed a more significant decline with a rate of 0.08 mg/L per year (~3.5% of the climatological mean of the bottom DO minimums). It was also noted that the intensity of deoxygenation varied between subregions (Figured 5c-h). As for the bottom DO concentrations, the most significant decrease was found in the eastern 235 subregion (with a deoxygenation rate of 0.05 mg/L per year, Figure 5h), while the most significant decline in the DO minimum appeared in the southern subregion (with a rate of 0.08 mg/L per year, Figure 5f). Likewise, significant increasing trends were also found for the areas of low oxygen and oxygen deficiency (Figure 5i-j), showing an annual growth rate at 1.95×10 4 m 2 and 4.75×10 3 m 2 , respectively. Regarding the changes in LOI, it had a growth rate of 0.20 per year, which corresponds to an increasing rate of 1.99×10 4 m 2 in the low-oxygen area and a declining rate of 0.07 mg/L in the DO minimum. 240 Furthermore, the long-term oxygen changes varied between different months of the summer season as well ( Figure 6). It could be seen that the decreasing magnitude of the averaged DO concentration was close to each other for all the summer months, while the decline in the DO minimum was most pronounced in July (with a decreasing rate of 0.10 mg/L per year, Figure 6c), followed by that in August (0.06 mg/L per year, Figure 6e). In fact, the long-term changes in the DO minimum had different patterns in July and August. As for July, the DO minimum generally showed a consecutive decrease over the past 25 245 years (Figure 6d). While in August, the DO minimum experienced a rapid decline with a rate of 0.14 mg/L per year during 1994-2011, which was higher than that in July during the same period (0.11 mg/L per year), but subsequently undertook a recovery from the hypoxic conditions since 2012 (Figure 6f). Along with such distinctive intra-seasonal patterns, an interesting phenomenon was also noticed: hypoxic events were present mostly in August prior to 2012 (e.g. in 2007 and 2010-2011; no hypoxia was found in July during the same period) but only in July instead since 2012 (e.g. in 2014 and 2016-2017), as shown 250 in Figure 6. This finding implied a potential shift in the onset of hypoxia generation from August to July, i.e. an earlier timing for the arrival of the summertime hypoxia. Accordingly, distinct changes were found for the areas affected by hypoxia in the two periods around 2012. The hypoxic area estimated in July increased from zero during 1994-2011 to (5.42±8.77)×10 3 m 2 during 2012-2018, whereas the hypoxic area in August decreased from (0.89±2.82)×10 4 m 2 to zero.

Discussion 255
4.1 Primary factors controlling the interannual variabilities in low-oxygen conditions As shown above, significant interannual variabilities were observed in the spatial extent (e.g. cross-sectional area) and intensity of oxygen conditions (e.g. the mean bottom DO concentrations). Such variabilities were largely influenced by multiple physical and biogeochemical factors, including wind forcing, freshwater discharge, water temperature and nutrient loads. These processes jointly act to affect density stratification (Yu et al., 2015), water residence time (Li et al., 2021) and 260 temporal and spatial distributions of eutrophication parameters (Cui et al., 2018). As described in section 2.2, four important influential factors (i.e. WS, flow, T, and DIN) were used to predict the interannual variations in LOI by the multiple regression models, in which there have been 56,010 cases (~12% of the total, Figure 7) with R 2 ≥ 0.6 both in the training dataset (mean R 2 of 0.64) and the testing dataset (mean R 2 of 0.70). The standardized coefficients (mean±standard deviation) for these wellperforming regression cases were given as follows: 265 LOI = −(0.39 ± 0.12) × WS − (0.14 ± 0.12) × flow − (0.11 ± 0.08) × T + (0.49 ± 0.12) × DIN As denoted by the regression coefficients, wind forcing has exerted a significant impact on the interannual changes in LOI, with a relative contribution of 39%±12% to the LOI variability explained. Its importance could also be evidenced by the significant negative correlation between WS and LOI (r = -0.67, p < 0.01, Figure 8), suggesting that calm winds were beneficial to low-oxygen conditions. In most cases (e.g. in Modaomen Bay in the PRE and the northern Gulf of Mexico), strong winds 270 could break down stratification in the water column (Rabalais et al., 1991;Feng et al., 2012), which was conducive to water mixing and atmospheric reoxygenation (Rabalais et al., 1991). However, the weak correlation between WS and Δρ ( Figure 8) indicated that the wind forcing may control hypoxia through other alternative mechanisms. Actually, weak winds in combination with flow convergence induced by wind-driven circulation could contribute to long water residence time and nutrient accumulation in the eastern PRE and thus favored the phytoplankton blooms (Li et al., 2021). This could be supported 275 by the significant negative correlation between WS and Chl a (r = -0.62, p < 0.01). In contrast, the wind direction showed less significant effect on the interannual variability in low-oxygen conditions, as suggested by the comparatively poor performance in the LOI fitting and the weaker correlations of the wind direction-related variables with LOI (Table A3). It was noted that the monthly average wind direction in summer were generally southerly with small changes (mostly varying between 150° and 200°, Figure A1A2). Overall, our results indicated that the wind speed played a more important role in regulating the low-280 oxygen conditions in the coastal waters off Hong Kong from an interannual perspective, although the wind direction could significantly influence the short-term generation and development of low-oxygen conditions by modulating the Pearl River plume and material fluxes (Yin et al., 2004;Li et al., 2021). With respect to the DIN concentrations, it played a vital role in determining the interannual variabilities of the oxygen conditions, with a contribution up to 49%±12%. It has been widely recognized that eutrophication stimulated by anthropogenic nutrient inputs could provide a large quantity of depositing detritus 285 and subsequently led to substantial oxygen depletion and occurrence of low-oxygen events (Rabalais et al., 2010;Fennel and Testa, 2018); for example, in the northern Gulf of Mexico (Feng et al., 2012;Forrest et al., 2011) and Chesapeake Bay (Wang et al., 2015), the interannual hypoxic areas in summer were directly regulated by the nutrient levels. Similar situation was found in the PRE  and Hong Kong waters, as confirmed by the significant positive correlation between DIN and LOI (r = 0.65, p < 0.01). Collectively, DIN and WS were identified as the two key factors controlling the interannual 290 changes in low-oxygen conditions. Compared to WS and DIN, the freshwater discharges (flow) had a much smaller contribution (~14%±12%) to the variations in LOI. Generally speaking, large freshwater inputs tend to enhance the intensity of water stratification and facilitate the generation of hypoxia (Rabalais et al., 1991). However, we found a negative correlation between flow and LOI (r = -0.45, p < 0.05, Figure 8), implying that the effect of freshwater discharges on low-oxygen conditions might involve more complex 295 mechanisms and act through indirect pathways. Due to its long distance from the river outlets of the Pearl River, the coastal waters off Hong Kong were relatively less influenced by terrestrial inputs  and the effect of freshwater discharge and its carrying organic matter in this area was not as significant as that in other subregions (e.g. the upper reach of Lingdingyang Bay and the western PRE). Nevertheless, freshwater discharge in combination with the wind-driven circulation could significantly affect the water residence time (Sun et al., 2014) and nutrients accumulation in the Hong Kong waters (Li 300 et al., 2021). Specifically, the weakened discharge could prolong the retention of nutrients and thereby stimulate local productions of organic matter in the region (Li et al., 2021), which ultimately promoted oxygen depletion. Regarding the water temperature, previous studies have showed that it could exert significant influence on coastal hypoxia largely by regulating water stratification intensity, oxygen solubility, and microbial respiration rate (Breitburg et al., 2018). However, our results showed that the contribution of water temperature (T) to the LOI changes (~11%±8%) was not significant in the Hong Kong 305 waters, as revealed by its weak correlation with LOI as well. Given the fact that the Hong Kong waters are a region heavily affected by human activities, the effect of temperature (e.g. global warming) might be more significant in the region with larger geography scale. Overall, the role of temperature and freshwater discharges in regulating the interannual oxygen variability in the Hong Kong waters appeared to be secondary.

Drivers of the long-term deoxygenation trend 310
The data over the past 25 years showed that the coastal waters off Hong Kong has experienced a notable long-term oxygen decline, especially for the DO minimum in the bottom waters. Based on the observed deoxygenation rate, the bottom DO minimum was expected to decrease by approximately 15%-70% in 5-20 years (reaching a level of 0.4-1.6 mg/L) compared to the climatological mean of 1994-2018. The impacts of influential factors on the long-term deoxygenation trend were then evaluated using the regression models mentioned in section 4.1 and quantified by the relative changes of LOI in the sensitive 315 experiments (see details in section 2.2) compared to the original one. It was noted that WS exhibited a decreasing trend of 0.03 m/s per year (p < 0.05, Figure 9a) within the coastal regions off Hong Kong over the past 25 years, while similar situation was also found in the Pearl River Basin  and the northern South China Sea (Gao et al., 2020) due to the longterm climate changes (Xu et al., 2006;Zhang et al., 2009;Chen et al., 2020). Meanwhile, DIN showed an increasing trende with a rate of 0.01 mg/L per year (p < 0.01, Figure 9e). The growth in DIN and decline in WS have led to a 56%±10% and 320 39%±14% increase in LOI (Table 1), respectively, indicating that DIN and WS were the main driving factors for the long-term deoxygenation. On the other hand, significant long-term trends were also found for the freshwater discharges (with a decreasing rate of 4.19×10 2 m 3 /s per year, Figure 9b) and water temperature (with an increasing rate of 0.06 °C per year, Figure 9d), but their impacts were relatively small, resulting in a 16%±14% increase and 11%±9% decrease in LOI, respectively.
Despite the different influences of the factors mentioned above, they were likely to impose synergetic impacts on the low-325 oxygen conditions by aggravating eutrophication as discussed earlier; it could be observed that the long-term growth in Chl a (with a rate of 0.15 μg/L per year, Figure 9f) matched well with the increase in LOI. Specifically, the significant increase in phytoplankton biomass was primarily due to the combined effects of more stable water-column condition and longer residence time facilitated by the weaken wind forcing and river discharges, higher nutrient levels, and lower water turbidity (Figure 9g) in recent years. Consequently, the elevated organic matter through phytoplankton primary production would lead to strong 330 oxygen consumption, thereby contributing to an expansion of low-oxygen conditions in terms of areal extent and intensity.
Moreover, with massive algal fragments provided by primary production, the composition of organic matter in the coastal waters off Hong Kong has probably changed and would cause substantial changes in the timing of hypoxia generation. As noted in section 3.2, the onset of hypoxia was observed to shift from August to July around 2012. To explore this issue, we first used the ratio of TOC to TN measured in the sediments to estimate the main source of organic matter, with values of 14-335 30 pointing to a terrestrial source and values of 4-10 indicating a marine source from in-situ production (Bordovskiy, 1965;Meyers and Ishiwatari, 1993). It is clear that the TOC:TN showed a significant decreasing trend and was mostly below 10 since 2012 (Figure 9h). This implied a shift in the dominant source of organic matter from terrestrial inputs to local production (marine-sourced). As such, oxygen consumption became faster because the marine-sourced organic matter was fresher and more active (Raymond and Bauer, 2001) and therefore the time required to reach hypoxia would be shortened. Furthermore, 340 changes in the physical conditions provided sufficient time for more thorough decomposition of organic matter in July, which left less organic matter for August and thus weakened the deoxygenation therein.
Similarly, the long-term oxygen changes in terms of the areal extents and arrival timing of hypoxia have also been found in other coastal systems. For example, in Chesapeake Bay, sea level rise and elevated freshwater discharges would lead to an approximately 10%-30% increase in hypoxic volume between the late 20 th and the mid-21 st centuries (Ni et al., 2019), while 345 the increase in water temperature would cause hypoxia to develop 5-10 days earlier in ~30 years (Ni et al., 2020). In the northern Gulf of Mexico, the growth in riverine nutrient inputs would result in an increase in the frequency of hypoxia occurrence by 37% (Justić et al., 2003). While in the Hong Kong waters, low-oxygen conditions would develop into hypoxic conditions in two decades with larger areal extent and earlier arrival ascribed to the ongoing alterations in physical conditions and nutrients as mentioned earlier. This inference was based on the assumption that the external factors (e.g. wind speed, DIN, 350 discharges) would change at the same rates as those in the past 25 years. Although the real situation would be more complicated and compounded by factors such as the implementation of management and non-linear changes in climatic factors, our findings still served as an alarming signal that changes in wind and freshwater discharges could cancel out potential benefits of nutrient management. To this end, it is of great importance to conduct long-term and more intensive control on nutrient inputs in order to mitigate the low-oxygen conditions in the region. 355

Conclusion
We have comprehensively investigated the spatiotemporal characteristics of DO and various related water quality variables in the coastal waters off Hong Kong and found that low-oxygen conditions occurred mostly in the bottom waters of summer, with significant interannual variability and an apparent deoxygenation trend over the past 25 years. We have also quantified the contribution of each primary factor by statistic methods and found that the increasing DIN levels and the 360 decreasing wind speeds, both of which would eventually lead to the intensification of eutrophication, contributed most to the interannual variations and long-term trend in LOI. Therefore, more marine-sourced organic matter was produced by the elevated primary production, leading to an exacerbation in low-oxygen conditions with larger areal extents as well as a potential earlier onset of the summertime hypoxia. By comparison, the freshwater inputs and water temperature had relatively small impacts on the long-term changes in LOI. To sum up, this study has shown that oxygen conditions in the coastal waters off 365 Hong Kong have been deteriorating under the interactions of altered physical forcing (e.g. winds) and aggravated eutrophication and it would develop into a severe hypoxic state within the next two decades. Lastly, given the significant intraseasonal variability in low-oxygen conditions during summer, it is of great importance to conduct more cruise surveys to collect estuary-wide observations on a longer time scale in order to fully capture the generation and development of hypoxia and to confirm the change in the timing of its arrival. 370 Data availability. The marine water quality and the sediment data during 1994-2018 from the HKEPD are available at https://www.epd.gov.hk/epd/epic/english/epichome.html, while the daily wind observation data from Waglan Island automatic weather station are available at https://www.hko.gov.hk/sc/cis/climat.htm. Daily discharge data of hydrological stations (i.e. Gaoyao, Shijiao, and Boluo) can be collected at http://www.zwswj.com/cms/webfile/waterInfo/index.html. 375 Author contributions. Under the conceptualization of JH, ZC completed the data analysis and graphic visualization. This work was supervised by JH and SL. ZC wrote the paper with contributions from all co-authors and all co-authors contributed to the review and editing of manuscript, especially for JH and BW.

380
Competing interests. The authors declare that they have no conflict of interest. Table 1 Long-term trends in the fitted LOI on average for the selected regression cases with R 2 > 0.6 (baseline) and for the sensitive experiments with respect to the effects of wind speeds (b), freshwater discharges (c), water temperature (d), and surface DIN concentrations (e).       Cross-sectional area of low-oxygen (DO<4 mg/L) of each year during 1994-2018 Area3 Cross-sectional area of oxygen-deficiency (DO<3 mg/L) of each year during 1994-2018 Thickness4 Cross-sectional thickness of low-oxygen of each year during 1994-2018 Thickness3 Cross-sectional thickness of oxygen-deficiency of each year during 1994-2018 Low-oxygen Index (LOI) First principal component of PCA dimension (86.40% of variation) for measuring interannual variations in scope and intensity of oxygen conditions  Table A3. Coefficients of determination (R 2 ) for average wind speed (WS), southwestern wind duration (SWWD), southwestern wind cumulative stress (SWCS), southeastern wind cumulative stress (SECS) in fitting LOI; Pearson correlation coefficient of WS, SWWD, SWCS, and SECS with LOI. Note that the symbols * and ** represents the significant level at p < 0.05 and p < 0.01, respectively.

610
Note that the negative values of SWCS and SECS represent southwestern and southeastern wind, respectively. The trends and significant p values were shown in the title of each subgraph. Figure A53. Combined fitting results of the regression models with different combinations of training and testing datasets. R2 in (a)

615
were greater than or equal to 0.6 both in training and testing datasets. R2 in (b) were less than 0.6 both in training and testing datasets. Fitting results of total samples were in (c). Note that the red hollow dots denote the LOI estimated based on observational data, while the blue solid dots and the gray patch represent the mean values and ranges of the fitted LOI in the selected regression cases, respectively. Figure A42. Flowchart describing the fitting of LOI and the cases sampling used for analysis.