Seasonal Variation of Dissolved Oxygen in the Southeast of the Pearl River Estuary

Dissolved oxygen (DO) concentration in estuaries is highly variable at different spatial and temporal scales, which is affected by physical, chemical and biological processes. This study analyzed the spatial–temporal distributions of dissolved oxygen concentration and bottom hypoxia in the southeast of the Pearl River Estuary (PRE) using monthly water quality monitoring and hydrographic data covering the period 2000–2017. The seasonal spatial–temporal variation of DO concentration was studied using various methods, such as rotated empirical orthogonal functions, harmonic analysis, and correlation analysis. The results showed that DO stratification was significant in summer, but it was not distinct in winter, during which DO concentration peaked. DO stratification exhibited a significantly positive correlation with water stratification. In the south and west of Hong Kong (SHK and WHK, respectively), DO concentration fields exhibited distinct seasonal changes in the recent 18 years. In SHK, the main periods of the surface DO variation were 24, 12, and 6 months, whereas the main period was 12 months in WHK. The main period of the bottom DO variation was 12 months in both SHK and WHK. In SHK, the spatial–temporal variations in surface and bottom DO were highly related to the variations of salinity, dissolved inorganic nitrogen (DIN), and active phosphorus, and the variation of surface DO was also connected to the variation of temperature and chlorophyll a. In WHK, the variations in surface and bottom DO were highly related to the variations of salinity and temperature, and the variation of surface DO was also connected to the variation of DIN. The river discharge and wind had a different important influence on the temporal variability of DO in WHK and SHK. These findings suggested that the variations of DO may be controlled by coupled physical and biochemical processes in the southeast of PRE. From 2000 to 2017, bottom hypoxia in the southeast of PRE occurred in the summers of 7 years. SHK appeared to be more vulnerable to hypoxia than WHK.


Introduction
Oxygen (O 2 ) is fundamental to life and a critical constraint on marine ecosystems [1,2]. Dissolved oxygen (DO) concentration affects marine biogeochemical processes as well as the survival and distribution of marine organisms [3,4]. Under the combined actions of physical, chemical, and biological factors, the distribution of DO concentration shows obvious spatial-temporal differences in coastal and estuarine areas [5,6]. At present, influenced by human activities and climate change, the contents of DO are decreasing sharply, with increasing hypoxia in coastal and estuarine areas [2,4,7]. Hypoxia generally refers to a DO concentration below 2 mg/L [8][9][10], which is extremely hazardous to marine ecosystems [2,7,11]. The variations of DO and hypoxia have been documented in many coastal and estuarine systems worldwide such as the Gulf of Mexico, Chesapeake Bay and the Changjiang (Yangtze River) Estuary [2,12,13].
The Pearl River Estuary (PRE) is an important estuarine ecosystem, which is located in the northern shelf of the South China Sea, with Hong Kong at its eastern coast [10,14,15]. The annual average river discharge is approximately 3.26 × 10 11 m 3 , with 20% occurring during the dry season in the period October-March and 80% during the wet season in the period April-September [10,16]. PRE is a complex estuary owing to its complicated topography and multiple forcing factors, such as river outflows, tides, monsoons, and coastal currents [17][18][19]. In winter, the northeast monsoon prevails and the coastal current dominates the southeast of PRE, while in summer, the interaction of the estuarine plume and oceanic waters dominates due to the southwest monsoon [9,15].
During the last two decades, PRE has been experiencing extremely rapid environmental changes caused by anthropogenic inputs, such as discharges of industrial wastewater and domestic sewage [16,20,21]. Consequently, DO concentration began to decrease in 2000 [9,22]. During the last 25 years, bottom oxygen in PRE decreased by~2 µmol kg −1 year −1 [23]. Since episodic hypoxia was first reported in the 1970s to 1980s, it has increasingly become an important issue of concern for the scientific community and general public [16,24,25]. Hypoxia mainly occurred after 2000 and has aggravated in recent years in the lower reach of PRE [22,23]. Yin et al. [9] investigated the effect of estuarine circulation, sediment oxygen demand and P (phosphorus) limitation on the distribution of DO and the formation of hypoxia. Modeling studies have identified that the occurrence and development of hypoxia are mainly controlled by the salinity front [26], stratification [26,27], sediment oxygen demand [28][29][30], and re-aeration [29]. The evolution of hypoxia is also strongly associated with wind and river discharge [31][32][33]. He et al. [21] explored hypoxia as the combined effect of aerobic respiration and nitrification, which contribute almost equally to oxygen consumption in the upper reaches of PRE. In the lower reach of PRE, terrestrial organic matter significantly contributed to the formation and maintenance of hypoxia [34,35]. During the summer of 2015, a low-DO (<4 mg/L) zone covering an area of~1500 km 2 was observed [36]. Cui et al. [37] identified a mid-depth barrier layer inhibiting vertical exchange between the river plume and shelf salt wedge in the hypoxia zone.
These studies have analyzed the characteristics of DO concentration and provided evidence that both biological and physical processes play important roles in controlling hypoxia in PRE. However, the variability of DO was regulated by various factors in different regions of PRE [22]. Although seasonal variations in DO in some marine waters of Hong Kong were clearly present [24,38], the spatial-temporal characteristics of the seasonal cycles of DO remain unclear in the southeast of PRE. In addition, little information is available on the seasonal characteristics of hypoxia in the southeast of PRE. Accordingly, this study aimed to identify the spatial-temporal characteristics of the seasonal cycles of DO and seasonal characteristics of hypoxia in the southeast of PRE.

