In this study, the Regional Ocean Modeling System (ROMS) (Shchepetkin and Mcwilliams, 2005) was configured with the model domain covering the Pearl River Estuary and its adjacent shelf (Figure 1A). The horizontal grid resolution was 500–800 m in the Pearl River Estuary and increased to 3–6 km in the open ocean. The eastern and western boundaries were perpendicular to the isobaths, while the southern boundary was located near the shelf break. The observed water depth data provided by Southern Marine Science and Engineering Guangdong Laboratory (Zhuhai) combined with the ETOPO2 data were interpolated to the model grids. A terrain-following S-coordinate (Song and Haidvogel, 1994) with 20 vertical levels was configured, and the surface and bottom layers were refined. A third-order HSIMT advection scheme developed by Wu and Zhu (2010) was selected for tracer advection simulation. Mellor and Yamada’s level 2.5 turbulent closure scheme (Mellor and Yamada, 1982) and Smagorinsky scheme (Smagorinsky, 1963) were used to parameterize vertical mixing and horizontal mixing. The model was initialized with climatologic temperature and salinity from the Ocean General Circulation Model for the Earth Simulator (OFES) data. Open boundary conditions of temperature, salinity and shelf current were also obtained from OFES data. Thirteen tide constituents (M₂, S₂, N₂, K₂, K₁, O₁, P₁, Q₁, MF, MM, MN₄, M₄, and MS₄) extracted from Oregon State University (OSU) global inverse tidal model of TPXO8 (Egbert and Erofeeva, 2002) were used to provide tidal current and elevation for the open boundary. According to the diversion ratio suggested by Lu et al. (2013), the Pearl River runoff was distributed at eight major outlets. The atmospheric parameters (wind, air temperature, relative humidity, surface air pressure, cloud cover and net shortwave radiation) were provided by the European Centre for Medium-Range Weather Forecast (ECMWF) products (ERA-5 for realistic experiment and ERA-Interim for climatological experiment). The climatological experiment was driven by climatological monthly mean data of runoff, open boundary conditions, and atmospheric parameters. It was spun up for one year, and results of the second year were output for analysis. The configuration of the realistic experiment for model validation was the same as that of the climatological model except it was driven by 6-h atmospheric parameters.

Observed elevation and salinity data were used to validate the model. The elevation data used for comparison was observed during January 2008 at tide gauge stations that are marked with green triangles in Figure 1B. The salinity data were measured by a cruise organized by the Southern Marine Science and Engineering Guangdong Laboratory (Zhuhai). The cruise was conducted in the Pearl River Estuary and the adjacent shelf region from July 6 to 13, 2020. The sampling stations (as shown in Figure 1B) spanned from the Pearl River Estuary to the Huangmaohai and reached offshore to ∼50-m isobath.

Tides in the Pearl River Estuary and the adjacent water feature irregular semi-diurnal characteristics. K₁, O₁, M₂ and S₂ are the dominant tidal constituents in the NSCS, of which M₂ is the strongest. Harmonic analysis was conducted for observed and modeled elevation data to obtain the amplitude and phase of different tidal constituents. The vectorial difference (Foreman et al., 1993) was then used to compare the modeled results with the observation:

in which a_0 and g_0 are the observed harmonic tidal amplitude and phase, and a_m and g_m are the modeled counterparts. To normalize the results, the vectorial difference was further divided by the amplitude:

The results of vectorial difference are listed in Table 1. The mean Diff of M₂ of all stations was 8.08 cm, and the mean Rdiff was 15.91%. The model simulated the M₂ constituent reasonably. The Diff of S₂ and K₁ were 5.73 cm and 7.35 cm, respectively, and the associated Rdiff did not exceed 25%. Rdiff of O₁ was 26.53% higher than other tidal constituents, but its absolute value was small given the small O₁ amplitude. Among the six tide gauge stations, the simulation at Shanbanzhou, located at the upper estuary inside the Humen, showed the largest error. The Rdiff of three harmonic tides (O₁, M₂ and S₂) exceeded 30% at this tide station. The tide in the river mouth is deformed to a certain extent, increasing the difficulty in accurate simulations. Given that the focus of this study was on the plume region, we did not further improve the validations. The overall model accuracy of simulating the tide was comparable with previous studies (e.g., Zu et al., 2014).

