Characteristics and Driving Mechanisms of Salinity Stratification during the Wet Season in the Pearl River Estuary, China

Abstract

In an estuary, stratification processes play a major role in inhibiting estuarine circulation, sediment transport, and the estuarine ecosystem. A detailed examination of the salinity stratification through the gradient Richardson number and the potential energy anomaly equation has been undertaken along the West Channel of the Pearl River Estuary, China. The results show that the estuarine circulation within the West Channel is much weaker on a spring tide than that on a neap tide, exhibiting apparent spring–neap tidal variability. The calculated gradient Richardson number displays its intratidal and spring–neap tidal variability within the West Channel, indicating the existence of intratidal and spring–neap tidal variability of stratification. In addition, the tidally averaged change rate of total potential energy anomaly within the West Channel suggests more than a 4.53 × 10⁻³ W·m⁻³ increase from spring to neap tides, demonstrating strong stratification on a neap tide. The longitudinal advection and the longitudinal depth-mean straining are the leading physical mechanisms contributing to intratidal and spring–neap variability of salinity stratification within the West Channel. However, the effects of the lateral terms cannot be ignored especially on a neap tide.

1. Introduction

In estuarine waters, stratification plays a unique and fundamental part in altering estuarine circulation [1,2], sediment transport [3,4] and bottom hypoxia [5,6,7,8]. Thus, understanding the physical processes underlying salinity stratification is of great importance to estuarine science.

Owing to the changes in tidal, riverine dynamics and estuarine bathymetry, stratification may exhibit remarkable temporal and spatial variations [8,9]. Several metrics have been suggested by previous works to quantify estuarine stratification, such as Froude number, stratification number, gradient Richardson number, and estuarine Richardson number, etc. [10,11]. As proposed by Simpson [12,13], the potential energy anomaly is one of the widely used metrics. Based on the potential energy anomaly, Simpson et al. [13] quantified the estuarine stratification and highlighted the important role of tidal straining. This method has been used successfully in numerous estuarine research in recent years [14,15,16,17,18,19].

In most estuaries, an increase in river discharge results in a variation in vertical density that intensifies the stratification, whereas the interaction between tides and river discharge affects mixing [9,18]. Generally speaking, the strengthened mixing during the flood tide destructs the stratification, while the weakened mixing and tidal straining make the water column of the estuary well stratified during the ebb tide [14,20]. The impacts of longitudinal straining on estuarine circulation and stratification have been widely established [21,22,23,24]. The effect of lateral dynamics on stratification has been recorded and studied in recent works [25,26,27].

Located in the central part of the Guangdong–Hong Kong–Macao Greater Bay Area, the Pearl River Estuary (PRE) is a typical well-developed estuary in China. Previous research in the PRE mainly concentrated on the vertical and longitudinal estuarine dynamic processes, whereas the effect of lateral dynamics was mostly neglected (e.g., You [28]). Numerical studies have uncovered that the stratification was highly strengthened over the lower reach of the estuary [29,30]. However, the variation in 3D hydrodynamics as well as their influence on stratification is barely known. In this paper, a well-validated numerical model is applied to quantify the stratification variations in the PRE and to investigate the relative importance of advection, straining, and mixing on salinity stratification by using the potential energy anomaly.

The remainder of this paper is structured as follows. Section 2 introduces the study area and the model simulation implementation. In Section 3, the model results are exhibited including the descriptions of the velocity and salinity structures, and stratification alterations on inter- and intratidal time scales. Section 4 discusses the principal physical mechanisms controlling estuarine salinity stratification. Finally, conclusions are summarized in Section 5.

2. Study Area and Methodology

2.1. Study Area

The Pearl River Delta is a complicated river network system, with three major rivers (i.e., North River, West River, and East River) and eight outlets (i.e., Humen, Jiaomen, Hongqimen, Hengmen, Modaomen, Jitimen, Hutiaomen, and Yamen, Figure 1). The Pearl River Estuary (PRE) is a funnel-shaped estuary with a width of approximately 21 km at the mouth and a width of about 4 km at the head. The water depth in the PRE is shallower than 5 m, except for the two main channels (i.e., the West and East Channels) with depths of 5–20 m (Figure 1).

The PRE receives an annual average freshwater runoff of 3 × 10¹¹ m³ from the Pearl River, with strong seasonal alterations from ~3400 to ~21,000 m³·s⁻¹ [31]. The PRE is characterized by a micro-tidal regime with an average tidal range of about 1.0 m and the M₂ semi-diurnal tide is the dominant tidal constituent followed by K₁, O₁, and S₂ [32,33]. The tidal range increases in the landward direction with the maximum mean tidal range of 1.7 m at the Humen outlet near the estuary head. Under the combined effect of river discharge and the tides, the tidal currents in the PRE are relatively strong with a maximum current velocity up to 1.5 m·s⁻¹ [32]. The estuary is a partially mixed estuary in the wet season; however, it may change into a highly-stratified or well-mixed estuary relying on river runoff and tidal ranges.

2.2. Numerical Modelling