Data
The Environmental Protection Department of the Hong Kong Government has been maintaining a comprehensive sampling program for over 30 years to monitor water quality in the territorial waters southeast of PRE. There are 76 water quality monitoring stations covering the marine territorial waters of Hong Kong [39]. In this study, 10 stations representing the southeast of PRE were selected (DM5-MM8, Figure 1). They were located in PRE (DM5-NM8) and in southern Hong Kong (SM20-MM13). Therefore, these 10 stations are representative of the estuarine influence (DM5-NM8) and the estuarine coastal plume (SM20-MM13). Monthly data of DO, temperature, salinity, chlorophyll a (Chl a), dissolved inorganic nitrogen (DIN, sum of NH 4 , NO 2 and NO 3 ), and active phosphorus (PO 4 ) covering the period 2000-2017 were used in this study. A SEACAT19 CTD (Sea-Bird Scientific, Bellevue, WA, USA) was used to obtain vertical profiles of salinity, temperature and other parameters [39]. DO was measured using the SBE23Y DO (membrane electrode) detector (Sea-Bird Scientific, Bellevue, WA, USA) assembled on the CTD [39]. NH 4 , NO 2 , and NO 3 were measured by the flow injection analysis method of the American Society for Testing and Materials (ASTM) D3590-89 Bstandard, the American Public Health Association (APHA) 20ed. 4500-NO2 -B standard, and APHA 20ed. 4500-NO3 -F & I standard, respectively [39][40][41]. PO 4 was measured by the flow injection analysis method of ASTM D515-88 A standard, and Chl a was determined spectrophotometry using APHA 20ed. 10200H 2 standard [39][40][41]. This dataset has been used in several studies over the past ten years [11,[42][43][44].  Figure 1. Sampling stations (dots in black) in the southeast of the Pearl River Estuary. PRE and N represent the Pearl River Estuary and Nei Lak Shan, respectively. Blue and red triangles represent the river discharge gauges station (Gaoyao, Shijiao, and Boluo) and automatic weather station, respectively.
The mean monthly river discharge data were obtained from gauge stations of the Pearl River Water Resources Commission of the Ministry of Water Resources (PRWRC) at Gaoyao, Shijiao, and Boluo ( Figure 1). The sum of the river discharges at Gaoyao, Shijiao, and Boluo was used to represent the total river discharge entering PRE. The mean monthly wind speed and direction data were obtained from the Nei Lak Shan automatic weather station of the Hong Kong Observatory (Figure 1). The alongshore wind speed was calculated by rotation into a coordinate system aligned with the local coastline with a rotation angle of 25 • (between due east and the alongshore direction measured in a counter-clockwise direction).

Methods
In order to analyze the spatial-temporal characteristics of DO concentration, the data for each sampling month during the period 2000-2017 were averaged to obtain monthly and annual means of DO, salinity, temperature, Chl a, DIN, and PO 4 . Subsequently, the spatial-temporal variations of their characteristics on seasonal cycle were analyzed using rotated empirical orthogonal functions (REOF). From among various types of REOF schemes [45], the Varimax REOF analysis [46], which linearly transforms spatial patterns derived from EOF analysis into rotated spatial patterns (REOFs), was employed because it maximizes the variance of squared correlation coefficients between each rotated principal component (RPC, rotated principal component) and original EOF mode. After rotation, high loading factors correlated to principal components concentrate only in a small area, and the loading factors in most areas are as close to zero as possible. Therefore, significant anomalies appear where regional phenomena are dominant, which makes the spatial structures easier to detect and interpret [47]. Finally, the main periods and amplitudes of the RPCs were calculated by the maximum entropy spectral analysis and the least squares method, respectively. The correlation of RPCs was analyzed using the method of delay correlation analysis. Figure 2 shows changes in the annual and monthly mean of water temperature and salinity in the southeast of PRE. The annual mean temperatures at the surface and bottom were 24.0 and 22.9 • C, respectively. The temperature difference between the surface and bottom increased with increasing depth (Figure 2a). Surface water temperature was higher in summer, reaching a maximum value of 28.6 • C in July, and lower in winter, with a minimum value of 17.2 • C in February, which could be attributed to solar radiation. The water temperature at the bottom was 27.1 • C in October (autumn) and 17.0 • C in February (winter). The largest difference between the surface and bottom water temperatures was observed in July (summer) at 3.5 • C and the smallest was observed in winter, indicating obvious water stratification in summer ( Figure 2b). From the outside of Deep Bay to the south of Hong Kong (SHK), surface salinity increased from 23.7 to 31.7, and bottom salinity increased from 27.4 to 33.7 ( Figure 2c). Salinity was lowest in June (summer), with values of 20.2 and 29.2 at the surface and bottom, respectively, and highest in December (winter), with values of 32.0 and 32.5 at the surface and bottom, respectively. The difference between surface and bottom salinity was largest in June (summer) at 9.0 and smallest in winter at 0.4, indicating obvious salinity stratification in summer ( Figure 2d). The station MM5 is located near the channel and has a water depth of 20 m. In summer, bottom temperature and salinity at station MM5 were lower and higher than those at the adjacent station, respectively. The salinity decreased with an obvious stratification in summer, which can be attributed to the rainy season. The bottom temperature decreased due to upwelling induced by the southwest monsoon in summer.

Characteristics of Main Environmental Factors
Water 2020, 12, 2475 4 of 18 linearly transforms spatial patterns derived from EOF analysis into rotated spatial patterns (REOFs), was employed because it maximizes the variance of squared correlation coefficients between each rotated principal component (RPC, rotated principal component) and original EOF mode. After rotation, high loading factors correlated to principal components concentrate only in a small area, and the loading factors in most areas are as close to zero as possible. Therefore, significant anomalies appear where regional phenomena are dominant, which makes the spatial structures easier to detect and interpret [47]. Finally, the main periods and amplitudes of the RPCs were calculated by the maximum entropy spectral analysis and the least squares method, respectively. The correlation of RPCs was analyzed using the method of delay correlation analysis.  (Figure 2d). The station MM5 is located near the channel and has a water depth of 20 m. In summer, bottom temperature and salinity at station MM5 were lower and higher than those at the adjacent station, respectively. The salinity decreased with an obvious stratification in summer, which can be attributed to the rainy season. The bottom temperature decreased due to upwelling induced by the southwest monsoon in summer.  The annual mean values of Chl a at the surface and bottom were 4.5 and 2.3 µg/L, respectively. Chl a concentration was relatively low at the bottom of SHK (Figure 3a). In summer, surface Chl a concentration was higher, with the highest value of 9.8 µg/L in July, while in other seasons, both surface and bottom Chl a concentrations were lower, with the lowest value at the surface in winter (Figure 3b). Chl a concentration was lower in winter because low water temperatures were not conducive to phytoplankton growth; the rising temperature in spring promoted the growth of phytoplankton, thus increasing Chl a concentration. In summer, the higher solar radiation and rich nutrients promoted phytoplankton growth, resulting in a peak of Chl a, while in autumn, with reduced temperature and radiation, the phytoplankton died off rapidly, resulting in decreased Chl a concentration.

Characteristics of Main Environmental Factors
Water 2020, 12, 2475 5 of 18 The annual mean values of Chl a at the surface and bottom were 4.5 and 2.3 μg/L, respectively. Chl a concentration was relatively low at the bottom of SHK (Figure 3a). In summer, surface Chl a concentration was higher, with the highest value of 9.8 μg/L in July, while in other seasons, both surface and bottom Chl a concentrations were lower, with the lowest value at the surface in winter (Figure 3b). Chl a concentration was lower in winter because low water temperatures were not conducive to phytoplankton growth; the rising temperature in spring promoted the growth of phytoplankton, thus increasing Chl a concentration. In summer, the higher solar radiation and rich nutrients promoted phytoplankton growth, resulting in a peak of Chl a, while in autumn, with reduced temperature and radiation, the phytoplankton died off rapidly, resulting in decreased Chl a concentration. Generally, nitrogen and phosphorus are absorbed by phytoplankton according to the specific Redfield ratio (16:1) [48][49][50]. The Redfield ratio represents the average value with an optimum proportion for phytoplankton absorbing and using nutrient substances and is widely applied when investigating nutrient limitations of marine phytoplankton. An nitrogen to phosphorus ratio (N:P ratio) higher than 22:1 in seawater suggests a potential P limitation, whereas a ratio below 10:1 suggests a potential N limitation [51]. Figure 4e,f show the N:P ratio of DIN:PO4. The surface N:P ratio in WHK was higher than 22:1, indicating an imbalanced nutrient proportion with excessive N and a potential P limitation. In SHK (SM18-MM13), the bottom N:P ratio was below 10:1, suggesting excessive P and a potential N limitation. From SM20 to the east, this ratio decreased gradually, indicating a shift from P limitation to N limitation. From April to September, the surface N:P ratio was higher than 22:1, indicating excessive nitrogen and a potential P limitation in summer. Except in May and June, the bottom N:P ratio was only slightly higher than 22:1. The surface N:P ratio was higher than 22:1 and the excessive N resulted in an antecedent consumption of P, with phytoplankton growth under P limitation. The considerable water stratification hindered the Generally, nitrogen and phosphorus are absorbed by phytoplankton according to the specific Redfield ratio (16:1) [48][49][50]. The Redfield ratio represents the average value with an optimum proportion for phytoplankton absorbing and using nutrient substances and is widely applied when investigating nutrient limitations of marine phytoplankton. An nitrogen to phosphorus ratio (N:P ratio) higher than 22:1 in seawater suggests a potential P limitation, whereas a ratio below 10:1 suggests a potential N limitation [51]. Figure 4e,f show the N:P ratio of DIN:PO 4 . The surface N:P ratio in WHK was higher than 22:1, indicating an imbalanced nutrient proportion with excessive N and a potential P limitation. In SHK (SM18-MM13), the bottom N:P ratio was below 10:1, suggesting excessive P and a potential N limitation. From SM20 to the east, this ratio decreased gradually, indicating a shift from P limitation to N limitation. From April to September, the surface N:P ratio was higher than 22:1, indicating excessive nitrogen and a potential P limitation in summer. Except in May and June, the bottom N:P ratio was only slightly higher than 22:1. The surface N:P ratio was higher than 22:1 and the excessive N resulted in an antecedent consumption of P, with phytoplankton growth under P limitation. The considerable water stratification hindered the nutrient exchange between the surface and bottom, leading to a lower concentration of PO 4 at the surface. nutrient exchange between the surface and bottom, leading to a lower concentration of PO4 at the surface.

Characteristic of Dissolved Oxygen
The annual mean DO concentrations at the surface and bottom were 6.7 and 6.0 mg/L, respectively (Figure 5a). Surface DO concentration was higher in SHK than in WHK (Figure 5a). In winter, bottom DO concentration was slightly higher than surface DO concentration, whereas in summer, surface DO concentration was higher than bottom DO concentration (Figure 5b). In winter, DO concentration reached the highest values in February at 7.6 mg/L (surface) and 7.7 mg/L (bottom). Surface DO concentration was lowest in autumn at 5.9 mg/L in October, while bottom DO concentration was lowest in summer at 4.2 mg/L in July. The maximum difference between surface and bottom DO concentration was 2.9 mg/L, and the vertical stratification of DO concentration was significant in summer. In winter, surface water temperature was the lowest and surface DO concentration was the highest. Surface DO concentration also reached the highest value in winter due to the strong vertical mixing induced by the northeast monsoon. In spring, with increased solar radiation, surface water temperature increased and the dissolvability of oxygen decreased. Hence, DO concentration was lower than that in winter. In summer, surface water temperature, DIN, PO4 and Chl a reached the highest value for the year, and surface DO increased. The low-oxygen bottom waters on the continental shelf from the South China Sea moved towards the shore to compensate

Characteristic of Dissolved Oxygen
The annual mean DO concentrations at the surface and bottom were 6.7 and 6.0 mg/L, respectively ( Figure 5a). Surface DO concentration was higher in SHK than in WHK (Figure 5a). In winter, bottom DO concentration was slightly higher than surface DO concentration, whereas in summer, surface DO concentration was higher than bottom DO concentration (Figure 5b). In winter, DO concentration reached the highest values in February at 7.6 mg/L (surface) and 7.7 mg/L (bottom). Surface DO concentration was lowest in autumn at 5.9 mg/L in October, while bottom DO concentration was lowest in summer at 4.2 mg/L in July. The maximum difference between surface and bottom DO concentration was 2.9 mg/L, and the vertical stratification of DO concentration was significant in summer. In winter, surface water temperature was the lowest and surface DO concentration was the highest. Surface DO concentration also reached the highest value in winter due to the strong vertical mixing induced by the northeast monsoon. In spring, with increased solar radiation, surface water temperature increased and the dissolvability of oxygen decreased. Hence, DO concentration was lower than that in winter. In summer, surface water temperature, DIN, PO 4 and Chl a reached the highest value for the year, and surface DO increased. The low-oxygen bottom waters on the continental shelf from the South China Sea moved towards the shore to compensate for surface outflow due to the summer upwelling [42,52]. Because of the low-oxygen bottom waters, the strong stratification, and more oxygen consumption by the bottom phytoplankton debris, bottom DO concentration decreased to the lowest value. In autumn, during the monsoon transition period, with decreasing solar radiation, surface water temperature decreased gradually, while the dissolvability of oxygen increased gradually. However, because of the relatively high water temperature, DO concentration was only slightly higher than that in summer.
Water 2020, 12, 2475 7 of 18 for surface outflow due to the summer upwelling [42,52]. Because of the low-oxygen bottom waters, the strong stratification, and more oxygen consumption by the bottom phytoplankton debris, bottom DO concentration decreased to the lowest value. In autumn, during the monsoon transition period, with decreasing solar radiation, surface water temperature decreased gradually, while the dissolvability of oxygen increased gradually. However, because of the relatively high water temperature, DO concentration was only slightly higher than that in summer. In WHK and SHK, the annual means of salinity, DIN, and PO4 showed significant differences, whereas those of temperature, Chl a, and DO were similar. The monthly means of temperature, salinity, DO, and DIN showed obvious seasonal variations, whereas PO4 showed only small changes. The monthly means of surface Chl a showed obvious seasonal variations, whereas bottom Chl a did not show any significant changes. DO showed seasonal changes with temperature, salinity and DIN, and it was also influenced by the Chl a at the surface.
Differences in the monthly means of DO between the surface and bottom presented a significant negative correlation with salinity difference between the surface and bottom (correlation coefficient r = −0.98, significance level p < 0.01) and a significant positive correlation with temperature difference between the surface and bottom (r = 0.99, p < 0.01), indicating that water stratification is correlated with DO concentration.

Spatial-Temporal Mode of Dissolved Oxygen
The Shapiro filter [53] was applied to the monthly data of DO concentration and other environmental factors for surface and bottom waters. Then, REOF analysis was conducted for the anomaly field of each factor, and the Gaussian low pass filter was applied to the modes of RPCs of DO concentration.
EOFs of a space-time physical process can represent mutually orthogonal space modes where the data variance is concentrated, with the first mode being responsible for the largest part of the variance, the second for the largest part of the remaining variance, and so on [54]. Table 1 shows the variance contributions from the EOF analysis (pre-rotation) and REOF analysis (post-rotation) for environmental factors. The first two modes occupied more than 85% of the total variance. Considering their significant contributions and accumulated contribution to the total square error, these two modes could be selected to analyze the characteristics of regional environmental factors [55]. The first two EOF modes of surface DO concentration accounted for 69.2% and 18.0% of the total variance, respectively, and those of bottom DO concentration accounted for 92.4% and 3.3%, respectively. The first two REOF modes of surface DO concentration explained 45.6% and 41.6% of the total variance, respectively, and those of bottom DO concentration explained 62.3% and 33.4%, respectively. The distributions of the contribution of loading factors for REOF were more uniform than those for EOF, resulting in a more distinct structure of spatial variation. The variance contributions of the first mode for both the surface and bottom were considerably larger than those In WHK and SHK, the annual means of salinity, DIN, and PO 4 showed significant differences, whereas those of temperature, Chl a, and DO were similar. The monthly means of temperature, salinity, DO, and DIN showed obvious seasonal variations, whereas PO 4 showed only small changes. The monthly means of surface Chl a showed obvious seasonal variations, whereas bottom Chl a did not show any significant changes. DO showed seasonal changes with temperature, salinity and DIN, and it was also influenced by the Chl a at the surface.
Differences in the monthly means of DO between the surface and bottom presented a significant negative correlation with salinity difference between the surface and bottom (correlation coefficient r = −0.98, significance level p < 0.01) and a significant positive correlation with temperature difference between the surface and bottom (r = 0.99, p < 0.01), indicating that water stratification is correlated with DO concentration.

Spatial-Temporal Mode of Dissolved Oxygen
The Shapiro filter [53] was applied to the monthly data of DO concentration and other environmental factors for surface and bottom waters. Then, REOF analysis was conducted for the anomaly field of each factor, and the Gaussian low pass filter was applied to the modes of RPCs of DO concentration.
EOFs of a space-time physical process can represent mutually orthogonal space modes where the data variance is concentrated, with the first mode being responsible for the largest part of the variance, the second for the largest part of the remaining variance, and so on [54]. Table 1 shows the variance contributions from the EOF analysis (pre-rotation) and REOF analysis (post-rotation) for environmental factors. The first two modes occupied more than 85% of the total variance. Considering their significant contributions and accumulated contribution to the total square error, these two modes could be selected to analyze the characteristics of regional environmental factors [55]. The first two EOF modes of surface DO concentration accounted for 69.2% and 18.0% of the total variance, respectively, and those of bottom DO concentration accounted for 92.4% and 3.3%, respectively. The first two REOF modes of surface DO concentration explained 45.6% and 41.6% of the total variance, respectively, and those of bottom DO concentration explained 62.3% and 33.4%, respectively. The distributions of the contribution of loading factors for REOF were more uniform than those for EOF, resulting in a more distinct structure of spatial variation. The variance contributions of the first mode for both the surface and bottom were considerably larger than those of other modes. The changes indicated by the first mode represented the major characteristics of environmental factors.  Figure 6 displays the REOFs of DO. The first spatial patterns (REOF1) of surface and bottom DO exhibited positive loading factors in SHK, whereas the second spatial patterns (REOF2) exhibited positive loading factors in WHK (Figure 6), reflecting the difference in DO variation between these two areas.     Table 2 shows the correlation coefficients between the REOFs of DO and those of other environmental factors. The absolute value of correlation coefficient for surface waters ranged from 0.84 to 0.99, and that for bottom waters ranged from 0.16 to 0.87. The correlation was higher for surface waters than for bottom waters. In surface waters, the REOF1 and REOF2 of DO was similar to those of temperature and Chl a but opposite to those of salinity, DIN and PO4. The REOF1 of DO and Chl a had the highest correlation coefficient in surface waters, reaching 0.99 (p < 0.01). In bottom waters, the REOF1 of DO was opposite to those of salinity, DIN and PO4. The absolute value of the correlation coefficient between the REOF1 of DO and Chl a was the lowest in bottom waters at −0.48 (p > 0.05). In bottom waters, the REOF2 of DO was opposite to that of temperature, salinity and PO4. The absolute value of the correlation coefficient between the REOF2 of DO and DIN was the lowest in bottom waters at only 0.16 (p > 0.05).    Table 2 shows the correlation coefficients between the REOFs of DO and those of other environmental factors. The absolute value of correlation coefficient for surface waters ranged from 0.84 to 0.99, and that for bottom waters ranged from 0.16 to 0.87. The correlation was higher for surface waters than for bottom waters. In surface waters, the REOF1 and REOF2 of DO was similar to those of temperature and Chl a but opposite to those of salinity, DIN and PO 4 . The REOF1 of DO and Chl a had the highest correlation coefficient in surface waters, reaching 0.99 (p < 0.01). In bottom waters, the REOF1 of DO was opposite to those of salinity, DIN and PO 4 . The absolute value of the correlation coefficient between the REOF1 of DO and Chl a was the lowest in bottom waters at −0.48 (p > 0.05). In bottom waters, the REOF2 of DO was opposite to that of temperature, salinity and PO 4 . The absolute value of the correlation coefficient between the REOF2 of DO and DIN was the lowest in bottom waters at only 0.16 (p > 0.05).
The maximum entropy spectral analysis showed that the main periods of the first RPC (RPC1) of DO for surface waters were 24, 12, and 6 months, and the main period of the second RPC (RPC2) was 12 months (Table 3). The main period of the RPC1 and RPC2 of DO for bottom waters was 12 months. Except for the RPC1 of DO, the main period of the RPC1 for other environmental factors was 12 months.
The main period of the RPC2 for other environmental factors was 12 months or 6 and 12 months. The RPC1 of Chl a for bottom waters, the RPC2 of Chl a and PO 4 for surface and bottom waters did not significantly correspond to any period.   Harmonic constants in significant cycles of RPCs for surface and bottom DO were calculated by harmonic analysis using the least squares method. The annual amplitude of RPC1 for surface DO was the largest (1.03 mg/L), followed by the semiannual amplitude (0.66 mg/L) and the biennial amplitude (0.63 mg/L).
Except for the main periodic components, both the RPC1 and RPC2 of surface and bottom DO exhibited short-term disturbances. As we focused on the characteristics of regular seasonal changes for factors, we discussed the fitting sequence of RPCs obtained by the least squares method. Table 4 shows the correlation coefficients between the RPCs of environmental factors and their fitted values. The correlation coefficients of bottom DO were higher than those of surface DO, indicating only small disturbance of bottom DO. The correlation coefficient was the highest for temperature but limited for PO 4 (r < 0.6). Fitting data showed that the RPC1 of surface DO had a bimodal distribution, with peaks in February and June and valleys in April and October (Figure 7a). The RPC2 of surface DO had a unimodal distribution with the peak and valley in February and August, respectively (Figure 7b). The RPC1 and RPC2 of bottom DO were unimodally distributed with the peak and valley in January and July, respectively (Figure 7c  Delayed correlation analysis indicated that the RPC1 of surface DO was significantly correlated with the RPC2 of surface DO leading by 2 months (r = 0.75, p < 0.01) and also significantly correlated with the RPC1 of bottom DO leading by 3 months (r = 0.75, p < 0.01). The RPC1 and RPC2 of bottom DO showed a significant synchronous correlation (r = 1.00, p < 0.01), whereas the RPC2 of surface DO had a significant positive correlation with the RPC2 of bottom DO leading by 1 month (r = 1.00, p < 0.01). Therefore, the seasonal change between the RPC1 and RPC2 of surface DO showed a phase difference of 2 months, whereas the variability between the RPC1 and RPC2 of bottom DO were in phase. The first mode phase took approximately 3 months to propagate vertically from the surface to the bottom, whereas the second mode phase took approximately one month. The seasonal disturbance of surface DO concentration could be mainly attributed to oxygen exchange between water and air and oxygen consumption by bioactive mitogenetic effects. For bottom DO concentration, it could be mainly attributed to oxygen consumption by bioactive mitogenetic effects and sediments, and the variation of bottom waters lagged that of surface waters by 3 months.
The delayed correlation coefficients between the fitting sequences of RPCs of DO and those of other environmental factors are shown in Table 5. In surface waters, the RPC1 of DO had a significant negative correlation with the RPC1 of salinity, and a significant positive correlation with the RPC1 of temperature, Chl a, DIN and PO 4 , leading by 2~3 months. The absolute values of correlation coefficients were greater than 0.70 (p < 0.01). The RPC2 of surface DO showed a significantly negative correlation with the RPC2 of temperature and DIN, and a significantly positive correlation with the RPC2 of salinity lagging by 1~2 months; the correlation coefficients were −0.99, −0.91 and 0.90, respectively. In bottom waters, the RPC1 of DO had a significant negative correlation with the RPC1 of DIN and PO 4 with a 1 month lag or lead, and a significantly positive correlation with the RPC2 of salinity. The RPC1 and RPC2 of DO were both negatively correlated with those of temperature. The correlation coefficients between the RPC1 of DO and those of DIN and PO 4 in bottom waters were the highest, reaching −1.00 (p < 0.01). The RPC2 of DO was significantly negatively correlated with the RPC2 of salinity and DIN, lagging by 1 month, and the correlation coefficients were both −0.98 (p <0.01).  There was no significant correlation between the RPC1 of DO and that of Chl a in bottom waters, as well as between the RPC2 of DO and that of Chl a and PO 4 in surface and bottom waters, because their cycles were not significant. The REOF2 of DO was only weakly correlated with the REOF2 of DIN in bottom waters. Therefore, considering the correlation of REOFs and RPCs between DO and other environmental factors, in SHK, the spatial-temporal variations in surface and bottom DO were highly related to the variations of salinity, DIN, and PO 4 , and the variation of surface DO was also connected to the variation of temperature and Chl a. In WHK, the variations in surface and bottom DO were highly related to the variations of salinity and temperature, and the variation of surface DO was also connected to the variation of DIN. Thus, the variation of DO was highly related to the variation of salinity in the southeast of PRE. The influence of Chl a and PO 4 on DO was stronger in SHK than those in WHK.
The variability of salinity is obviously affected by river discharges and winds [19,56]. Figure 8 shows that the river discharge and alongshore wind speed varied greatly over the years but higher river discharge often coincided with a predominant southwesterly wind during the wet season in PRE. To examine the effects of river discharges and wind on the variation of DO, we calculated the correlations between the fitted RPCs of DO and the river discharge and the alongshore wind speed at Nei Lak Shan ( Table 6). The RPCs of DO were significantly correlated with the river discharge and the alongshore wind speed, whereas the correlation between the RPC1 of surface DO and river discharge and the correlation between the RPC2 of surface DO and alongshore wind speed were relatively weaker ( Table 6). The RPC1 of surface DO was positively correlated with river discharge and alongshore wind speed, and the correlation between the RPC1 of surface DO and river discharge was greater than the correlation between the RPC1 of surface DO and alongshore wind speed ( Table 6). In contrast, the RPC2 of surface DO was negatively correlated with river discharge and alongshore wind speed, and the correlation between the RPC1 of surface DO and alongshore wind speed was greater than the correlation between the RPC1 of surface DO and river discharge (Table 6). Thus, the temporal variability of surface DO in WHK was more affected by river discharge than that in SHK, whereas the temporal variability of DO in SHK was more affected by alongshore wind speed than that in WHK. The RPCs of bottom DO were negatively correlated with river discharge and alongshore wind speed. The correlation between the RPCs of bottom DO and river discharge was greater than the correlation between the RPCs of bottom DO and alongshore wind speed. The temporal variability of bottom DO was more affected by river discharge than by alongshore wind speed.
Water 2020, 12, 2475 12 of 18 river discharge often coincided with a predominant southwesterly wind during the wet season in PRE. To examine the effects of river discharges and wind on the variation of DO, we calculated the correlations between the fitted RPCs of DO and the river discharge and the alongshore wind speed at Nei Lak Shan ( Table 6). The RPCs of DO were significantly correlated with the river discharge and the alongshore wind speed, whereas the correlation between the RPC1 of surface DO and river discharge and the correlation between the RPC2 of surface DO and alongshore wind speed were relatively weaker ( Table 6). The RPC1 of surface DO was positively correlated with river discharge and alongshore wind speed, and the correlation between the RPC1 of surface DO and river discharge was greater than the correlation between the RPC1 of surface DO and alongshore wind speed ( Table 6). In contrast, the RPC2 of surface DO was negatively correlated with river discharge and alongshore wind speed, and the correlation between the RPC1 of surface DO and alongshore wind speed was greater than the correlation between the RPC1 of surface DO and river discharge (Table 6). Thus, the temporal variability of surface DO in WHK was more affected by river discharge than that in SHK, whereas the temporal variability of DO in SHK was more affected by alongshore wind speed than that in WHK. The RPCs of bottom DO were negatively correlated with river discharge and alongshore wind speed. The correlation between the RPCs of bottom DO and river discharge was greater than the correlation between the RPCs of bottom DO and alongshore wind speed. The temporal variability of bottom DO was more affected by river discharge than by alongshore wind speed.

Bottom Hypoxia
From 2000 to 2017, bottom hypoxia was observed in the summers of 7 years at 11 stations in the southeast of PRE, with the maximum in August (Figure 9a, Table 7). Except at station SM20, the maximum difference of DO between surface and bottom waters was 7.7 mg/L, whereas the minimum difference was 4.4 mg/L (Table 7). Regarding temperature, the maximum difference between surface and bottom waters was 6.7 • C, and the minimum difference was 3.4 • C. The maximum difference of salinity between surface and bottom waters was −24.1, and the minimum difference was −4.4. Thus, the water stratification was obvious. There was no obvious stratification at station SM20 may be due to the shallow water depth (7 m). Hypoxia occurred twice in WHK and nine times in SHK, which corresponded to 18.2% and 81.8% of the total occurrences, respectively (Figure 9b). Hypoxia was observed the most frequently at station SM18 with up to four events, accounting for a proportion of 36.4%, and the lowest concentration of DO was 0.40 mg/L in bottom waters. Recent research has found that there was a hypoxic zone near station SM18, which occurred in summer of 2011 [22]. At all hypoxic stations, the N:P ratios of surface waters were beyond 22:1, indicating excessive nitrogen levels and a potential P limitation. This suggested that large amounts of nutrients were consumed by phytoplankton and the PO 4 in surface waters was used up. After the plankton died, organic debris sank to the bottom and degraded, resulting in the generation of PO 4 , which consumed the DO at the bottom [25,57,58]. The strong water stratification prevented the vertical exchange of oxygen between bottom and surface waters, resulting in hypoxia. At station SM20, hypoxia occurred in August 2010, but the concentration of surface DO was 2.50 mg/L, which may be attributable to the close proximity of the station to the shore and strong vertical mixing.

Bottom Hypoxia
From 2000 to 2017, bottom hypoxia was observed in the summers of 7 years at 11 stations in the southeast of PRE, with the maximum in August (Figure 9a, Table 7). Except at station SM20, the maximum difference of DO between surface and bottom waters was 7.7 mg/L, whereas the minimum difference was 4.4 mg/L (Table 7). Regarding temperature, the maximum difference between surface and bottom waters was 6.7 °C, and the minimum difference was 3.4 °C. The maximum difference of salinity between surface and bottom waters was −24.1, and the minimum difference was −4.4. Thus, the water stratification was obvious. There was no obvious stratification at station SM20 may be due to the shallow water depth (7 m). Hypoxia occurred twice in WHK and nine times in SHK, which corresponded to 18.2% and 81.8% of the total occurrences, respectively (Figure 9b). Hypoxia was observed the most frequently at station SM18 with up to four events, accounting for a proportion of 36.4%, and the lowest concentration of DO was 0.40 mg/L in bottom waters. Recent research has found that there was a hypoxic zone near station SM18, which occurred in summer of 2011 [22]. At all hypoxic stations, the N:P ratios of surface waters were beyond 22:1, indicating excessive nitrogen levels and a potential P limitation. This suggested that large amounts of nutrients were consumed by phytoplankton and the PO4 in surface waters was used up. After the plankton died, organic debris sank to the bottom and degraded, resulting in the generation of PO4, which consumed the DO at the bottom [25,57,58]. The strong water stratification prevented the vertical exchange of oxygen between bottom and surface waters, resulting in hypoxia. At station SM20, hypoxia occurred in August 2010, but the concentration of surface DO was 2.50 mg/L, which may be attributable to the close proximity of the station to the shore and strong vertical mixing.  Based on the existing studies, the entire pattern in PRE between the DO and various factors (temperature, salinity, nutrients, and Chl a) did not exhibit a simple linear relationship for the complex biochemical and physical processes [22]. The variability of DO was regulated by different factors in SHK and WHK. In SHK, surface DO was always oversaturated in summer, because of the high transparency, the phytoplankton bloom driven by nutrients carried by river fresh water [22,59,60]. Bottom DO was always low in summer, because of the strong stratification, more oxygen consumption by the bottom phytoplankton debris and the low-oxygen bottom waters from the South China Sea due to upwelling. This was in agreement with the results of previous studies [22,60,61]. Bottom waters experienced some seasonal hypoxia. In WHK, due to the greater turbidity, low residence time, and shallow topography, hypoxia rarely occurred. SHK appeared to be more vulnerable to hypoxia than WHK.

Conclusions
The spatial-temporal characteristics of seasonal variations in DO concentration in the southeast of PRE were analyzed using multiple statistical analysis methods for monthly data covering 18 years. Surface DO concentration was higher in SHK than in WHK. DO stratification was significant in summer, but it was not distinct in winter, during which DO concentration was the highest. DO stratification showed a significantly positive correlation with water stratification.
DO concentration fields in surface and bottom waters in SHK and WHK exhibited distinct characteristics of seasonal change in the recent 18 years. The main periods of the surface DO variation were 24, 12, and 6 months in SHK, while main period was 12 months in WHK. The main period of the bottom DO variation was 12 months in the southeast of PRE. In SHK, the spatial-temporal variations in surface and bottom DO were highly related to the variations of salinity, DIN, and PO 4 , and the variation of surface DO was also connected to the variation of temperature and Chl a. In WHK, the variations in surface and bottom DO were highly related to the variations of salinity and temperature, and the variation of surface DO was also connected to the variation of DIN. The temporal variability of surface DO in WHK and SHK was greatly affected by the river discharge and the alongshore wind speed, respectively. However, the river discharge had a greater influence than the alongshore wind speed on the temporal variability of bottom DO in WHK and SHK. These findings demonstrated that coupled physical and biochemical processes were important in driving the variations of DO in the southeast of PRE.
From 2000 to 2017, bottom hypoxia was observed in the southeast of PRE in the summers of 7 years. SHK appeared to be more vulnerable to hypoxia than WHK. A potential P limitation appeared in surface waters at stations experiencing hypoxia. Bottom hypoxia in the estuary could be mainly attributed to oxygen consumption by organic debris and water stratification.