Observed data during July 2020 showed there was a low-salinity tongue outside the Pearl River Estuary (Figure 3A) and some researchers signify that the phenomenon is correlated to the West Channel geometry, estuarine outflow, tide, and wind (Pan et al., 2014; Liu et al., 2020). The surface plume extended in both directions to the east and west simultaneously (Figure 3A). The output data from the realistic forcing was compared with the observation (Figures 3B,D), showing reasonable skill in simulating the structure of the Pearl River plume both at the surface and the bottom. Inside the Pearl River Estuary, the shelf water intrusion simulated by the model was slightly weaker. Also, the modeled salinity near the Huangmaohai was relatively lower than the observation (Figures 3C,D). The reason could be that the model’s runoff releasing position was too close to the ocean, resulting in insufficient mixing before leaving the estuary. Other than these two areas, the modeled salinity distribution was close to the observed data. Model performance was quantified via the correlation coefficient (CC), the root-mean-square error (RMSE), and the skill score (SS). These metrics were calculated as follows:

in which X is the variable of interest and X¯ is its mean. Model performance was categorized as excellent when SS ≥ 0.65, very good when 0.65 > SS ≥ 0.5, good when 0.5 > SS ≥ 0.2, and poor when SS < 0.2 (Maréchal, 2004; Allen et al., 2007). Although RMSE was greater than 3, the CC exceeded 0.9 and the SS was 0.8 in both the surface and bottom layers (Figures 3E,F). The relatively large RMSE but excellent CC and SS were caused by the strong spatial variability of salinity in the Pearl River Estuary. Therefore, the model was reliable and can be used for further investigation.

According to previous studies (Ou et al., 2009; Zu et al., 2014), the 32-psu isohaline can roughly denote the horizontal boundary of the Pearl River plume. We used the 75-h tidal-averaged salinity during the spring tides in different seasons to represent the seasonal variation of salinity distributions (Figure 4). In summer, as the runoff increased, the freshwater gradually extended outside the river mouth, and the prevailing southwesterly wind drove the coastal current to the northeast. The plume was thus brought to the eastern side of Pearl River Estuary. However, a portion of the plume extended toward the west coast simultaneously. The bottom plume, however, did not turn as did the surface. The bottom plume front existed only on the west coast of the Pearl River Estuary in summer, indicating that the plume extension to the east was surface-advected and that the extension to the west was bottom-trapped. The plume completely extended to the west in winter when the wind turned northeasterly. Its width was narrower than in summer. There were also seasonal variations inside the Pearl River Estuary. In summer, the surface plume filled the entire estuary, and the bottom shelf water intruded landward along two deep channels. In winter, the plume accumulated to the western side of the Pearl River Estuary, and the bottom shelf water intrusion that occupied the east side of the estuary was even more robust.

Whitney and Garvine (2005) suggested a wind index Ws≡UwindUbuoyancy, in which U_wind≈2.65×10⁻²U, Ubuoyancy=1K(2gr′Qf)1/4, U is along-shelf wind speed (up-shelf direction is positive in this study), K is internal Kelvin number, gr′=g(ρ0-ρr)/ρ0 is the reduced gravity, ρ₀ is ambient seawater density, ρ_r is the freshwater density, Q is the runoff, and f is the Coriolis parameter. W_s quantifies the relative importance of wind stress and the down-shelf buoyancy-driven current. |W_s| less than one indicates that the buoyancy-driven current is stronger than the wind-driven current. The calculated |W_s| was ∼0.3 in summer but > 1 in winter. The existence of plume water along the west coast and the fact |W_s| < 1 in summer indicated that the extension of the Pearl River plume was controlled not only by the wind, but also by its baroclinicity.

