Effects of Turbulence Parameterization on Hydrodynamics and Sediment Transport in Tidal Channels: A Case Study of Yamen Channel in the Northern South China Sea
https://doi.org/10.1007/s11804023003279

Abstract
In this study, we conducted numerical experiments to examine the effects of turbulence parameterization on temporal and spatial variations of suspended sediment dynamics. Then, we applied the numerical model to the Yamen Channel, one of the main eight outfalls in the Pearl River Delta. For the field application, we implemented the kε scheme with a reasonable stability function using the continuous deposition formula during the erosion process near the watersediment interface. We further validated and analyzed the temporalspatial suspended sediment concentrations (SSCs). The experimental results show that under specified initial and boundary conditions, turbulence parameterization with stability functions can lead to different vertical profiles of the velocity and SSC. The kε predicts stronger mixing with a maximum value of approximately twice the kkl. The kkl results in smaller SSCs near the surface layer and a larger vertical gradient than the kε. In the Yamen Channel, though the turbulent dissipation, turbulent viscosity and turbulence kinetic energy exhibit similar trends, SSCs differ significantly between those at low water and high water due to the tidal asymmetry and settling lag mechanisms. The results can provide significant insights into environmental protection and estuarine management in the Pearl River Delta.Article Highlights• Numerical experiments were conducted to examine the effects of turbulence parameterization on suspended sediment dynamics.• The k−ε scheme predicts stronger mixing with a maximum value of approximately twice the k−ε. The k−ε results in smaller SSCs near the surface layer and a larger vertical gradient than the k−ε.• Suspended sediment concentrations (SSCs) differ significantly between those at low water and high water due to the tidal asymmetry and settling lag mechanisms in the Yamen Channel. 
1 Introduction
Turbulence parameterization is highly uncertain in both physical and sediment numerical models (Geyer and Ralston, 2015). In the past decades, the major focus on turbulence modelling is to suggest equations closures at various sophistication levels. Numerous modifications to these models have been proposed (Rodi, 1984; Chen and Beardsley, 1998; Umlauf et al., 2003; Burchard et al., 2014; Liu et al., 2022).
The k−kl turbulence closure scheme (Mellor and Yamada level2.5, MY2.5, MY 1982) has been used widely in estuaries and coastal oceans. In the MY2.5, the vertical turbulent viscosity is calculated by solving the turbulence kinetic energy (TKE) equation and turbulence length scale equation. This scheme seems to work well for tidalinduced mixing in shallow waters (Chen and Beardsley, 1998), but further validation is needed for the estuarine regions. Several studies have been carried out to improve this turbulent model (Galperin et al., 1988; Burchard et al., 1998; Burchard and Bolding, 2001). A modification of the macroscale equation has been introduced, together with a new wall proximity function to investigate the vertical structure of open channel flows (Blumberg et al., 1992). Rodi (1984) significantly advanced the k−ε. Recent contributions paid particular attention to the buoyancy parameterization and model comparison (Burchard et al., 1998; Burchard and Bolding, 2001; Burchard et al., 2014). Umlauf et al. (2003) suggested a generalization of a class of differential lengthscale equations typically used in secondorder turbulence models for oceanic flows. Generally used models, such as the k−ε and the MY2.5, can be treated as special cases of this generic model. Recently, some studies have investigated physical and numerical turbulentmixing behavior through numerical experiments (Reffray et al., 2015; Ralston et al., 2017; Costa et al., 2017; Tu et al., 2019; Choi et al., 2021).
Since vertical diffusion is one of the key processes for suspended sediment concentration (SSC) in estuaries, the uncertainty due to various turbulence parameterization schemes should be considered in the assessment of the reliability of the sediment transport model. Warner et al. (2005) evaluated the performance of four turbulence closure models and their influences on suspended sediment transport using a generic length scale method with several ideal cases based on ROMS. In this contribution, we carry out numerical experiments to investigate the impacts of turbulence parameterization on the suspendedsediment profile, together with different stability functions. We conducted the experiments by running FVCOM with various turbulence closure modules implemented in the General Ocean Turbulence Model (GOTM).
In recent years, since the reclamation on both sides of the waterway and the evolution of shoals, the deposition in local thalweg areas has aggravated to obstruct navigation (Lima et al., 2021; Wu et al., 2022). Yang and Wang (1994), and Chen et al. (2003) have conducted numerical simulations of the tidal wave and sediment transport in the Pearl River, but without systematic analyses in the Yamen Channel, which is one of the eight outfalls of the Pearl River. It is of great significance to investigate the hydrodynamics and suspended sediment transport in this area. Recently, some studies have been carried out in the Pearl River estuary (Zhang et al., 2019; Yang et al., 2022). To further investigate the turbulence simulation and better understand its impact on the sediment dynamics in the Pearl River Estuary, we conduct numerical simulations to investigate the sensitivity of turbulence parameterization on vertical distributions of velocity and SSC in the Yamen Channel. This paper is structured as follows: after describing the study area and observations (Section 2), we present three main modules of the numerical model together with model configuration (Section 3). Section 4 gives the sensitivity analysis of turbulence parameterization and model validation, followed by a detailed investigation of the spatialtemporal variations of SSC and intratidal variations of turbulent parameters, vertical velocity, and SSC. The conclusions are listed in Section 5.
2 Study area and observations
Yamen Channel, located at the head of the Huangmao Sea (HM Sea), is the westernmost one of the eight outfalls (YM: Yamen, HTM: Hutiaomen, JTM: Jitimen, MDM: Modaomen, HM: Hengmen, HQM: Hongqimen, JM: Jiaomen, HuM: Humen) of the Pearl River Delta (Figure 1). In Figure 1, Z2, Z3 and Z4 are sites for water levels. A3, A4, and A5 are for velocity and SSCs. The runoff of the Tanjiang River (TR) flows into HM Sea through the Yamen Channel. This typical tidaldominated channel is an important navigable waterway, which is surrounded by extensive intertidal salt marshes and numerous tidal creeks. The annual discharge is about 196 × 10^{8} m^{3}, which covers 6% of the total runoff of the Pearl River (Figure 1). Compared to the strong tidal current, the river flow appears weak. The average tidal range is approximately 1.24 m. The model domain is up to the Tanjiang River, and reaches the HM Sea downstream. A comparison between the Yamen and the Hutiaomen (HTM) at the HM Sea is shown in Table 1.
Items Branches HTM YM Total amount of 8 branches Annual discharge (10^{8} m^{3}) 202.0 196.0 3 260.0 Percent (%) (the discharge of the Pearl River) 6.2 6.0 100.0 Annual sediment transport rate (10^{4} t) 509.0 363.0 8 872.0 Percent (%) (the sediment of the Pearl River) 7.2 5.1 100.0 Annualaveraged tidal range (m) 1.20 1.24 The maximum tidal range 2.66 2.95 Annual flood tide volume (10^{8} m^{3}) 57.0 636.0 3 763.0 The ratio of total runoff and flood tide volume 3.54 0.31 The medium diameter of the suspended load is less than 0.1 mm and the diameter of sediment particles from the upstream during ebb tide is larger than that from the outer sea during flood tide. In flood season, flow with high SSC enters the Yamen channel. The mean grain size is about 5 μm. The sediment is mainly transported in suspension. The particle diameter of bed sediment between 0.0050.05 mm covers about 50%, and the median particle size is less than 0.05 mm.
3 Numerical Model
3.1 Hydrodynamic model
We adopted an unstructuredgrid, threedimensional primitive equation community ocean model FVCOM (Finite Volume Community Ocean Model) in this study (Chen et al., 2007). It solves the governing equations in unstructured triangular mesh by finitevolume method, which provides accurate conservations of mass, momentum, heat, and salt. This model has been applied to a number of estuaries and coastal oceans characterized by highly irregular geometry, large intertidal salt marshes, and steeply sloping bottom topography (Ding et al., 2013; Lai et al., 2015; Sun et al., 2016).
3.2 Suspended sediment model
Suspended sediment is simulated by the threedimensional suspended sediment equation in σcoordinate:
$$ \begin{array}{r} \frac{\partial C}{\partial t}+\frac{\partial u C}{\partial x}+\frac{\partial v C}{\partial y}+\frac{1}{H} \frac{\partial w C}{\partial \sigma}=\frac{\partial}{\partial x}\left(A_h \frac{\partial C}{\partial x}\right)+ \\ \quad \frac{\partial}{\partial y}\left(A_h \frac{\partial C}{\partial y}\right)+\frac{1}{H} \frac{\partial}{\partial \sigma}\left(\frac{K_h}{H} \frac{\partial C}{\partial \sigma}\right)+\frac{1}{H} \frac{\partial C \omega_s}{\partial \sigma} \end{array} $$ (1) where A_{h} is the horizontal diffusion coefficient, K_{h} is the vertical diffusion coefficient, C is the suspended sediment concentration, ω_{s} is the settling velocity, H is total water depth, and (u, v, w) represent horizontal and vertical velocity, respectively.
Sediment particles near the bed can be resuspended when the bottom bed stress exceeds critical stress. The bottom boundary condition is determined by nearbed sediment flux (van Rijn, 2005; Li et al., 2008):
$$ F_s=ED $$ (2) where F_{s} is nearbed sediment flux, E is erosion flux, and D is deposition flux. The deposition rate (D in kg/(m^{2}s)) is determined as the mass that is removed from the suspended load and integrated into the bed. It depends on the particle's sinking velocity ω_{s}, and the near bottom concentration C_{b} in kg/m^{3}:
$$ D=\omega_s C_b $$ (3) Erosion is an instantaneous process that is initiated when the bed shear stress exceeds a critical value τ_{ce}, which is derived from the critical Shields parameter (Warner et al., 2008).
$$ E=\left\{\begin{array}{cl} Q\left(1P_b\right) F_b\left(\frac{\tau_b}{\tau_{c e}}1\right) & \text { if } \tau_b>\tau_{c e} \\ 0 & \text { if } \tau_b \leqslant \tau_{c e} \end{array}\right. $$ (4) where Q is the erosion flux, P_{b} is the bottom porosity, F_{b} is the fraction of sediment and τ_{ce} represents critical shear stress.
In terms of the initial conditions, the suspended sediment concentration is zero for the normal component on the fixed boundary. At the outflow open boundary, it is determined as $ \frac{\partial C}{\partial t}+v_n \frac{\partial C}{\partial n}=0$, vn is the normal velocity component, and n is the normal direction to the open boundary. For the inflow, C = C_{0}, where C_{0} is the given values.
The sediment bed is represented by 3D arrays with a userspecified, constant number of layers beneath each horizontal model cell. At the beginning of each time step, activelayer thickness z_{a} is calculated based on the relation (Warner et al., 2005, W2005)
$$ z_a=\max \left[k_1\left(\tau_b\tau_{c e}\right) \rho_0, 0\right]+k_2 D_{50} $$ (5) where τ_{b} is bottom stress, τ_{ce} is the critical stress for erosion, D_{50} is the median particle size, and k_{1} and k_{2} are empirical coefficients (0.007 and 6.0, respectively). The thickness of the top bed layer has a minimum thickness equivalent to z_{a}.
3.3 Turbulent closure
The Reynolds stress tensor and the turbulent tracer flux need to be parameterized to close the hydrostatic primitive equations that resulted from the Reynoldsaveraging NavierStokes equations. The transport equations are parameterized by approximations of the thirdorder moments and pressure correlation terms. Additional scaling and boundary layer assumptions further reduce the set of equations of the Reynolds stresses and turbulent scalar fluxes in the form as follows:
$$ \overline{u^{\prime} w^{\prime}}=K_M \frac{\partial U}{\partial z}, \overline{v^{\prime} w^{\prime}}=K_M \frac{\partial V}{\partial z}, \overline{w^{\prime} \rho^{\prime}}=K_H \frac{\partial \rho}{\partial z} $$ (6) with
$$ K_M=c \sqrt{2 k} l S_M+v, K_H=c \sqrt{2 k} l S_H+v_\theta $$ (7) where ρ is total density, and U and V are velocity of the mean components in the horizontal directions, respectively. u', v', and w' are the turbulent components of velocity in the horizontal and vertical directions, respectively. K_{M} is the turbulent viscosity, K_{H} is the turbulent diffusivity, and S_{M} and S_{H} are stability functions that describe the effects of shear and stratification, respectively. υ and υ_{θ} represent molecular viscosity and diffusivity, respectively. l is turbulent length scale.
The GOTM developed by Umlauf and Burchard (2003)is a community turbulence module and is continuously upgraded (Burchard et al., 2008; Reissmann et al., 2009; Umlauf and Burchard, 2011). This module includes two typical turbulence closure models, (1) k−kl equation and (2) k−ε equation. Both turbulence module groups include the original code with a Richardson number cut off between 0.2 and 1.0. The GOTM has been coupled to FVCOM through an interface library, thus one can select different turbulence schemes to replace the default setup of the MY2.5 turbulence closure scheme.
The turbulent closure schemes rely on the stability functions parameterizing the pressurestrain correlation terms of the Reynolds stress (Burchard and Bolding, 2001). Each stability function considers a different approximate equilibrium of the five pressurestrain effects, including isotropy, shear production, buoyancy production, nonisotropic contribution and vorticity contribution. Five different stability functions have been explored: the modified MY2.5 which only contains isotropy and shear productions (Galperin et al., 1988), KC (Kantha and Clayson, 1994) which includes both isotropic and nonisotropic contributions; BB which contains the isotropy shear and buoyancy turbulence productions (Burchard and Baumert, 1995); CA (Canuto et al, 2001) which includes all five pressurestrain terms; and CB (Burchard and Baumert, 1995; Canuto et al., 2001) which was deduced for an equilibrium state that turbulence dissipation balances turbulence productions. We implemented the original MY2.5, modified MY2.5 (G1988), and the KC stability functions in the k−kl equation and the other three stability functions (BB, CA, and CB) for the k−ε equation model (Table 2).
Turbulent scheme Stability functions Parameterizations Isotropy Shear production Buoyancy Production Nonisotropic Vorticity k − ε BB √ √ CA √ √ √ √ √ CB Deduced by the quadrivium of turbulent dissipation and production. k − kl MY2.5 Functions of the gradient Richardson number. G1988 √ √ KC √ √ 3.4 Model configuration
The model area is discretized by nonoverlapping triangular cells with a horizontal resolution of 2080 m. Twenty even sigma layers are separated in the vertical direction, which allows for a slightly smooth representation of finiteamplitude irregular bottom bathymetry. The total numbers of triangular elements and nodes are 6 295 and 3 563, respectively. The time step for the external mode is 0.1 s.
The model is forced by river flow at the upstream riverine boundary and by water level at the southern open boundary. Tanjiang River, Sanjiang Mouth, HuKeng, and HTM channel are all input sources of riverine discharges and SSCs. Hushan is the downstream control boundary with water level extracted from the tidal dataset TPXO7.2 (Egbert and Erofeeva, 2002). Topography is based on the nautical maps published Dec. 2004 as the model terrain(Figure 1).
The wind is weak during the simulation period, thus windinduced wave effects and wind forcing are not included in these simulations. We chose the Smagorinsky turbulence parameterization method for the horizontal diffusion coefficient.
Since the SSCs are relatively small (0.020.4 kg/m^{3}) along the channel, and the influence of turbid water on particle settling is neglected, the settling velocity is set as a constant of 0.75 mm/s for simplicity. The critical erosion shear stress is set to 0.04~0.06 N/m^{2}, and a mean value of 0.05 N/m^{2} is used in the simulations. The erosion constant is 5×10^{−}^{5} kg/(m^{2} s) (Table 3).
Horizontal resolution (m) Nodes Cells Time step (s) Time step (s) τ_{ce} (N/m^{2}) Erosion constant (kg/m^{2}s) 2080 3 563 6 295 0.1 0.75 0.05 5×10^{5} 4 Results and analysis
4.1 The sensitivity analysis of turbulence closures
We conducted sensitivity analysis based on turbulent closure methods using the k − ε equations with stability functions of CA, CB, and BB, k − kl with KC, MY2.5, and G1988, and Analysis (ANA). The ANA is an algebraic expression for the turbulent viscosity asK_{M} = ku_{*} z (1 − z/H), where k is a coefficient, u_{*}, frictional velocity, z, distance above the sea bed, and H, total water depth. Z4 is located near the Yamen Bridge with a depth of ~10 m. We chose this site to investigate the influence of turbulence closure schemes on the vertical distributions of velocity and SSC.
4.1.1 Vertical profile of turbulent variables
As a measure of turbulence intensity, TKE is one of the most important variables related to the momentum transport through the bottom boundary layer (BBL). Figure 2a shows that the k − kl with the three stability functions of KC, MY2.5, and G1988 generates almost identical patterns. The k − ε produces a notable difference among the three different stability functions, especially in the lower half of water depth. The values of CB and BB are ~2 times larger than that of CA. ANA would become more turbulent than the other methods.
Turbulent dissipations computed by different stability functions keep in good consistency, except that ANA would dissipate strongly near the BBL (Figure 2b). Generally, ANA would overestimate both TKE and turbulent dissipation especially near the BBL.
The k − kl generally underestimates turbulent viscosity by ~40% compared to that of k − ε and ANA (Figure 2c). The maximum values of k − ε and ANA are almost identical. The profile of turbulent viscosity demonstrates asymmetrical characteristics for all turbulent schemes. However, ANA and k − kl show an asymmetry biased towards the lower layer of water depth, whereas k − ε simulates an asymmetry towards the upper layer.
4.1.2 Vertical profiles of velocity and SSC
The vertical velocity profile shows two typical different patterns, especially in the surface layer (Figure 3a). The k − kl with the three different stability functions and ANA produce anomalously high surface currents. This behavior has been documented (Blumberg et al., 1992; Burchard et al., 1998; Burchard and Bolding, 2001; Warner et al., 2005). Simulations by Warner et al. (2005) in open channel show that the original MY2.5 and k − kl with parabolic wall function have abnormally high surface currents. This is consistent with the results of other turbulence schemes, including k − ε, k −ω, ANA, and k − kl with improved open channel wall function (Warner et al., 2005). In this paper, k − ε schemes (CB, CA, BB) produce results that agree closely, while the surface flow velocities obtained using the k − kl schemes including G1988, original MY2.5, KC, and ANA are excessively large. Warner et al. (2005) pointed out that the wrong wall function produces an incorrect mixing length scale, which leads to the underestimation of turbulent visibility far from the boundary. Moreover, other important consequences, besides the resultant higher depthaveraged velocities in the open channel, would also exist in the estuary. For instance, the mixing of fresh water and saline water in the estuary has a critical impact on the estuarine gravitational circulation (Hansen and Rattray, 1965). The gravitational circulation is suppressed when the vertical mixing is strong, and vice versa. Moreover, when surface velocity in the estuary is enhanced artificially, more water flux flows seaward in the upper layer. Due to the continuity, additional seawater should enter the estuary in the lower layer. This process will not only falsely enhance the estuarine circulation, but also affect the exchange capacity and time scale of material transport. In W2005, two wall functions, i. e., the parabolic wall function and the open channel wall function, were introduced into the k − kl scheme. The results show that the former made no improvement, but the open channel function eliminated the abnormally high surface current and produced results similar to those by k − ε, k − ω, etc. (Warner et al., 2005). In this paper, we implemented three stable functions in the k − kl and k − ε. The results show that an abnormally high surface current remains without introducing the correct wall function regardless of the stability functions.
Different stability functions impose distinct influences on the vertical distribution of SSCs (Figure 3b). The k − ε appears less gradient in the vertical direction than k − kl near the upper layer (05 m). Suspended sediment distributes more uniformly in the vertical due to a larger diffusivity coefficient and an asymmetry towards the upper layer. Simulations by k − ε may lead to wide variations in the magnitude of verticalaveraged SSCs. In W2005, ANA and k − ω generated the largest SSC. The k − kl with parabolic wall function and the original MY2.5 produce the lowest values. The k − ε and k − kl with open channel wall functions are in between. Nevertheless, their vertical distributions of SSC display extremely similar patterns. In our simulations for the Yamen Channel, SSC based on ANA and k − kl schemes are relatively large. Results by KC and CB are very close except obvious decrease by KC due to the small turbulent diffusion in the surface layer. The results by CA and BB are smaller. This might be attributed to the relatively large turbulent viscosities of the CA and BB, which weakens the estuarine circulation and reduces upstream sediment transport.
Generally, the k − kl closure scheme depends on wall function, which should be chosen with caution. Several alternate wall functions are proposed to produce an apparently correct parabolic turbulent viscosity profile. Given certain boundary and initial conditions, k − ε produces more intensive vertical mixing than k − kl. On the other hand, the k − kl tends to generate stratification. According to the vertical profile of turbulent parameters, vertical velocity, and SSCs, the CB function based on k − ε generates reasonable patterns. Thus, we apply this parameterization to the model validation and further analysis in the Yamen Channel.
4.2 Model validation in the Yamen Channel
Based on the above sensitivity analysis, the k − ε with CB stability function is chosen to validate the simulated results. Three metrics are used to quantify the agreements between observed and modeled data of water level, velocity, and SSC. A skill that has been widely used is written as (Willmott, 1981; Li et al., 2005):
$$ \text { Skill }=1\frac{\sum\limits_{i=0}^n\leftx_{o i}x_{m i}\right^2}{\sum\limits_{i=0}^n\left(\leftx_{o i}\overline{x_o}\right+\leftx_{m i}\overline{x_o}\right\right)^2} $$ (8) where the subscript'o' and 'm' indicate observed and modeled data respectively, and'n' is the number of total time series. Another index to evaluate the absolute deviation of velocity and suspended sediment is mean absolute error (MAE):
$$ \mathrm{MAE}=\frac{\sum\leftx_{o i}x_{m i}\right}{n} $$ (9) The deviation of the current direction is computed using mean absolute error coefficient (MAEC) defined in equa‐ tion (10). Compared to equation (9), equation (10) empha‐ sizes the contributions in the main directions and weakens the influences that arise from the slow slack flow which might be measured less precisely. The weights are chosen as square velocity amplitudes in the sense of kinetic energy.
$$ \text { MAEC }=\frac{\sum\limits_{i=0}^n\left(x_{o i}x_{m i}\right) u_i^2}{\sum\limits_{i=0}^n u_i^2} $$ (10) where u is the velocity, and others are the same as the above.
We validated the modeled results against the observa‐ tions during April 1112 in 2005 (Figure 4). Though the discrepancy increases slightly toward the downstream, the predictions show good agreement with the observations. We collected the statistics of water elevation with two metrics, MAE and Skill (Table 4). Skills are 0.99, and MAEs are less than 0.08 m at the three stations.
Z2 Z3 Z4 MAE 0.021 0.033 0.085 Skill 0.99 0.99 0.99 The comparisons between the modeled and observed velocities at the three stations are shown in Table 5. MAEs of the velocity are 0.030.09 m/s. As the modeled and observed data are not at the exactly same location, and combined with grid resolution, some large deviations emerge. However, the Skills show good performance. MAECs of current direction are less than 11°, and Skills are larger than 0.9.
A3 A4 A5 0.2 h 0.8 h 0.2 h 0.8 h 0.2 h 0.8 h Velocity MAE 0.072 0.064 0.092 0.054 0.035 0.067 Skill 0.97 0.97 0.93 0.98 0.99 0.98 SSC MAE 0.018 0.048 0.045 0.07 0.016 0.039 Skill 0.91 0.84 0.68 0.63 0.79 0.85 A3 A4 A5 MAEC 10.3 9.4 8.6 Skill 0.93 0.96 0.99 The comparisons of SSCs between the measured and predicted are shown in Figure 6 during the spring tide. Circle ("○") represents the measured and the thin line ("—") is the computed. The SSCs are lower than 0.4 kg/m^{3}; in contrast, MAEs are less than 0.07 kg/m^{3}. We concluded that the trend of SSC is considerably identical. According to the comparisons, the simulations are acceptable.
4.3 Spatial variations of SSC
The SSC fields of 0.2 h and 0.8 h layers during the flood and ebb tide are presented in Figures 6 and 7, respectively. The SSCs to the south of Yamen Bridge are higher than that to the north during the flood tide, which vary from 0.1 to 0.3 kg/m^{3} and are basically lower than 0.1 kg/m^{3} to the north of the YB.
The SSCs of each layer during the flood tide are higher than that during the ebb tide, which can be attributed to the tidal asymmetry in the HM Sea. The mean duration of the ebb tide is longer than that of the flood tide, and tidal currents adversely. In the dry season, the duration of the ebb tide is longer than that of the flood tide for 12 hours. The bidirectional current flows mainly along the deep channel. Due to the tidal asymmetry in the channel, the flow velocity is smaller during the ebb tide than during the flood tide. Therefore, suspended particles tend to settle and deposit to the bed, and fewer sediment particles on the bed can be resuspended correspondingly.
The SSCs near the estuarine bar from the HTM to the middle channel are relatively high. At the same station, SSCs during the flood tide are higher than that during the ebb tide. The SSCs vary from 0.03 kg/m^{3} to 0.14 kg/m^{3} in the HM Sea. Due to the complicated interactions between the freshwater and shelf saline water and the effect of local terrain, the SSCs exhibit obvious variations in the vertical.
4.4 Intratidal variations
According to the above analysis, we chose the results based on k − ε and CB stability function to investigate the temporal variations of turbulent viscosity, turbulent kinetic energy, dissipation, velocity and SSC at site A3, which is located near the thalweg of the upstream channel in the Yamen Channel. The tidal currents exhibit a typical flooddominated asymmetry. At high and low water, the velocity is extremely weak (Figure 8a). Flood current is larger than the ebb current of ~0.2 m/s, which is referred to the tidal velocity asymmetry (Song et al., 2011; Gong et al., 2016). The turbulent viscosity is extremely small at both high and low tide levels due to slack water. On the contrary, it is much greater at the maximum ebb and flood tide, which shows a positive correlation between turbulent viscosity and tidal velocity. The turbulent kinetic energy also shows similar variations. However, there is no obvious trend of dissipation. In W2005, the dissipations based on different schemes are nearly identical except MY2.5. Contrastingly, the results of turbulent kinetic energy are different. Moreover, the comparisons in W2005 are only at a timing of maximum ebb. Nonetheless, the simulated SSCs are considerably consistent except for that by MY2.5. This issue deserves further investigation: why does the abnormally high surface velocity have almost no impact on SSC?
In Figure 8d, the dissipations scarcely change corresponding to the tidal velocity similar to turbulent kinetic energy and turbulent viscosity. At high water levels, the dissipation is the smallest, while the largest value is at maximum flood tide. The turbulent dissipation at maximum ebb tide and low water level fall in between. We speculate that this temporal variation is due to the local instantaneous water depth, which may restrict the development of turbulence, and ultimately determine the turbulent dissipation $ \left(\frac{k^{1 / 2}}{l}\right)$ together with the tidal velocity. At high water levels, due to small velocity and large water depth, the turbulent dissipation is significantly weak. At maximum flood tide, the flow velocity is large, and the water level is close to the mean tide level; therefore, the turbulent dissipation maintains large.
Due to the resuspension and increased sedimentcarrying capacity, SSC is the largest at maximum flood tide. At high water levels, in pace with the decrease of velocity, the sediment gradually settles and the SSC decreases. Settling/ scour lag is an interesting issue calling for special attention. As an agent, the sediment response time can be expected to be approximately equal to the vertical distance which scales the fall off sediment concentration divided by the suspended sediment fall velocity (Friedrichs et al., 1998). Analogously, the sediment response time is about 20 min in the Yamen Channel. Therefore, it can be prudently considered that the lag effect is not significant. After the high water levels, the sediment keeps settling. The maximum ebb velocity is generally less than 0.5 ms^{−}^{1}, which is not strong enough for sediment incipience corresponding to a median particle size of 0.05 mm. Therefore, during the whole ebb period, the SSC would not increase. This process can also be seen in Figure 5, where the period of SSC is about semidiurnal rather than quarterdiurnal. Although the turbulent viscosity is the greatest at maximum flood tide (Figure 8d), a certain vertical gradient of SSC still exists due to the resuspension from the seabed acting as a sediment source. At high water levels, the vertical distribution of SSC is more uniform than that at maximum flood tide. Since the weakened turbulent mixing and the sediment settling result in more suspended sediment loss in the upper layer than in the lower layer, the SSC in the surface layer is slightly smaller consequently. At maximum ebb tide, although the suspended sediment cannot be supplemented by the resuspension, the SSC is uniform in the vertical due to the enhanced turbulent mixing. At the low water level, the SSC is extremely low, and the vertical variation is insignificant.
5 Conclusions
Based on the coupled SedFVCOM system with two kinds of turbulence closure schemes with different stability functions implemented in the General Ocean Turbulence Model (GOTM), we carried out numerical experiments to examine the effects of turbulence parameterization on the temporal and spatial distribution of suspended sediment in a natural tidal channel.
The k−kl scheme with three typical stability functions produces anomalously high surface currents. It can be attributed to the improper wall functions embedded in the original MY2.5 closure scheme, which produces an incorrect length scale and underestimation of the turbulent viscosity near the surface. One should use the original MY2.5 with caution and select the proper wall function. The effects of each stability function on vertical dissipation are extremely small. By contrast, the k−kl and k−ε can produce significantly different turbulent viscosity in the vertical, and thus turbulent diffusivity. We concluded that under specified boundary and initial conditions, k−ε turbulence schemes generally produce more intensive vertical mixing than the k−kl does.
Allowing for the simultaneous depositional process during erosion near the watersediment interface, we applied k−ε with CB stability function in the Yamen channel, one of the eight outfalls of the Pearl River. During a tidal period, the turbulent dissipation show highly similar trends while the patterns of the SSCs differ significantly. Although velocity, turbulent viscosity, and TKE keep similar trends and magnitudes respectively, the SSCs show different patterns between low and high water due to complex mechanisms.
This paper illuminates the influence of two kinds of turbulent closure schemes with several different wall functions on velocity and SSCs in the Yamen channel, which is important to the GuangdongHong KongMarco Greater Bay Area. Our findings provide insights into, and references for, turbulent and sediment parameterization for researchers and estuarine management.
Competing interestThe authors have no competing interests to declare that are relevant to the content of this article. 
Table 1 Comparison between the watersediment distribution and tidal range between YM and HTM at the HM Sea
Items Branches HTM YM Total amount of 8 branches Annual discharge (10^{8} m^{3}) 202.0 196.0 3 260.0 Percent (%) (the discharge of the Pearl River) 6.2 6.0 100.0 Annual sediment transport rate (10^{4} t) 509.0 363.0 8 872.0 Percent (%) (the sediment of the Pearl River) 7.2 5.1 100.0 Annualaveraged tidal range (m) 1.20 1.24 The maximum tidal range 2.66 2.95 Annual flood tide volume (10^{8} m^{3}) 57.0 636.0 3 763.0 The ratio of total runoff and flood tide volume 3.54 0.31 Table 2 Turbulence parametrization including the stability functions
Turbulent scheme Stability functions Parameterizations Isotropy Shear production Buoyancy Production Nonisotropic Vorticity k − ε BB √ √ CA √ √ √ √ √ CB Deduced by the quadrivium of turbulent dissipation and production. k − kl MY2.5 Functions of the gradient Richardson number. G1988 √ √ KC √ √ Table 3 Parameters of the hydrodynamic and sediment model
Horizontal resolution (m) Nodes Cells Time step (s) Time step (s) τ_{ce} (N/m^{2}) Erosion constant (kg/m^{2}s) 2080 3 563 6 295 0.1 0.75 0.05 5×10^{5} Table 4 jmsa2222844
Z2 Z3 Z4 MAE 0.021 0.033 0.085 Skill 0.99 0.99 0.99 Table 5 Statistics of tidal velocity (m/s) and SSC (kg/m^{3})
A3 A4 A5 0.2 h 0.8 h 0.2 h 0.8 h 0.2 h 0.8 h Velocity MAE 0.072 0.064 0.092 0.054 0.035 0.067 Skill 0.97 0.97 0.93 0.98 0.99 0.98 SSC MAE 0.018 0.048 0.045 0.07 0.016 0.039 Skill 0.91 0.84 0.68 0.63 0.79 0.85 Table 6 Statistics of vertical average flow direction
A3 A4 A5 MAEC 10.3 9.4 8.6 Skill 0.93 0.96 0.99 
Blumberg AF, Galperin B, O'Connor DJ (1992) Modeling vertical structure of openchannel flows. Journal of Hydraulic Engineering, 118(H8), 11191134. https://doi.org/10.1061/(ASCE)07339429(1992)118:8(1119) Burchard H, Baumert H (1995) On the performance of a mixedlayer model based on the kε turbulence closure. Journal of Geophysical Research 100, 85238540. https://doi.org/10.1029/94JC03229 Burchard H, Bolding K (2001) Comparative analysis of four secondmoment turbulence closure models for the oceanic mixed layer. Journal of Physical Oceanography, 31, 19431968. https://doi.org/10.1175/15200485(2001)031<1943:CAOFSM>2.0.CO;2 Burchard H, Craig PD, Gemmrich JR, van Haren H, Mathieu PP, Meier HEM, Smith WAMN, Prandke H, Rippeth TP, Skyllingstad ED, Smyth WD, Welsh DJS, Wijesekera HW (2008) Observational and numerical modelling methods for quantifying coastal ocean turbulence and mixing. Progress in Oceanography 76, 399442. https://doi.org/10.1016/j.pocean.2007.09.005 Burchard H, Gräwe U, Holtermann P, Klingbeil K, Umlauf L (2014) Turbulence Closure Modelling in Coastal Waters. Die Küste 81, 6987 Burchard H, Petersen O, Rippeth T (1998) Comparing the performance of the MellorYamada and the kε twoequation turbulence models. Journal of Geophysical Research, 103, 1054310554. https://doi.org/10.1029/98JC00261 Canuto VM, Howard A, Cheng Y, Dubovikov MS (2001) Ocean turbulence. Part Ⅰ: onepoint closure modelmomentum and heat vertical diffusivities. Journal of Physical Oceanography 31, 14131426. https://doi.org/10.1175/15200485(2001)031<1413:OTPIOP>2.0.CO;2 Chen C, Beardsley RC (1998) Tidal mixing and crossfrontal particle exchange over a finite amplitude asymmetric bank: a model study with application of Georges Bank. Journal of Marine Research 56, 11631201. https://doi.org/10.1357/002224098765093607 Chen X, Chen Y, Lai G (2003) Modeling of the transport of suspended solids in the estuary of Zhujiang River. Acta Oceanological Sinica, 25(2), 120127. https://doi.org/10.1353/crt.2004.0011 Chen C, Huang H, Beardsley RC, Liu H, Xu Q, Cowles G (2007) A finitevolume numerical approach for coastal ocean circulation studies: comparisons with finite difference models. Journal of Geophysical Research, 112, https://doi.org/10.1029/2006JC003485 Choi Y, Park Y, Choi M, Jung KT, Kim KO (2021) A fine grid tidewaveocean circulation coupled model for the Yellow Sea: comparison of turbulence closure schemes in reproducing temperature distributions. Journal of Marine Science and Engineering 9(21): 1460. https://doi.org/10.3390/jmse9121460 Costa A, Doglioli AM, Marsaleix P, Petrenko AA (2017) Comparison of in situ microstructure measurements to different turbulence closure schemes in a 3D numerical ocean circulation model. Ocean Modelling 120: 117. https://doi.org/10.1016/j.ocemod.2017.10.002 Ding Y, Chen C, Beardsley RC, Bao X, Shi M, Zhang Y, Lai Z, Li R, Lin H, Viet NT (2013) Observational and model studies of the circulation in the Gulf of Tonkin, South China Sea. Journal of Geophysical Research: Oceans 118, 64956510. https://doi.org/10.1002/2013JC009455 Egbert GD, Erofeeva SY (2002) Efficient inverse modeling of barotropic ocean tides. Journal of Atmospheric and Oceanic Technology 19(2), 183204, 2002. https://doi.org/10.1175/15200426(2002)019<0183:EIMOBO>2.0.CO;2 Friedrichs CT, Armbrust BD, Swart HED (1998) Hydrodynamics and equilibrium sediment dynamics of shallow, funnelshaped tidal estuaries. Physics of Estuaries and Coastal Seas, Dronkers & Scheffers (Eds). Balkema, Rotterdam, 315327 http://web.vims.edu/~cfried/cv/1998/Friedrichs_etal_1998_PECS_Volume.pdf Galperin B, Kantha LH, Hassid S, Rosati A (1988) A quasiequilibrium turbulent energy model for geophysical flows. Journal of the Atmospheric Sciences 45(1), 5562. https://doi.org/10.1175/15200469(1988)045<0055:AQETEM>2.0.CO;2 Geyer WR, Ralston DK (2015) Estuarine frontogenesis. Journal of Physical Oceanography, 45(2): 546561. https://doi.org/10.1175/JPOD140082.1 Gong W, Zhang H, Schuttelaars HM (2016) Tidal asymmetry in a funnelshaped estuary with mixed semidiurnal tides. Ocean Dynamics 66: 637658. https://doi.org/10.1007/s1023601609431 Hansen DV, Rattray M (1965) Gravitational Circulation in Straits and Estuaries. Journal of Marine Research 23(2): 104122. https://doi.org/10.1357/002224021834614399 Kantha LH, Clayson CA (1994) An improved mixed layer model for geophysical applications. Journal of Geophysical Research, 99, 2523525266. https://doi.org//10.1029/94JC02257 Lai Z, Ma R, Gao G, Chen C, Beardsley RC (2015) Impact of multichannel river network on the plume dynamics in the Pearl River estuary. Journal of Geophysical Research: Oceans, 120(8), 57665789. https://doi.org/10.1002/2014JC010490 Li R, Luo F, Zhu W (2008) The suspended sediment transport equation and its nearbed sediment flux. Science in China Series E: Technological Sciences, 52: 387391. https://doi.org/10.1007/s1143100801759 Li M, Zhong L, Boicourt WC (2005) Simulations of Chesapeake Bay estuary: Sensitivity to turbulence mixing parameterizations and comparison with observations. Journal of Geophysical Research Oceans 110(12): 347356. https://doi.org/10.1029/2004JC002585 Lima HP, Dias FJS, Teixeira C, GodoiVA, Torres Jr. AR, Araujo RS (2021) Implications of turbulence in a macrotidal estuary in northeastern BrazilThe São Marcos Estuarine Complex. Regional Studies in Marine Science, 47. https://doi.org/10.1016/j.rsma.2021.101947 Liu J, Yuan J, Liang J (2022) An evaluation of vertical mixing parameterization of ocean boundary layer turbulence for cohesive sediments, Deep Sea Research Part Ⅱ: Topical Studies in Oceanography 204. https://doi.org/10.1016/j.dsr2.2022.105168 Ralston DK, Cowles GW, Geyer WR, Holleman RC (2017) Turbulent and numerical mixing in a salt wedge estuary: dependence on grid resolution, bottom roughness, and turbulence closure. Journal of Geophysical Research: Ocean 122: 692712. https://doi.org/10.1002/2016JC011738 Reffray G, Bourdallé Badie R (2015) Modelling turbulent vertical mixing sensitivity using a 1D version of NEMO. Geoscientific Model Development, 8: 6986. https://doi.org/10.5194/gmd8692015 Reissmann J, Burchard H, Feistel R, Hagen E, Lass H, Volker, Mohrholz, Nausch G, Umlauf L, Wieczorek G (2009) Stateoftheart review on vertical mixing in the Baltic Sea and consequences for eutrophication. Progress in Oceanography, 82: 4780. https://doi.org/10.1016/j.pocean.2007.10.004 Rodi W (1984) Turbulence models and their application in hydraulics. A state of the art review. United States: N.P., Web Song D, Wang XH, Kiss AE, Bao X (2011) The contribution to tidal asymmetry by different combination of tidal constituents. Journal of Geophysical Research, 116: 112. https://doi.org/10.1029/2011JC007270 Sun Y, Chen C, Beardsley RC, Ullman D, Butman B, Lin H (2016) Surface circulation in Block Island Sound and adjacent coastal and shelf regions: A FVCOMCODAR comparison. Progress in Oceanography, 143: 2645. https://doi.org/10.1016/j.pocean.2016.02.005 Tu J, Fan D, Zhang Y, Voulgaris G (2019) Turbulence, sedimentinduced stratification, and mixing under macrotidal estuarine conditions (Qiantang Estuary, China). Journal of Geophysical Research: Oceans, 124, 40584077. https://doi.org/10.1029/2018JC014281 Umlauf L, Burchard H (2011) Diapycnal transport and mixing efficiency in stratified boundary layers near sloping topography. Jouranl of physical oceanography, 41, 329345. https://doi.org/10.1175/2010JPO4438.1 Umlauf L, Burchard H, Bolding K (2003) Extending the kω turbulence model towards oceanic applications. Ocean Modeling, 5, 195218. https://doi.org/10.1016/S14635003(02)000392 Van Rijn (2005) Principles of sedimentation and erosion engineering in rivers, estuaries and coastal seas. AQUA Publications Warner JC, Sherwood CR, Arango HG, Signell RP (2005) Performance of four turbulence closure models implemented using a generic length scale method. Ocean Modelling, 8(12): 81113. https://doi.org/10.1016/j.ocemod.2003.12.003 Warner JC, Sherwood CR, Signell RP, Harris CK, Arango HG (2008) Development of a threedimensional, regional, coupled wave, current, and sedimenttransport model. Computers & Geoscience 34: 12841306. https://doi.org/10.1016/j.cageo.2008.02.012 Willmott CJ (1981) On the Validation of Models. Physical Geography, 2: 184194. https://doi.org/10.1080/02723646.1981.10642213 Wu H, Wang Y, Gao S, Xing F, Tang J, Chen D (2022) Fluid mud dynamics in a tidedominated estuary: A case study from the Yangtze River. Continental Shelf Research 232, https://doi.org/10.1016/j.csr.2021.104623 Yang X, Wang W (1994) A statistical study on hydrodynamic sedimentation and sand movement Huangmaohai Estuary Bay (in Chinese). The Ocean Engineering, 12(4): 4551 http://en.cnki.com.cn/article_en/cjfdtotalhygc404.005.htm Yang Y, Guan W, Deleersnijder E, He Z (2022) Hydrodynamic and sediment transport modelling in the Pearl River Estuary and adjacent Chinese coastal zone during Typhoon Mangkhut. Continental Shelf Research 233. https://doi.org/10.1016/j.csr.2022.104645 Zhang G, Cheng W, Chen L, Zhang H, Gong W (2019) Transport of riverine sediment from different outlets in the Pearl River Estuary during the wet season. Marine Geology, 415, https://doi.org/10.1016/j.margeo.2019.06.002