Sensitivity analysis of a three-dimensional simulation of turbidity currents in a sloping flume

Sensitivity analysis is necessary in numerical models to establish the optimal model configuration and simulate turbidity currents accurately. The sensitivity of turbidity currents simulated by a three-dimensional numerical model named the semi-implicit cross-scale hydroscience integrated system model (SCHISM) was assessed using domain discretization, the turbulence closure model, and the transport scheme. The results show that a higher horizontal resolution with an appropriate aspect ratio, localized sigma coordinates with shaved cells (LSC2) grid with sufficient layers, a two-phase mixture turbulence closure model, and transport schemes with critical depth ratio lower than 1 could improve the model performance on the turbidity currents. The study can provide helpful guidance to establish accurate turbidity current simulations in complex water environments.


Introduction
When muddy water flows into stagnant water, it usually plunges to the bottom due to the presence of sediment solids with heavier density than the ambient water. This common phenomenon is called a turbidity current, occurring when muddy water enters an ocean, estuary, reservoir or lake (Choi, 1998;Garcia & Parker, 1993;Huang et al., 2009;Luchi et al., 2018;Nourmohammadi et al., 2011;Qian & Wan, 1983;Simpson, 1997). The turbidity current is the main pattern for sediment transport, entrainment and deposition in reservoirs, significantly affecting the reservoir operation and environment.
Numerical models have become an attractive tool to study turbidity currents with the rapid development of computational fluid mechanics. Numerical models can be divided into layeraveraged and depth-resolving models. Layer-averaged models include single layer-averaged models (Lai & Wu, 2013;Parker et al., 1986;Winslow, 2001) and double layer-averaged models (Bonnecaze et al., 1993;Cao et al., 2015;La Rocca et al., 2012), and can briefly describe the propagation speed and deposition patterns of turbidity currents. Depth-resolving models include three-dimensional and vertical two-dimensional models, and can capture the formation and evolution of turbidity currents, as well as the vertical structure of velocity and density (An & Julien, 2014;Bournet et al., 1999;De Cesare et al., 2001;Kassem & Imran, 2001;Kassem et al., 2003).
In three-dimensional (3D) models, the vertical structure of the turbidity current is mainly solved by different degrees of simplification of the Navier-Stokes equation, such as applying Reynolds approximation. A turbulence closure model (TCM) coupled with wall functions is used to estimate the Reynolds stress (Pérez-Díaz et al., 2019). A mixing-length model was used for TCM by Bowen (1988a, 1988b), and other TCMs, such as the k-ε model, have been widely used (Bournet et al., 1999;Choi & Garcia, 2002;Farrell & Stefan, 1988). Warner et al. (2005) tested the performance of the original k-kl, k-ε and k-ω models, and new models suggested by Umlauf and Burchard (2003) for an oscillatory stratified pressure-gradient driven flow in a rectangular channel. The salt intrusion and salinity distribution of the original k-kl model were distinctly different from those of the others. Abd El-Gawad (2011), Li et al. (2005) and Ye et al. (2018) obtained similar responses of the k-kl, k-ε and k-ω models for vertical stratification simulations. Recently, a new turbulence closure model named the two-phase mixture turbulence model (TPM model), derived from the two-phase mixture theory, was proposed for open channel sediment-laden flows (Huang et al., 2019;Zhong et al., 2011Zhong et al., , 2014. However, the TPM model has not been applied to turbidity currents, and the differences between the TPM model and the previous TCMs have not been discussed. In addition, the previous sensitivity analyses paid more attention to the effects of TCMs on the stable stratification after the plunging process, but neglected the potential effects on the plunging process of turbidity currents. Another key influencing factor is the advection scheme for the tracer, which affects the accuracy of the vertical structure of turbidity currents. The improved upwind and 2nd-order total variation diminishing (TVD) schemes have been widely used in numerical models (Bourban, 2013;Kurganov & Petrova, 2007;Luchi et al., 2018;Yeh et al., 2013). The 2nd-order TVD scheme treats all fluxes explicitly, especially vertical fluxes with smaller grid sizes in the vertical dimension. Therefore, a large number of sub-cycles are required with each time step, and TVD is much more expensive than upwind schemes. Ye et al. (2018) found that upwind and explicit TVD schemes under-predicted intrusion and stratification in Chesapeake Bay. To address this, they introduced an implicit transport solver in time and space with two total variation diminishing limiters (TVD 2 ) of second-order accuracy, showing that this scheme could improve the simulation results in saltwater intrusion and density stratification. To enhance the ability of the model to resolve the baroclinic circulation in an eddy ocean, a third-order weighted essentially non-oscillatory (WENO) transport scheme for horizontal transport was applied and shown to outperform the TVD 2 scheme for better control of numerical dissipation, maintenance of sharp frontal instabilities, and better resolution of baroclinic circulation and eddy regime in the Gulf Stream (Ye et al., 2019). However, the performance of the TVD 2 and WENO schemes on the plunging and vertical profiles of the turbidity currents in the reservoirs has not been investigated, which leads to a lack of advice on the appropriate transport scheme for 3D numerical simulation of turbidity currents.
The computational grid should be designed appropriately in horizontal and vertical dimensions to simulate real physical processes, so that artificial numerical effects can be avoided, especially for applications involving steep bathymetry and strong stratification. The impacts of grid resolution on the propagation rate, dilution and vertical profiles of density currents were analysed by Choi (1999) and Pérez-Díaz et al. (2019). These studies only focused on the effects of grid resolutions on the density currents with stable stratification because they used submerged jets to create the inflow. Therefore, for turbidity currents in reservoirs with obvious plunging processes, the effects of grid resolution on the plunging points are still unclear.
Due to the limited focus on the effects of the TCM, the transport scheme and domain discretization on turbidity currents in reservoirs with a plunging process, this study aims to establish rules for simulating this kind of turbidity current by sensitivity analysis. A 3D numerical model named the semi-implicit cross-scale hydroscience integrated system model (SCHISM) is utilized to simulate the laboratory turbidity currents in a sloping flume of Lee and Yu (1997) under different conditions. The sensitivity analyses of the TCM, transport scheme and domain discretization are conducted by comparing the plunging points, propagation speed, and vertical profiles between the simulated and measured turbidity currents. The proposed rules can provide advice for more turbidity current simulations under different conditions by SCHISM.
In the following sections, a description of the SCHISM model and the analytical methods are presented first. Second, the sensitivity tests for domain discretization, TCM, and transport scheme are discussed to propose the appropriate advice for model set-up. Third, the performance of the proposed model set-up to reproduce the turbidity currents is assessed. Finally, the model setup rules of the turbidity currents based on the sensitivity analysis are summarized.

Model description
SCHISM is a widely used 3D numerical model and is distributed as an open-source community-supported model mostly via github.com (Zhang et al., 2004). SCHISM solves Reynoldsaveraged Navier-Stokes (RANS) equations in hydrostatic form with the Boussinesq approximation, using a hybrid finite element and finite volume approach.
The continuity equation is: The momentum equation is: The sediment transport equation is: and the free surface equation is derived from the surface tracking method: where (x, y) represent horizontal Cartesian coordinates; z represents the vertical coordinate with positive upward; (u, v, w) represent velocity in three directions; t is time; η is the free-surface elevation; h is the bathymetric depth; p A is the atmospheric pressure at the free surface; ρ 0 is reference water density; ρ is the water density considering temperature (T), hydrostatic pressure (p) and sediment concentration (C); g is the acceleration of gravity; K mv is the vertical eddy viscosity coefficient, which is closed with the TCM; K sv is the vertical eddy diffusivity for tracers; ω s is the grain settling velocity; and F mx , F my and F h are the horizontal diffusion for momentum and transport equations. SCHISM is grounded on unstructured grids in the horizontal dimension and provides two different vertical grids, including a hybrid SZ grid (i.e. terrain-following generated S-and shaved Z-coordinates) (Song & Haidvogel, 1994;Zhang & Baptista, 2008) and a highly flexible hybrid "localized sigma coordinates with shaved cells" (LSC 2 ) grid (Zhang et al., 2015). The LSC 2 grid can be generated by the vanishing quasi sigma (VQS) of Dukhovskoy et al. (2009) to define specific vertical layers for each horizontal node to minimize non-smoothness among neighbouring nodes (Zhang et al., 2016a).
Equations (1-6) are closed by the TCM named the generic length-scale model (GLS, Umlauf & Burchard, 2003) or the new TPM model (Huang et al., 2019). The GLS model not only includes the k-ε (Rodi, 1984), MY25 (original k-kl, Mellor & Yamada, 1982), and k-ω (Wilcox, 1988) models but also includes the new k-kl model developed from MY25 by Warner et al. (2005) and the UB model generated by Umlauf and Burchard (2003). The TPM model is established based on the k-ε model (Huang et al., 2019). The GLS model contains two standard equations for transport k and a generic parameter ψ: where σ k and σ ψ are the Schmidt numbers for k and ψ; P and B represent production by shear and buoyancy, respectively; F wall is a wall proximity function; c 1 and c 2 are coefficients selected to be consistent with von Karman's constant and with experimental observations for decaying homogeneous, isotropic turbulence (Wilcox, 1988); c 3 is associated with the buoyancy term and approaches the value of c -3 in stably stratified flows, and c + 3 in unstably stratified flows. The value of c -3 depends on the stability function of Kantha and Clayson (1994) by a relationship of c − 3 = 5.08c 1 − 4.08c 2 for the k-ε, k-kl, k-ω and UB models. The c -3 values of the TPM and MY25 models are specifically defined in Huang et al. (2019) and Mellor and Yamada (1982). The value of c + 3 was set to 1 according to Warner et al. (2005), except for the c + 3 of MY25, where c + 3 = 0.9 (Mellor & Yamada, 1982). The generic parameter ψ, the dissipation rate ε and the frequency of the turbulence decay process ω (Saffman, 1970) are defined as: where c 0 μ is the stability coefficient, and p, m, and n are modeldefined constants for different TCMs. Table 1 lists the values of the parameters in Eqs (7-10) and the corresponding forms of the generic parameter ψ for the k-ε, k-kl, MY25, k-ω, UB, and TPM models (Huang et al., 2019;Umlauf & Burchard, 2003;Warner et al., 2005).
The terms with stringent Courant-Friedrichs-Lewy (CFL) condition are treated implicitly and the implicit time-stepping scheme is adopted to obtain efficiency and robustness. The Eulerian-Lagrangian method (ELM) is used for momentum advection to balance the numerical diffusion and dispersion (Stanev et al., 2018). The TVD 2 and WENO methods are used to reach the 2nd-or 3rd-order accuracy for the advection schemes for transport equations, respectively (Ye et al., 2019;Zhang et al., 2016b).
The suspended sediment transport for non-cohesive sediment in SCHISM is coupled with the regional ocean modelling system (ROMS; Warner et al., 2008). The model implements this algorithm to calculate an unlimited number of non-cohesive sediment size classes. The mass balance between the suspended sediment and bedload is controlled by erosion and deposition algorithms. A multi-level bed framework is set up to track the mass, fractions, thickness, and surface properties of the sediment to characterize the morphology and stratigraphy.

Experimental database and model set-up
The laboratory experiments of turbidity currents by Lee and Yu (1997) were selected for this study because turbidity currents could be generated with various initial conditions, and sufficient Table 1 Parameters for different turbulence closure models   Lee and Yu (1997) investigated the migration of plunge points, the vertical profiles of velocity and sediment concentration, and the longitudinal variations in the thickness, depth-averaged velocity, and depth-averaged concentration. The inflow conditions of the turbidity currents and the stable plunge locations are described in Table 2, where q 0 is the inlet discharge per unit width, C 0 is the inflow sediment concentration, T is the temperature of the inflow and the ambient water, and X p and h p are the plunge locations from the inlet and plunge depths, respectively. The initial and boundary conditions of the SCHISM applications were defined to numerically reproduce the real experimental set-up. The sand used in the experiments has a specific gravity of 2.65 and a mean particle size of 0.0068 mm. The initial water level of the model was set to maintain the 4 m length of the open channel region, as shown in Fig. 1, and the initial velocity field and suspended sediment concentration field for the computational domain were both set to zero. The initial temperature of the stagnant water and the inflow was the same, so the effects of temperature on the density stratification were neglected in this study.
At the free-surface boundary, SCHISM considers zero gradients perpendicular to the boundary when the effects of wind stress are neglected in rivers or reservoirs. At the bottom Table 2 Initial inflow conditions and plunge locations of different experimental cases from Lee and Yu (1997) No. boundary, the no-slip condition is replaced by a balance between the internal Reynolds stress and the bottom frictional stress due to the presence of the bottom boundary layer. The bottom stress here depends on the turbulent boundary layer (Blumberg & Mellor, 1987), in which the velocity profile obeys the logarithmic law and is matched to the exterior flow at the top of the boundary layer. The vertical boundary conditions of the transport equation ensure consistency with the source/sinks. At the free surface and at the bottom of rivers, the turbulent kinetic energy κ and the mixing length l are specified as Direchlet boundary conditions.

Analytic approach
The TC14 case in Table 2 was chosen as the fixed inflow condition for the sensitivity analysis. This analysis includes a series of simulations with different key numerical aspects, including domain discretization, TCM selection, and tracer advection schemes. In different simulations, one numerical aspect was varied while the rest were unchanged, as summarized in Table 3. After determination of the most appropriate rules of the model set-up based on the sensitivity analysis, the other cases in Table 2 were all simulated by the proposed model configuration to evaluate the effectiveness. The optimal numerical aspects were assessed and determined by comparing the model and measured results of the plunging locations of X p and h p , the vertical profiles of velocity (U) and sediment concentration (C) between the bed and the interface of the turbidity current and ambient water, the thickness (h tc ), the depth-averaged velocity (U tc ) and the depth-averaged sediment concentration (C tc ), as shown in Fig. 2. In addition, the longitudinal profiles of the maximum velocity (U m ) and the averaged sediment concentration of the denser layer (C m ) in Fig. 2 and the propagation speed of the turbidity current head were also investigated. The values of h tc , U tc , C tc and Reynolds number (Re tc ) were computed according to Ellison and Turner (1959): where δ is the distance between the channel bed and the position where the velocity equals zero and ν t is the kinematic viscosity of water.
3 Sensitivity analysis of model parameters

Horizontal discretization
Considering the rectangular flume shape, the horizontal domain can be discretized with quadrangular meshes of different sizes, as given in Table 4, where x and y represent the grid size parallel to and perpendicular to the flow direction, respectively. The vertical discretization was temporarily fixed according to the setting of Zhang and Wu (2020), with 51 layers in the maximum depth region and four layers in the open channel region. The drag coefficient was set to 0.01, in the range of the estimated drag coefficient of 0.007-0.02 based on the measured velocity profiles in Lee and Yu (1997) by the calculation method of Garcia and Parker (1993). Table 4 shows the plunge locations and depths of the model results with different horizontal grid sizes. Comparing the X p and h p of H01-H05 with the same y but different x values, the low resolution of x (H05) resulted in farther plunge locations and larger plunge depths. At a sufficiently high resolution of x, the plunge points of H01-H03 converged to a stable position, whose X p and h p were close to the measured values in Table 2. As shown in Table 4, H04, H06, H07 and H08 had the same x but different y values, resulting in almost the same X p and h p values, which indicated that the resolution of y had few effects on the plunge points. Figure 3 shows the propagation speed of the turbidity current, and the longitudinal profiles of U m and C m at 210 s, when the turbidity currents reached the tail of the flume. As shown in Fig. 3a, turbidity currents with different horizontal grids initially exhibited distinctly different speeds, and the higher resolution could cause the lower propagation speed of H01-H03. As shown in Fig. 3b, only the U m of H02 and H03 converged to the same profiles. As shown in Fig. 3c, the C m of H01-H05 presented similar profiles, showing the limited effects of x on the sediment concentration. However, the C m of H04 and H06-H08 decreased along the flume more quickly with lower resolution, showing the obvious effects of y on the sediment concentration.
Considering the effects of horizontal resolutions on the plunging and transport process of turbidity currents together, only H02 and H03 could obtain the converged plunge depths, propagation speed of turbidity current, and longitudinal profiles of U m and C m , which indicated that the highest horizontal resolution of H01 did not present the best results for the turbidity   currents. Instead, the resolution of H02 or H03 was more appropriate for the turbidity current simulation in this study, with x/ y between 4 and 5. Therefore, H02 was chosen as the horizontal discretization size to simulate the turbidity current, allowing the balance between the plunging and transport processes of the turbidity currents with relatively high accuracy and efficiency.

Vertical discretization
For turbidity currents with density stratification, the vertical structure will be sensitive to the vertical grids and resolution. To test the sensitivity of the two vertical grids of the SZ and LSC 2 methods on the turbidity currents in the reservoir, cases of different vertical grid systems and layers were designed, as shown in Table 5. Cases V01-V06 were designed by the LSC 2 method, with four layers in the open channel region and more layers in the maximum depth region. Cases V07-V09 were designed by the pure S grid of the SZ method of θ f = 0.001, with the same layers from the open channel to the maximum depth region. The vertical depth for each case could be transformed to σ -coordinates with the range from − 1 to 0, meaning from the bottom to the surface. The σ values between the adjacent vertical layers in Table 5 indicate the thickness of each layer. Table 5 also shows the plunge locations and depths of the model results with different vertical grid sizes. The SZ grids with different layers have few effects on the plunge points, because the plunge points remained at a stable location approximately 11 m from the inlet of V07-V09. However, the plunge points were sensitive to the resolution of the LSC 2 grid of V01-V06. With a lower resolution of the LSC 2 grid, the plunge point was farther downstream and the plunge depth was larger. When the vertical layers were sufficient for cases V01-V03, the plunge points converged to a stable location approximately 11 m from the inlet. As shown in Fig. 4a, V01-V03 and V07-V08 exhibited almost converged propagation speed, but the other cases showed lower propagation speed due to the lower vertical resolution. In Fig. 4b, the abrupt decrease in the surface sediment concentration indicated the location of the plunge points, and the longitudinal profiles of V01-V06 before the plunge points were much smoother than those of V07-V09. The severe fluctuation of the surface sediment concentration of V07 and V08 before plunging could be due to too many vertical layers in the shallow region and extremely thin computational cell thickness, which caused the thin layer syndrome of small time steps to satisfy the CFL condition, inappropriate pressure gradient computations, exaggeration of the drag coefficient and so on (Yamazaki & Satomura, 2010;Zhang et al., 2015). Therefore, even though V07 and V08 showed converged propagation speed and stable plunge points, they were not suitable choices.
The U tc , C tc and h tc values computed by Eqs (12-14) were used to normalize the velocity and sediment concentration profiles in terms of U/U tc and C/C tc versus z/h tc . The normalized profiles of the velocity and sediment concentration at 180 s and a cross section of 15 m from the inlet are presented in Fig. 5. Cases V01-V03 converged to the same velocity and concentration profiles in Fig. 5a and 5b with enough vertical layers. For fewer vertical layers of V04-V06, the velocity profiles showed sharper forms with larger velocity peaks and lower locations, and the sediment concentration profiles showed thinner denser layers. The results revealed that V04-V06 were unable to accurately capture the real profiles of velocity and sediment concentration. The effects of different layers of SZ grids on the vertical structures of the velocity and sediment concentration in Fig. 5c and 5d were similar to those of LSC 2 grids.
Based on the plunge points, propagation speed, longitudinal profiles of the surface sediment concentration, and dimensionless vertical profiles of velocity and sediment concentration, V01-V03 were all considered as the ideal vertical discretization for the model. Considering the accuracy and the computational efficiency, V03 was selected as the vertical grid for the subsequent discussion.

Turbulence closure modelling
Turbulence mixing plays a significant role in determining the stratification and residual circulation in estuaries. This section presents the comparison between the applications of wellknown TCMs, including k-ε, k-kl, MY25, k-ω, UB and TPM, as introduced in the previous section. Based on the horizontal grids of H02 and vertical grids of V03, case TC14 was simulated with different TCMs.  Zhang et al. Journal of Hydraulic Research Vol. 61, No. 1 (2023)  The plunge depths and locations of the model results with different TCMs are shown in Table 6. Two methods of plunging criteria, zero-velocity and half-concentration, were used to estimate the location of the plunge point. For the zerovelocity method, the stable plunge point was defined as the point with zero streamwise velocity on the surface, and the interface between the turbidity current and the ambient water was also defined by zero streamwise velocity, as shown in Fig. 6a. According to this method, the plunge depths of UB, k-kl, and MY25 were much smaller than those of k-ε, TPM, and k-ω in Table 6. For the half-concentration method, the plunge point was defined where the surface sediment concentration decreased to half of the inlet sediment concentration, approximately 4 kg m -3 of TC14, as shown in Fig. 6b. According to this method, the plunge depths of different TCMs were slightly different. The plunge locations of the UB, k-kl and MY25 models by the half-concentration method were obviously farther than the locations by the zero-velocity method. The plunge locations of the k-ε, TPM and k-ω models remained stable by two different criteria, and were closer to the measured data in Table 2 than those of the other models.
The normalized profiles of the velocity and sediment concentration with different TCMs are presented in Fig. 7, showing that the vertical profiles were not extremely sensitive to the TCMs, as observed by Li et al. (2005) and Ye et al. (2018). For velocity profiles, the U m /U tc were all approximately 1.2, but the corresponding heights h um /h tc were slightly different. For concentration profiles, the h cm /h tc of M02, M04, M05 and M06 were similar, slightly thinner than the h cm /h tc values of M01 and M03. The longitudinal profiles of sediment concentration at 200 s of different TCMs are presented in Fig. 8. Figure 8a, 8b and 8c show similar longitudinal profiles of sediment concentration, with a concave profile around the plunge zone, consistent with the experiments of Lamb et al. (2010) and Sequeiros et al. (2009), and the model results of Schuch et al. (2018). However, Figure 8d, 8e and 8f show a gradual change in the plunge zone with a convex shape, similar to the observations of An and Li (2010) and Sequeiros et al. (2018) and the model results of Cao et al. (2015). Although the two different shapes of the plunge zone were both observed, the longitudinal profiles of M01-M03 were more stable with a constant thickness and smooth shapes along the flume than M04-M06. It should be noted that the free surface in Fig. 8 presented an unnatural oscillation due to the different scales of the x and z coordinates, which exaggerated the normal fluctuation in the vertical direction. Based on the analysis of the plunge points and longitudinal profiles, the k-ε, TPM, and k-ω models were all recommended for use in this model, and TPM was finally selected for subsequent s analysis.

Transport schemes
To reduce the computational cost of the model, a critical depth ratio is defined as the cut-off depth divided by the depth from the open channel region to the backwater region in SCHISM to control the transition between the upwind and TVD 2 or WENO schemes. In an area with no strong stratification, an upwind scheme is encouraged (Zhang et al., 2016b). When the depth ratio is larger than the critical depth ratio, the TVD 2 and WENO methods with higher accuracy are used. Based on the selected domain discretization and TPM model, TC14 was simulated and compared under the TVD 2 and WENO methods with different critical depth ratios. The effects of different transport schemes and critical depth ratios on the plunge points are presented in Table 7, indicating that the plunge depths of the WENO scheme were all smaller than those of the TVD 2 scheme. The smaller critical depth ratio caused the smaller plunge depth, but only limited differences were observed for T01-T08, which could be neglected.
As shown in Fig. 9a, the velocity profiles were not sensitive to the transport schemes and critical depth ratios. However, the sediment concentration profiles were sensitive to the critical depth ratios as shown in Fig. 9b. When the critical depth ratio was higher than 1 of T01, T02, T05 and T06, the sediment concentration profiles exhibited a smooth transition from the denser layer to the interface at z/h tc = 0.6-0.8. When the critical depth ratio was lower than 1 for T03, T04, T07 and T08, the transition was sharper and the denser layer was more obvious, which better captured the characteristics of the measured sediment concentration profiles (Lee & Yu, 1997). It can be concluded that TVD 2 and WENO can both be used in this model to accurately predict Journal of Hydraulic Research Vol. 61, No. 1 (2023)   the turbidity current in a sloping flume when the critical depth ratio is lower than 1.

Model validation and assessments
Based on the aforementioned sensitivity analysis, the optimal choices of the numerical aspects for the turbidity current simulation in a sloping flume are listed in Table 8. The calculated vertical profiles of velocity and sediment concentration of TC14 by the optimal set-ups were compared with the observed profiles in Fig. 10a and 10b. The high agreement between the simulated and observed profiles indicated the superiority of the proposed model configuration. To further assess the performance of the proposed setup rules, all cases in Table 2 were simulated using the configuration in Table 8. The comparison between the measured and simulated vertical profiles in dimensionless forms of all cases presented high agreement, as shown in Fig. 10c and 10d, except for TC01-TC05 and TC09, which showed unstable and sensitive profiles with lower inflow discharge.
In addition, the five variables including h p , U tc , h tc , C tc and Re tc calculated by Eqs (12-15) of the model results are compared with the measured data of Lee and Yu (1997) in a Taylor diagram as shown in Fig. 11. The principle of the Taylor   Table 8 Optimal numerical aspects for simulation of turbidity currents

Numerical aspects Optimum
Horizontal discretization 2 cm × 0.5 cm (H02) Vertical discretization LSC 2 with more than 60 layers in maximum depth region Vertical TCM TPM, k − ε, and k − ω models Tracer advection scheme TVD 2 or WENO with critical depth ratio lower than 1 diagram (Hofmann et al., 2008;Taylor, 2001) and the model results can be seen in the supplemental materials. The correlation coefficients of the five variables were all larger than 0.9, the normalized centred pattern RMSE was between 0.2 and 0.4, and the normalized standard deviations were between 0.9 and 1.1. The red points of the measured variables were all close to the black point of the observation, which illustrated that the proposed model set-up rules performed well to simulate the turbidity currents under different inflow conditions.

Plunging and transport process of turbidity currents
Simulations using the 3D model with the proposed set-up can accurately reproduce the phenomena of density currents, such as the plunging process and migration of the plunge points in 28 R. Zhang et al. Journal of Hydraulic Research Vol. 61, No. 1 (2023)  0.1 m s -1 , and the C m was 8 kg m -3 without the denser layer structure. The points of the maximum velocity and sediment concentration indicated the structure of the "nose" in the turbidity current head with intense turbulent mixing and large kinetic energy. With the movement of the turbidity current downstream, the point of the maximum velocity moved upward quickly until it reached the converged height of 0.08 m, while the maximum sediment concentration changed from a single point to a converged denser layer with a thickness of 0.1 m. It can be concluded that the vertical profiles of the turbidity current head are distinctly different from the profiles of the body. The head of the turbidity current was in an unstable state and only held for a short time, but the turbidity current body could transport downstream in a similar shape for a long distance.

Conclusions
This study provides comprehensive sensitivity tests for the optimal model set-up of the 3D hydrodynamic numerical model SCHISM to accurately reproduce the evolution of turbidity currents in a sloping flume. The main conclusions drawn from the study are as follows.
(1) Horizontal resolutions parallel and perpendicular to the flow direction affected the turbidity current simulations differently. Larger y had limited effects on the plunge locations, but obviously accelerated the decay process of the sediment concentration along the flume. Smaller x resulted in more converged plunge locations and longi-tudinal profiles of velocity and sediment concentration. In addition, the value of x/ y between 4 and 5 also improved the model results.
(2) The SZ vertical grid was not sensitive to the plunge locations, but caused unacceptable vertical profiles of the turbidity current and obtained fluctuant longitudinal profiles of sediment concentration in the shallow region. More layers of the LSC 2 grid could ensure converged plunge locations, stable vertical profiles of velocity and sediment concentration, and smooth longitudinal profiles in the shallow region.
(3) The turbulence closure models, including the k-ε, k-kl, MY25, k-ω, UB and TPM models, performed distinctly differently from each other. The superiority of the kε, k-ω and TPM models in the stable plunge points, reasonable vertical structures and longitudinal profiles of velocity and sediment concentration was observed, especially for the TPM model. (4) Different transport schemes of the TVD 2 and WENO methods could present consistent results of the plunge locations and vertical profiles of velocity. When the critical depth ratio between the upwind and TVD 2 or WENO schemes was lower than 1, the vertical profiles of sediment concentration were more accurate. (5) The proposed optimal model set-up could successfully reproduce the turbidity currents in a sloping flume with high accuracy and efficiency. The migration of the plunge points and the different vertical structures in the head and body of the turbidity currents could be described by the simulated turbidity currents in detail. The rules can provide helpful advice for other applications to set up the appropriate turbidity current simulation in complex water environments.

Acknowledgements
The 3D model SCHISM used in the study is distributed via github.com, and the source code is available from https://github. com/schism-dev. The numerical tasks were performed on the Explore 100 cluster system of Tsinghua National Laboratory for Information Science and Technology.

Disclosure statement
No potential conflict of interest was reported by the author(s).