The profile of Section A (Figures 5A–D) indicates that the river plume was bottom-trapped to the west of the Pearl River Estuary in summer and winter. The isohaline in winter was much sparser and more vertical than in summer. The isohaline during the spring tide was more vertical than during the neap tide. There was strong stratification in Section B in summer with a significant surface-to-bottom salinity difference. In winter, the isohaline in Section B was nearly vertical and the stratification was weakened significantly. In summer, the plume extended to the east coast under the upwelling-favorable summer monsoon (Figures 5I,J). A brackish water mass appeared in the offshore region as influenced by the surface-advected plume. In winter, Section C was occupied by high salinity water.

It is evident that the spatiotemporal variations of the Pearl River plume were very drastic. The surface-to-bottom salinity difference was also considerable in summer (Figure 5; later further shown in Figure 13A). Therefore, it is important to choose an appropriate method to describe the variability of the plume.

The gradient method is often used to identify fronts (Huang et al., 2010; He et al., 2016; Wu and Wu, 2018). In this study, the salinity gradient magnitude (SGM) was used to calculate the gradient in the surface and bottom layers of the Pearl River Estuary and its adjacent waters:

In which ΔS₁ and ΔS₂ are the diagonal salinity difference of a 3 × 3 grid, while ΔX₁ and ΔX₂ are the diagonal distance between corners of the corresponding grids. The hourly modeled salinity was used to calculate SGM and then averaged over the period of interest. The seasonal and spring-neap cycles of the front are the main focus of this study. Annual-mean salinity gradient was used to indicate the average position of the front (Figures 6A,B). The surface salinity gradient in most areas of the model domain was small except in the Modaomen, Huangmaohai and Pearl River Estuary (Figure 6A). The bottom salinity gradient in the Pearl River Estuary and its adjacent coastal water was concentrated in a narrow region shoreward of ∼10-m isobath with a value up to 2 psu/km (Figure 6B). The distribution of the high bottom salinity gradient basically followed along the coastline. The standard deviation of the surface salinity gradient was large because the surface plume was affected by wind stress. Whereas in the bottom, the standard deviation of the salinity gradient was small, indicating the bottom plume front was more stable.

Summer (June to August) and winter (December to the next February) are flood and dry seasons, respectively. Therefore, these two seasons were selected as typical seasons to explore the cross-shelf variation in salinity gradients in Sections A, B, and C (Figure 7). The bottom front did not exist in Section C but occurred in Section A and Section B in winter and summer. Both the surface front and the bottom front had variations in intensity and spatial position. The bottom salinity gradients in Section A and Section B were stronger in summer than in winter on the seasonal scale. The surface and bottom front extended seaward in summer but shrank shoreward in winter. The seasonal cross-shelf movement of the surface front has a migration scale of ∼2 km in Section A and ∼20 km in Section B. The seasonal cross-shelf movement of the bottom front was smaller than that of the surface one, with Section A being ∼1.2 km and Section B being ∼10 km. In spring-neap tide cycles, the cross-shelf movement of the surface front was not obvious in Section B, although the phenomenon was apparent in Section A. The bottom front had a stronger intensity during the neap tides than during the spring. This was due to the increase of the tide mixing in the spring tide inducing more vigorous shear dispersion that decreased the horizontal gradient. There were three extreme values on the bottom front at Section B and their positions correspond to the bathymetric variations. The bottom front moved seaward during the spring tides and landward during the neap tides. The migration scale of the bottom front in neap-spring cycles in Section A was ∼0.8 km in summer and ∼0.4 km in winter. The migration scale in neap-spring cycles in Section B was relatively larger than in Section A with ∼3 km in summer and ∼1 km in winter, respectively. The migration of the front depends not only on the strength of the external forces but also on the seafloor slope. According to the Yankovsky and Chapman (1997) theory, under the same external force, the greater the seafloor slope the smaller the horizontal migration distance of the front.