To simulate the alterations of the complicated hydrodynamics accurately in the Pearl River Delta, the Semi-implicit Cross-scale Hydroscience Integrated System Model (SCHISM) [34,35] is adopted in this study. The hydrodynamic model has a large model domain with finer grids in the inner delta to well resolve the complex bathymetry and geometry of the channel networks.

A mesh of unstructured triangular elements is generated with the total mesh element number being 56,911. The computational domain and the unstructured triangle mesh of the model are shown in Figure 2. The southern open ocean boundary extends to about 30 m water depth. The topological data include the field data of the river networks in the 2010s, and nautical charts of the PRE from 2015 published by the Maritime Safety Administration of China. Then, the raw bathymetry data are interpolated and gridded to each mesh point with the model grid size ranging between 20 m and 3 km. Twelve σ layers are set in the vertical with the vertical grid size ranging from 1 m (within the PRE) up to approximately 3 m (outside the PRE).

The upstream riverine boundary is driven by the measured daily river discharges at the Gaoyao, Shijiao, Boluo, Shizui, and Laoyagang Stations (positions shown in Figure 1 and Figure 2). Figure 3 indicates the annual cycle of river flows at three major stations and the river flows during the specific modeled period. The open ocean boundary is driven by the predicted hourly tidal levels derived from eight harmonic tidal constituents (M₂, S₂, N₂, K₂, O₁, K₁, P₁, Q₁) taken from the TPXO8 global model of ocean tides (http://volkov.oce.orst.edu/tides/tpxo8_atlas.html; Accessed on 28 August 2019; [36]) with a resolution of 1/30°. The salinity data at the open ocean boundary are downloaded from the HYbrid Coordinate Ocean Model with a resolution of 1/12° (http://www.hycom.org; Accessed on 3 November 2020). The roughness ranges from 0.008 to 0.037, with an increasing tendency from the estuary to the river networks. The model ran for August 2007 for validation. After validation, a 57-day simulation is implemented during the wet season from 15 May 2010 to 10 July 2010 for stratification analysis.

2.3. Evaluations of Stratification

2.3.1. The Gradient Richardson Number

The Brunt–Väisälä frequency (also the buoyancy frequency), N, is often utilized to present the intensity of stratification [23,37]. It is denoted using:

$$N^2=-\frac{g}{{\rho }_0}\frac{\partial \rho }{\partial z}$$

(1)

where g indicates the gravitational acceleration (9.8 m s⁻²), ρ₀ is the reference density (1027 kg·m⁻³), and ρ means the density in each water column layer. The smaller the buoyancy frequency, the smaller the density gradient, and the weaker the water column stratification, and vice versa.

In a tidal estuary, the alteration of density with pressure may be ignored, and quite often the water density is much more affected by salinity than temperature [1,38]. In this case, the estuarine water density ρ could be assessed using the state equation with temperature neglected [39,40]:

$$\rho ={\rho }_0\left(1+\beta S_w\right)$$

(2)

where β denotes the saline contraction coefficient (7.7 × 10⁻⁴; [41]), and Sw represents salinity of the water.

Defined as the ratio between the gravitational buoyancy force and the turbulent shear production, the dimensionless gradient Richardson number (Ri) can be used to evaluate the effect of stratification on the turbulent mixing. The expression of Ri can be written as [37]:

$$Ri=\frac{N^2}{S^2}$$

(3)

Assuming that the horizontal velocity gradient is much smaller than the vertical velocity gradient, the mean shear (S) can be defined as:

$$S^2=(\frac{\partial u}{\partial z}{)}^2+(\frac{\partial v}{\partial z}{)}^2$$

(4)

where the ∂u/∂z and ∂υ/∂z on the right-hand side are the vertical shears in along-channel and lateral flow.

According to the linear stability theory, the critical value of Ri for active mixing is 0.25, because this is the necessary condition for the shear instability to occur in estuaries [42,43]. Therefore, the modeled gradient Richardson numbers are illustrated in the form of log₁₀(Ri/0.25) in the present study.

2.3.2. The Potential Energy Anomaly

Defined as the energy required to overturn the water column to a well-mixed layer, the potential energy anomaly (ϕ) is applied as another index of the strength of stratification [12]:

$$\varphi =\frac{\mathrm{g}}{D}{\displaystyle {\int }_{-H}^{\eta }\left(\overline{\rho }-\rho \right)}zdz$$

(5)

where z denotes the vertical coordinate from the water surface at z = η to the sea bottom at z = −H, D means the total water depth (D = η + H), and $\overline{\rho }$ represents the mean density of the water column. Based on the salinity and temperature budget equations, Burchard and Hofmeister [14] derived the potential energy anomaly budget equations analytically, which have been widely applied in estuarine stratification studies, and are:

$$\begin{array}{c}\hfill \frac{\partial \varphi }{\partial t}=-{\nabla }_h\left(\overline{\overrightarrow{u}}\varphi \right)+\frac{g}{D}{\nabla }_h\overline{\rho }{\int }_{-H}^{\eta }z\left(\overrightarrow{u}-\overline{\overrightarrow{u}}\right)dz-\frac{g}{D}{\int }_{-H}^{\eta }\left(\eta -\frac{D}{2}-z\right)\left(\overrightarrow{u}-\overline{\overrightarrow{u}}\right){\nabla }_h\overline{\rho }dz\hfill \\ -\frac{g}{D}{{\displaystyle \int }}_{-H}^{\eta }\left(\eta -\frac{D}{2}-z\right)\left(w-\overline{w}\right)\frac{\partial \left(\rho -\overline{\rho }\right)}{\partial z}dz+\frac{{\rho }_0}{D}{{\displaystyle \int }}_{-H}^{\eta }{P}_{b}dz-\frac{{\rho }_0}{2}\left(P_b^s+P_b^b\right)\hfill \\ +\frac{g}{D}{{\displaystyle \int }}_{-H}^{\eta }\left(\eta -\frac{D}{2}-z\right)Qdz+\frac{g}{D}{{\displaystyle \int }}_{-H}^{\eta }\left(\eta -\frac{D}{2}-z\right){\nabla }_h\left(K_h{\nabla }_h\rho \right)dz\hfill \end{array}$$

(6)

The eight terms on the right-hand side of Equation (6) are advection, depth-mean straining, non-mean straining, vertical advection, vertical mixing, buoyancy fluxes, inner sources or sinks potential density, divergence of horizontal turbulent transport terms, respectively. Since the depth-mean straining, the advection, and the vertical mixing terms in the estuarine and coastal waters have greater influence on stratification than the other terms, the simplified 3D potential energy difference equation is as follows [44,45]:

$$\begin{array}{c}\hfill \frac{\partial \varphi }{\partial t}=\underset{\mathrm{longitudinal} \mathrm{advection}: {\varphi }_{Ax}}{\underbrace{\frac{g}{H}{\int }_{-h}^{\eta }\overline{u}\frac{\partial \left(\rho -\overline{\rho }\right)}{\partial x}zdz}}+\underset{\mathrm{lateral} \mathrm{advection}: {\varphi }_{Ay}}{\underbrace{\frac{g}{H}{\int }_{-h}^{\eta }\overline{v}\frac{\partial \left(\rho -\overline{\rho }\right)}{\partial y}zdz}}+\underset{\mathrm{longitudinal} \mathrm{depth}-\mathrm{mean} \mathrm{straining}: {\varphi }_{\mathrm{Sx}}}{\underbrace{\frac{g}{H}\frac{\partial \overline{\rho }}{\partial x}{\int }_{-h}^{\eta }\left(u-\overline{u}\right)zdz}}\hfill \\ \hfill +\underset{\mathrm{lateral} \mathrm{depth}-\mathrm{mean} \mathrm{straining}: {\varphi }_{\mathrm{Sy}}}{\underbrace{\frac{g}{H}\frac{\partial \overline{\rho }}{\partial y}{{\displaystyle \int }}_{-h}^{\eta }\left(v-\overline{v}\right)zdz}}+\underset{\mathrm{vertical} \mathrm{turbulent} \mathrm{buoyancy} \mathrm{flux}: {\varphi }_M}{\underbrace{\frac{g}{H}{{\displaystyle \int }}_{-h}^{\eta }{K}_v\frac{\partial \rho }{\partial z}dz}}\hfill \end{array}$$

(7)

3. Results

3.1. Model Validations

The model performance for the simulation of tidal elevation, velocity, and salinity is validated in this section. Observational data of tidal elevation, velocity, and salinity in the summer of 2007 are collected. The measurement locations are at the same sites for both spring and neap tides. The two in situ observations both lasted for more than 25 h with a one-hour interval. The surface layer velocity and salinity are observed at 0.5 m below the water surface, while the bottom layer velocity and salinity are measured at 0.5 m above the sea bottom. The instruments used for the tidal level, current speed, direction, and salinity observations were Tidal Gauge Recorder-2050, SLC9-2-type direct reading current meter, and a CTD (Conductance Temperature Depth) sensor, respectively.

The comparisons of tidal elevation, velocity, and salinity between the modeled results and measurements are displayed in Figure 4, Figure 5, Figure 6 and Figure 7. The model skill is evaluated quantitatively against the observed time series of tidal elevation, velocity, and salinity. The model skill score (SS) adopted by Allen et al. [46] is utilized here, which is

$$SS=1-\frac{{\displaystyle \sum _{i=1}^n{\left(M_i-O_i\right)}^2}}{{\displaystyle \sum _{i=1}^n{\left(\left|M_i-\overline{O}\right|+\left|O_i-\overline{O}\right|\right)}^2}}$$

(8)

where M and O represent modeled and observed data, respectively, and n indicates the number of samples. The skill score exhibits an indicator of model-data agreement, with a minimum value of zero demonstrating complete disagreement and a maximum value of one demonstrating perfect agreement. The performance levels of SS are categorized as follows: >0.65 “excellent”; 0.65–0.5 “very good”; 0.5–0.2 “good”; and <0.2 “poor” [47]. The quantified evaluation of the agreement is also assessed by root-mean-squared error (RMSE) and mean-normalized bias (MNB, Figure 4, Figure 5, Figure 6 and Figure 7).

The statistics of tidal elevation, velocity, and salinity are shown explicitly in Figure 4, Figure 5, Figure 6 and Figure 7. Most modeled results show acceptable agreements with observations. For tidal elevation, the model skill scores alter from 0.81 to 0.95 (Figure 4), which is regarded as “excellent” performance. The maximum values of RMSE and MNB are 0.10 m and −0.10 (Figure 4) indicating that the model is robust in simulating the tidal elevations in the PRE.

The modeled velocity follows the observational data reasonably well for both neap and spring tides, and the model skill scores are larger than 0.57 at all the locations (Figure 5), which are regarded as “excellent” or “very good” performance. Generally, the discrepancies of the surface layers are less than those of the bottom layers, especially for the flow velocity. The observed and modeled data are not taken exactly at the same depth. This impact combined with the grid resolution might result in some large discrepancies. However, the relative values of large skill scores indicate reasonable performance. Skill scores in the current direction are all larger than 0.87 (Figure 6). Consequently, the modeled results are conducted as reasonable. The maximum values of RMSE are 0.20 m/s and 10.89° (Figure 5) for flow velocity and direction respectively; the maximum values of MNB are −0.20 and 0.19 (Figure 6), respectively. These indicate that the model has reliable performance in simulating both the surface and bottom circulations in the PRE.

The modeled salinity varies in almost the same scope as the measured salinity within a flood–ebb tidal cycle. Some deviations occur during the peak flows, likely owing to the complex bathymetry of the tidal flats, deep channels, and small islands in the Pearl River Estuary (PRE). The skill scores of salinity simulation vary from 0.63 to 0.83 (Figure 7) and show a better skill near the water surface than the sea bottom. The values of RMSE range between 0.28 to 1.99 psu and the values of MNB range from −0.21 to 0.12 (Figure 7), indicating that model results are capable of capturing the variations in observed salinity.

In general, the model results reflect good agreement with the observed values with acceptable errors and allowed us to carry out numerical simulations to diagnose the hydrodynamics in the PRE. Based on the model results, variations in tidal current, salinity, and stratification at two different time scales (flood–ebb and spring–neap tidal cycles) are described and discussed in the following sections.

3.2. Along-Channel Velocity and Salinity Structure

To better examine the variations in 3D hydrodynamics, both tidal current and salinity along the selected longitudinal section (the location is shown in Figure 1) are analyzed. Longitudinal distributions of velocity and salinity are illustrated in Figure 8. It should be noted that the mean freshwater discharge is up to 2.3 × 10⁴ m³·s⁻¹ during the wet season study period and generates a major Region of Freshwater Influence (RoFI) over a large water area.

The tidal current and salinity profile exhibit remarkable alteration over a flood–ebb cycle. Station HM, located near the bay head, is mainly controlled by runoff in the wet season, with the values of salinity less than 0.1 psu (Figure 8). The flood flow is obviously smaller than the ebb flow due to the interaction between the river discharge and the longitudinal estuarine circulation. For example, on a spring tide, the maximum surface ebb flow appears at low slack with a current speed of 0.95 m·s⁻¹ (Figure 8d) while the maximum surface flood flow appears at high slack with a current speed of 0.41 m·s⁻¹ (Figure 8b). Similarly, on a neap tide, the maximum surface ebb and flood flows are 0.62 and 0.36 m·s⁻¹, respectively (Figure 8f,h). The ebb-to-flood difference of maximum surface flow on a spring tide is 0.54 m·s⁻¹, while that on a neap tide is 0.26 m·s⁻¹, demonstrating a much smaller ebb-to-flood difference on a neap tide than that on a spring tide due to the weakening of tidal action. The situation of the near-bottom water body is slightly different, the maximum flood/ebb current decreases. For example, on a spring tide, the maximum bottom ebb flow is 0.65 m·s⁻¹ (Figure 8d) while the maximum bottom flood flow is 0.28 m·s⁻¹, the ebb-to-flood difference is 0.37 m·s⁻¹ (Figure 8b). The gravitational circulation greatly enhances the surface seaward flow and finally leads to the increase in the maximum surface ebb current velocity. The topographic effect offsets the landward flow near the bottom, resulting in a slight decrease in the maximum flood current velocity at the bottom.

Station LD1 is located at the upper part of the inner Lingding Bay (Figure 8). The maximum salinity occurs in the bottom layer of the spring tide, which is 1.41 psu, accompanied by the maximum vertical salinity difference of 1.34 psu. Station LD1 is also dominated by freshwater in the wet season, except for the maximum flood of spring tide. On a spring tide, the maximum surface ebb flow appears at maximum ebb with a current speed of 0.78 m·s⁻¹ (Figure 8c) while the maximum surface flood flow appears at maximum flood with a current speed of 0.40 m·s⁻¹ (Figure 8a). The maximum bottom ebb and flood flows are 0.46 and 0.27 m·s⁻¹, respectively (Figure 8b,d). On a neap tide, the maximum surface ebb and flood flows are 0.36 and 0.24 m·s⁻¹, respectively (Figure 8f,h), and the maximum bottom ebb and flood flow are 0.23 m·s⁻¹ and 0.16 m·s⁻¹, respectively. Similar to station HM, the longitudinal gravitational circulation enhances the surface seaward flow and increases the surface ebb current velocity. At the same time, due to the influence of shoreline and topography, the influence of gravitational circulation near the bottom is offset, and the bottom current velocity is smaller. In addition, because of the larger channel width, the bottom and surface velocities are smaller than those at station HM.

As shown in Figure 8, Station LD2 is located at the upper part of the outer Lingding Bay with the maximum vertical salinity difference of 24.1 psu at the low slack on a neap tide. On a spring tide, the maximum surface ebb and flood flows are 1.04 m·s⁻¹ and 0.45 m·s⁻¹, respectively, while the maximum bottom ebb and flood flow are 0.35 m·s⁻¹ and 0.38 m·s⁻¹, respectively. The bottom current velocity of the ebb tide is smaller than that of the flood tide, which reflects the effect of gravitational circulation. Similarly, on a neap tide, the maximum surface ebb and flood flows are 0.46 m·s⁻¹ and 0.18 m·s⁻¹, respectively, while the maximum bottom ebb and flood flow are 0.12 m·s⁻¹ and 0.29 m·s⁻¹, respectively.

It can be seen from Figure 8 that Station DH is located in the middle part of the outer Lingding Bay with a maximum salinity of 30.86 psu at the high slack on a spring tide. However, the vertical salinity difference is only 13.01 psu at this point, while the maximum vertical salinity difference is 25.16 psu at maximum ebb on a neap tide. In general, the stratification is stronger on a neap tide than on a spring tide. On a spring tide, the maximum surface flood flow occurs at maximum ebb (1.35 m·s⁻¹), and the maximum surface flood flow at maximum flood (0.49 m·s⁻¹). Meanwhile, the maximum bottom ebb and flood flows are 0.5 m·s⁻¹ and 0.67 m·s⁻¹, respectively. On a neap tide, the ebb flow dominates in the surface layer, with a maximum value of 0.81 m·s⁻¹, and the maximum bottom flood and ebb current velocities are 0.41 m·s⁻¹ and 0.27 m·s⁻¹, respectively. It can be seen that the runoff effect is enhanced due to the shoreline contraction at station DH.

According to Figure 8, Station WH is located near the entrance of the estuary and the salinity of the surface and bottom layers are high and barely change with time. The average salinity of the surface and bottom layers is 18 psu and 31.2 psu, respectively. At this station, the ebb current is dominant during the spring and neap tides, and the maximum flow velocity is 0.74 m·s⁻¹, which occurs at low slack on a spring tide. In the bottom layer, on a spring tide, the maximum ebb flow is 0.22 m·s⁻¹ at maximum ebb, and the maximum flood flow is 0.22 m·s⁻¹ at low slack. On a neap tide, the maximum bottom flood and ebb flows are 0.14 m·s⁻¹ and 0.08 m·s⁻¹, respectively, which occur at maximum flood and maximum ebb, respectively. It can be seen that due to the topographic effect at station WH, runoff and gravitational circulation have a great influence in the wet season.

In general, the maximum surface-to-bottom salinity difference is 25.21 psu, displaying strong stratification in the wet season. The spring–neap tidal alteration in stratification is relatively small due to huge river discharge. The variation in stratification in an intratidal tidal cycle is also clearly displayed in Figure 8. For example, on a neap tide, the surface-to-bottom salinity difference at station LD2 changes from 20.64 to 24.02 psu. The vertical salinity difference increases during the ebb flow while decreases during the flood flow on a neap tide. However, stratification is relatively weak at low and high slack waters, but relatively strong at the maximum ebb and maximum flood tides on a spring tide. Moreover, the inclination of the salinity contours also demonstrates the intensity of stratification. As shown in Figure 8, during the flood on a spring tide, the horizontal length of the 10 psu contour is about 26.5 km at maximum flood and reduces to about 23.8 km at high slack, suggesting weakening stratification. During the ebb, the horizontal length of the 10 psu contour amounts to 33.5 km at maximum ebb, and then increases to 35.5 km at low slack, and stratification reaches its strongest at the same period.

3.3. Stratification Indicated by Richardson Number

Stratification plays an essential part in modulating vertical mixing, estuarine circulation, salt balance, and sediment transportation, which has been widely established by previous works [3,39,48]. To better evaluate the spatial and temporal alterations of active mixing, shear and stratification, the square of buoyancy frequency (N²), the square of shear (S²), and the gradient Richardson number (Ri) are calculated and illustrated in their logarithms to analyze the intratidal, spring–neap variations in salinity stratification in the West Channel.

The variations in Ri, N², and S² at HM station in the wet season are shown in Figure 9, Figure 10 and Figure 11. The log₁₀(Ri/0.25) at HM station is low throughout the study period with a value less than −4 due to the large runoff in the wet season. Similarly, the buoyancy frequency square is small and remains below 1 × 10⁻⁶ s⁻². In addition, the velocity shear term is also small with a maximum value of 1 × 10^−2.2 s⁻² at low slack on a spring tide. Moreover, the velocity shear tends to be smaller at the surface and larger at the bottom. For example, at low slack on a spring tide, the velocity shear changes from 1 × 10^−4.1 s⁻² at the surface to 1 × 10^−2.2 s⁻² at the bottom, which is closely related to the local topography and channel width. Overall, small N² makes small Ri, resulting in strong turbulent mixing.

It can be seen from Figure 9, Figure 10 and Figure 11 that the value of Ri at station LD1 has an obvious flood–ebb cycle. On a spring tide, Ri gradually increases from maximum flood to high slack, and the value of log₁₀(Ri/0.25) changes from below 0 to above 0, and reaches the maximum value of 1.67 at the surface layer. From high slack to maximum ebb, there is little change in log₁₀(Ri/0.25), basically maintaining above zero. From maximum ebb to low slack, it gradually decreases and finally reduces to a value between −2 and −3.4. It can be seen that the turbulent mixing is stronger at maximum flood and low slack, while weaker at high slack and maximum ebb. Vertically, the value of Ri is larger in the surface layer than that in the bottom layer, which is basically consistent with the shear term. This is mainly because the friction leads to greater shear in the bottom boundary. Similar to HM station, the value of N² is small and maintains below 1 × 10⁻⁵ s⁻² due to the impact of the huge amount of freshwater. The shear term has obvious tidal cycle changes, and the S² at low and high slack waters are 1~2 orders of magnitude higher than those at maximum ebb and maximum flood tides. At the same time, the shear term of the surface layer is 2~3 orders of magnitude higher than that of the bottom layer.

At LD2 station, Ri is basically larger than the critical value 0.25, and the value of Ri here does not change much during the spring–neap tidal cycle, and the water stability is relatively high. In addition, the values of N² in the bottom layer at maximum flood and high slack are small, which is about one order of magnitude smaller than those at maximum ebb and low slack. For example, on a spring tide, the values of N² are 1 × 10^−2.9 s⁻² and 1 × 10⁻⁴ s⁻² at maximum ebb and maximum flood, respectively. Moreover, the values of N² in the middle of the surface layer are relatively large with a maximum value of 1 × 10^−1.2 s⁻², indicating obvious stratification. The shear term at this station reaches the largest in the near-surface layer, and 1~2 orders of magnitude smaller in the surface and bottom layers.

At station DH, Ri is higher in the surface and middle layers while lower in the bottom layer. The values of log₁₀(Ri/0.25) vary from 2.1 to 2.47 in the surface layer, showing high stability and weak turbulent mixing. The values of log₁₀(Ri/0.25) in the bottom water vary from −0.92 to 0.63 on a spring tide, while changing from −0.32 to −1.19 on a neap tide, demonstrating low water stability. In addition, the values of N² show a similar pattern to log₁₀(Ri/0.25) in the vertical, with the values mainly below 1 × 10⁻² s⁻². The stratification is strong and remains stable in the surface and middle layers. The shear terms at the DH station do not change much at maximum ebb, maximum flood, and low slack. However, the values of S² are 2~3 orders of magnitude lower in the bottom layer at high slack than those at the other three moments, indicating that the vertical flow velocity of the bottom layer is stable at this time.

In the vertical, the value of Ri at the WH station is the largest in the surface layer, with log₁₀(Ri/0.25) varying from 2.38 to 2.80 demonstrating weak but stable turbulent mixing. The value of Ri of the middle layer is slightly lower than that of the surface layer while in the bottom layer is less than zero due to tidal stirring. In terms of N², the values in the surface and middle layers are about 1 × 10⁻² s⁻² showing strong stratification, while the bottom layer values are small with the minimum value of 1 × 10^−6.3 s⁻² at high slack on a neap tide indicating weak stratification due to the influence of topography and geographical location. In addition, the shear term is larger in the near-surface layer, ranging from 1 × 10⁻¹ s⁻² to 1 × 10^−2.5 s⁻², which is about two orders of magnitude larger than the other shear terms in the vertical.

In general, the stratification in the wet season is strong due to the huge amount of freshwater from the Pearl River. Under the control of freshwater discharge, N² is small in the upper part of the inner Lingding Bay and the weak stratification leads to a low value of Ri and high turbulent mixing. The relatively high-value area (Ri > 0.25) almost occupied the middle and upper layers of the remainder of the inner Lingding Bay and the outer Lingding Bay. The turbulent mixing of the middle and upper layers is low, except for the bottom layer where Ri is lower than the critical value due to tidal stirring. Shear instability occurs at the salinity interface and the bottom interface in the upper part of the inner Lingding Bay, resulting in water mixing. In addition, shear instability also occurs in the rest of the water layer at each time period, but excessive stratification suppresses shear instability, resulting in persistently weak turbulent mixing. During the wet season, the water stability is relatively high during the spring–neap tidal cycle. It is worth noting that the turbulent mixing near the bottom layer of LD2 is weak during the spring tide, which is mainly caused by the small shear term there.

4. Discussion

4.1. The Change Rate of Potential Energy Anomaly

To demonstrate the contributions of different physical mechanisms to the salinity stratification, the five contribution terms in Equation (7) are calculated and compared cautiously using modeled tidal level, velocity, and salinity after careful calibration and validation. Spring–neap and intratidal variabilities of the terms along the West Channel are depicted in Figure 12.

The time series of the differential change rate of 3D potential energy at station LD1 is shown in Figure 12a. Under the control of the river effect, the longitudinal and lateral advection as well as the longitudinal and lateral depth-mean straining terms are close to zero most of the time, with occasional fluctuations. The minimum and maximum values of the longitudinal advection term are about −4.3 × 10⁻³ and 2.2 × 10⁻³ W·m⁻³ on a spring tide, while they are about −1.15 × 10⁻³ and 0.74 × 10⁻³ W·m⁻³ on a neap tide. The mean values on the spring/neap tides are −0.015 × 10⁻³ W·m⁻³ and 0.031 × 10⁻³ W·m⁻³. For the longitudinal and lateral depth-mean straining terms, the average values in the spring and neap tides are 0.035 × 10⁻³ W·m⁻³ and −0.002 × 10⁻³ W·m⁻³, 0.046 × 10⁻³ W·m⁻³, and 0.045 × 10⁻³ W·m⁻³, respectively. Unlike the other terms, the values of the vertical mixing term are always negative. The maximum values of the term during the spring and neap tides are −6.9 × 10⁻³ W·m⁻³ and −1 × 10⁻³ W·m⁻³, respectively. On the whole, the vertical mixing term makes the largest contribution at LD1.

The time series of the differential change rate of 3D potential energy at station LD2 is displayed in Figure 12b. During the wet season, the longitudinal advection is greater than zero most of the time, and the average values of the longitudinal advection in the spring and neap tides are 6.6 × 10⁻³ W·m⁻³ and 5.6 × 10⁻³ W·m⁻³, respectively. The differential change rate of potential energy caused by lateral advection varies from −6.9 × 10⁻³ to 11.9 × 10⁻³ W·m⁻³ and −8.9 × 10⁻³ to 4.1 × 10⁻³ W·m⁻³ during the spring and neap tides, respectively. Similar to the longitudinal advection, the longitudinal depth-mean straining increases the water stratification at most of the periods, and the average values of the longitudinal depth-mean straining are 3.2 × 10⁻³ W·m⁻³ and 2.9 × 10⁻³ W·m⁻³ on the spring/neap tides, respectively. However, the overall fluctuation of vertical mixing and lateral straining is small, and the effect can be neglected. Under the influence of longitudinal advection, lateral advection, and longitudinal average straining, the differential change rate of the overall potential energy fluctuates dramatically during the tidal cycle, ranging from −7 × 10⁻³ to 24.9 × 10⁻³ W·m⁻³ and −8.1 × 10⁻³ to 18.7 × 10⁻³ W·m⁻³ in the spring and neap tidal periods, respectively. The mean values are 4.5 × 10⁻³ W·m⁻³ and 3.4 × 10⁻³ W·m⁻³, which promoted stratification as a whole.

The differential change rate of 3D potential energy at station DH is exhibited in Figure 12c. During the wet season, the longitudinal advection mainly promotes stratification, and this effect is more obvious on a spring tide, with the average values of −3.2 × 10⁻³ W·m⁻³ and −0.6 × 10⁻³ W·m⁻³ on the spring and neap tides, respectively. The lateral advection term fluctuates greatly at the DH point and plays a stratified role at most periods, varying from −8.3 × 10⁻³ to 30.15 × 10⁻³ W·m⁻³ and −9.5 × 10⁻³ to 23.16 × 10⁻³ W·m⁻³ during the spring and neap tides, respectively. The longitudinal and lateral depth-mean straining terms are mainly related to forcing stratification, but their values are smaller. The average values are 1.2 × 10⁻³ W·m⁻³ and 1.6 × 10⁻³ W·m⁻³ in the spring and neap tides, and the calculated values are 1.25 × 10⁻³ W·m⁻³ and 1.99 × 10⁻³ W·m⁻³. The mixing effect of the vertical mixing mechanism is relatively small, with average values of −1.1 × 10⁻³ W·m⁻³ and −1.2 × 10⁻³ W·m⁻³ during the spring and neap tides. On the whole, the differential change rate of the total potential energy during the wet season is mainly affected by the longitudinal and lateral advection terms, with the average values of 1.8 × 10⁻³ W·m⁻³ and 5.3 × 10⁻³ W·m⁻³ during the spring and neap tides.

The variation rate of potential energy difference at station WH can be seen in Figure 12d. Since this station is close to the open sea, the density gradient is small, and all values are small in the wet season. The average change rate of the total potential energy difference is only 0.3 × 10⁻³ W·m⁻³ and 0.03 × 10⁻³ W·m⁻³ on spring and neap tides.

4.2. Physical Mechanisms Contributing to Spring–Neap and Intratidal Variabilities of Salinity Stratification

The values of the total potential energy difference are mainly positive promoting stratification within the middle and lower parts of the West Channel. In general, the longitudinal depth-mean straining makes the most significant contribution to stratification within the West Channel with the maximum tidally averaged values of 2.09 × 10⁻² W·m⁻³ and 2.59 × 10⁻² W·m⁻³ on spring and neap tides at station LD2. As mentioned by Burchard and Hofmeister [14], the longitudinal depth-mean straining superposes the combined effect of gravitational circulation and tidal straining. Interacting with the shear flow induced by river discharge, longitudinal density gradient can enhance stratification persistently [18,49,50]. Thus, the large positive value of the longitudinal depth-mean straining mainly results from the combined effect of river and circulation.

It could be unveiled that longitudinal depth-mean straining plays a significant part in spring–neap tidal variability of stratification within the West Channel, because the tidally averaged longitudinal depth-mean straining on a neap tide is about 1.5–5.0 × 10⁻³ W·m⁻³ larger than that on a spring tide. This indicates that longitudinal depth-mean straining generates weaker stratification on a spring tide but stronger stratification on a neap tide. Owing to the very fact that the spring–neap tidal variability of freshwater runoff is not evident, spring–neap tidal variability of the tidally averaged longitudinal depth-mean straining mainly arises from spring–neap tidal variability of circulation.

It should be noted that longitudinal and lateral advections impose their influence as well, but the relative importance change remarkably in space. Generally, the longitudinal advection makes the second contribution with the maximum tidally averaged values of 1.02 × 10⁻² W·m⁻³ and 4.09 × 10⁻³ W·m⁻³ on spring and neap tides. However, the effects of the lateral advection, lateral depth-mean straining, and vertical mixing cannot be ignored owing to their obvious spring–neap variability. For example, the maximum tidally averaged values of lateral depth-mean straining and lateral advection at station LD2 are 2.21 × 10⁻³ W·m⁻³ and 2.96 × 10⁻³ W·m⁻³, respectively, and the lateral advection prevails on a neap tide at station DH. This could be universal in channelized estuaries since the longitudinal density gradient is generally much smaller than the lateral density gradient owing to differential advection. The vertical mixing terms are negative throughout time with the maximum tidally averaged value of −1.31 × 10⁻³ W·m⁻³ on a spring tide promoting mixing.

The values of the total potential energy difference demonstrate obvious intratidal variability. For instance, at station LD2, the average values of total potential energy difference at high and low slacks are 5.0 × 10⁻² and 5.0 × 10⁻³ W·m⁻³, respectively. In general, the longitudinal advection and longitudinal depth-mean straining are the most important contributors to intratidal variability of salinity stratification. For example, at station LD2, the longitudinal advection term makes the most contribution at high slack and maximum flood while the longitudinal depth-mean straining term makes the most contribution at maximum ebb and low slack. The vertical mixing demonstrates evident intratidal variability with larger negative values dominating within the West Channel at maximum ebb, whereas smaller values prevail at maximum flood. Such intratidal variability of the vertical mixing largely generates from the intratidal variability of tidal stirring owing to flood–ebb tidal current asymmetry.

5. Conclusions

An analysis of stratification alterations in the Pearl River Estuary (PRE), China, is conducted in this study. Based on the modeled hydrological data during neap and spring tides in the wet season, the gradient Richardson number and the potential energy anomaly are calculated to explore the flood–ebb and spring–neap variations in salinity stratification.

The temporal and spatial variations of the salt wedge are remarkable along the West Channel. The horizontal length of the 10 psu contour is much shorter on a neap tide than that on a spring tide. According to the analysis of the gradient Richardson number, shear instability typically appears at the fresh-water–salt-water interface along the middle and lower parts of the channel.

Based on the simplified 3D differential equation of potential energy, the physical mechanism of stratification in the PRE is further studied. The results show that the longitudinal and lateral advection, the longitudinal and lateral depth-mean straining, and the vertical mixing jointly determine the stratification trend of the water body. Generally speaking, the effect of longitudinal advection and longitudinal depth-mean straining is greater than that of lateral terms, as shown at LD2 of Figure 12. However, the lateral terms at some times are similar to the longitudinal terms as displayed at LD1 of Figure 12. The lateral action in this study is mainly induced by the lateral variations in bathymetry in Lingding Bay. As shown in Figure 1, there are four outlets, two deep channels, and three shoals in Lingding Bay. The lateral exchanges among the outlets and between the shoal–channel system are significant. Moreover, the lateral mechanisms might be more complicated if wind forcings are included in the model. In general, the lateral action is an indispensable link for the study of stratification in the PRE.

With much less river discharge from the Pearl River, the salt wedge might reach further landward, and the salinity stratification might be much weaker in the dry season than that in the wet season. Other aspects including the effects of waves, winds, human interventions, and sediment-induced stratification, were not further investigated and await future explorations.