The annual-mean salinity gradient and its standard deviation (Figure 6) provided basic information of the Pearl River plume front. However, averaging over a long-time span may obscure the actual intensity of the front. Therefore, we further calculated the occurrence probability of the plume front when its intensity exceeded a given threshold. Fedorov threshold (Fedorov, 1986)-i.e., ten times the average salinity gradient in the study area- was often used to detect the front (Fedorov, 1986). However, it should be noted that the Fedorov threshold is sensitive to the domain selected. To reflect the position of the front in different seasons more reasonably, we choose 0.8 psu/km as the threshold according to Figure 7. Counting the time of which the salinity gradient is greater than 0.8psu/km in the four seasons, and then dividing it by the total time span, we got the occurrence probability of the salinity front (Figure 8) in four seasons (i.e., spring, summer, autumn and winter). In summer, the surface front existed in the east coast and moved offshore significantly; hence, there was a large area with a noticeable occurrence probability. The surface front also existed in the offshore region in spring and autumn; however, the probability dropped to ∼10%. In winter, the surface salinity front entirely extended to the west coast. In this period, the coastal extension of the front is the shortest, shrinking to the position around 113°E.

Whereas for the bottom front, the high occurrence probability area was distributed in a more narrow region and was nearly identical throughout four seasons (Figure 8). This indicated that the bottom front was strongly “rooted” on the seafloor. The occurrence of the surface front in the west region opposite to the wind direction in summer could be closely linked to the formation of the bottom front, as indicated by previous studies (e.g., Wu and Wu, 2018).

To infer the formation dynamics of such a stable bottom front, the variability of the bottom salinity gradient was further analyzed by the EOF method. Before doing the EOF, the bottom salinity gradient was low-pass filtered (with a cut-off window of 72 h) to remove the strong intra-tidal signal, and then the long-term mean was further subtracted from the data. Hence, EOF analysis was performed on the anomaly field of the bottom salinity gradient. The results are shown in Figure 9. The three leading EOF modes contributed 58.7%, 13.9% and 7.8% of the total variance, respectively, totaling 80.4%, which represented most of the variance in the bottom salinity gradient in the Pearl River Estuary.

Mode 1 showed a notable seasonal variation of the bottom front. During the spring and winter seasons, the bottom front shrank toward the coast, but it moved toward the sea during the summer and autumn seasons. This mode was mainly related to runoff and along-shelf wind with a correlation coefficient of -0.95 and -0.87, respectively. In summer and autumn, the runoff was large; meanwhile, the bottom front moved toward the sea. Correspondingly, the runoff was small in spring and winter, and the bottom front moved toward the land.

Mode 2 indicated a significant spring-neap cross-shelf movement. The correlation coefficient between Mode 2 and the tidal range was -0.57. The bottom front moved seaward during the spring tide, whereas it moved shoreward during the neap tides. It can also be seen that the influence of tide on the bottom front was the strongest in August and September, in consistence with the amplification of the tidal range around the autumnal equinox. West of ∼113.5°E, the spatial mode was different from that on the east side because there were freshwater sources from other tributaries of the estuary.

Mode 3, which contributed to the total variance less significantly, peaked during September and October. Its correlations with wind, runoff or tidal forcing were not as apparent as previous two modes. We inferred it was caused by the transition of the summer plume to the winter plume. In summer, a large area of river plume was surface advected due to the large river discharge and the upwelling-favorable summer monsoon. Then in autumn, the remainder of this plume was transported to the west coast, which supplied excess plume water that favored the formation of a strong bottom front.

Momentum terms from the model were used to diagnose the dynamic characteristics of the Pearl River plume. The average momentum terms over 75 h were used to remove tidal oscillations. The coordinate was rotated to represent the along-shelf (x) and cross-shore (y) directions.

The tidal-averaged momentum equation across the shelf reads:

in which the tidal-averaging operator was omitted. The terms on the right-hand side are the advection, Coriolis force, barotropic pressure gradient force, baroclinic pressure gradient force, and vertical friction force, respectively. The horizontal diffusion term was neglected because it was much smaller than the vertical diffusion.

Results across Section A during the spring tide in summer and winter seasons are shown in Figure 10. The barotropic pressure gradient force was seaward (positive value) and the baroclinic pressure gradient force was landward (negative value). The total pressure gradient force and Coriolis force were largely balanced in the cross-shelf direction. This momentum balance in winter was more in line with the thermal wind balance than in summer. The advection term also has larger values on the surface layers and at the top of the bottom boundary layer, indicating strong nonlinearity. The vertical friction of the surface layer was strong in summer, and the friction of the bottom boundary layer was strong in winter. Overall, the momentum balance was more geostrophic in the winter season, while the nonlinear momentum advection and friction became important in summer. This was apparently related to the different front shapes in these two seasons. Nevertheless, from Figure 10, it is clear that the locations associated with maximum bottom salinity gradient were similar in winter and summer, regardless of contrast dynamic balances. This seems consistent with the dynamics of EOF mode 2. The cross-shelf location was controlled by tidal forcing, which remains largely unchanged during different seasons.

The EOF analysis, as well as the above momentum balance analysis, indicated that tide was a dominant forcing determining the cross-shelf position of the bottom front. This, however, was in contrast to the previous conclusion that was derived from a non-tidal river plume (Chapman and Lentz, 1994; Yankovsky and Chapman, 1997). In these studies, the cross-shelf position was controlled by river discharge and wind through the bottom Ekman dynamics.

Hence, a sensitivity experiment with tide excluded was configured to further explore the role of the tide in controlling the bottom front in the Pearl River Estuary. The average surface salinity in April and the entire spring season are shown in Figure 11. The model results showed that dramatic differences can be found with tidal forcing removed. In the absence of tide, the Pearl River plume also extended to the upstream, in consistence with the findings in other estuaries (e.g., Guo and Valle-Levinson, 2007; Wu et al., 2011). With tide included, the surface plume extended further to the sea.

The difference between tidal and non-tidal simulations in all four seasons is shown in Figure 12. The bottom 32-psu contour was roughly located along the line connecting Macao and Hong Kong in the tide’s presence. In the absence of tide, a large amount of shelf water intruded into the Pearl River Estuary, hence, the bottom front moved shoreward. This was due to the enhanced estuarine circulation under weak mixing, thus strong stratification (Gong et al., 2018). The difference in the 32-psu contours in the shelf region was small between the two sensitivity experiments, indicating that the impact of the tide on the front weakened in the offshore region.

The spring-neap migration of the bottom front implied that it could be linked to the tidal front due to differential mixing over topography (Simpson and Hunter, 1974; Wu and Wu, 2018). Simpson-Hunter (SH) number can be used to compare the relative strength of the mixing caused by tidal current and by stratification (Simpson and Hunter, 1974),

in which h is the water depth, U_T is the maximum tidal velocity during the spring tide. The smaller the SH number is, the more likely the tidal mixing can overcome vertical stratification. An SH number between 1.7 and 2 was commonly taken as the position of the tidal-induced mixing front (Beardsley et al., 1985; Loder and Greenberg, 1986; Wu and Wu, 2018). As shown in Figure 13B, the SH number in the Pearl River Estuary was less than 2, indicating that tide could thoroughly mix the entire water column. In the coastal area from the entrance of the Pearl River Estuary to the west coast, SH number value was between 1.7 and 2. Thus, there was a boundary between stratified and mixed waters at the entrance of the Pearl River Estuary. This place was consistent with the region with a high occurrence probability of the bottom front (Figure 7). This indicated that the formation and maintenance of the bottom plume front were closely linked to the tidal mixing, which was similar to the Changjiang River plume (Wu and Wu, 2018).