Investigating Suspended-Sediment Transport in a Shallow Lake Using a Three-Dimensional Model

A three-dimensional, unstructured grid, hydrodynamic and suspended-sediment transport model (i.e., SELFE-SED) was developed to simulate temporal and spatial variations of suspended sediment and was applied to the subtropical subalpine Tsuei-Feng Lake (TFL) of Taiwan. The model was validated with measured water level and suspended-sediment concentration in 2009, 2010, and 2011. The overall model simulation results are in quantitative agreement with the observational data. The validated model was then applied to explore the most important parameter that affects the suspended-sediment concentration and to investigate the effect of wind stress on the mean current and suspended-sediment distribution in this shallow lake. Modeling results of sensitivity analysis reveal that the settling velocity is a crucial parameter and erosion rate is less important in the suspended-sediment transport model. Remarkable lake circulation was found based on the strength of wind speed and wind direction. Strong wind would result in higher mean current in the top layer and suspended-sediment distribution in the top and bottom layers. This study demonstrated that the wind stress played a s ignificant influence on mean circulation and suspended-sediment transport in a shallow lake.


Introduction
Sediment has been identified as one of the important nonpoint source pollutants [1].Sediment-water interactions in shallow lakes may become enormously important with time, since bed sediments constituent the ultimate repository for nutrients and contaminants, and fine particles and nutrients in the sediment bed may repeatedly be recycled [2,3].Nutrient-rich sediment layers in the lakes participate in chemical and biological processes and exchange mass with in the water column.The s ediment layers affect the nutrient cycles via the diffusive fluxes of nutrients from sediment and via sediment resuspension [4][5][6].The resuspension of toxic bottom sediments is an ecological concern, since toxic contaminants can accumulate in the biota and food web of a lake [7].Suspended sediment can also reduce the light intensity under water and affect the growth of phytoplankton.
The dynamic processes involved in suspended-sediment transport such as flocculation, deposition, and erosion (resuspension) have been studied by many scientists [8][9][10][11][12][13][14][15][16].A number of laboratory studies have improved our knowledge of t he resuspension behavior of fine-grained material.However, transferring the results of laboratory experiments to field studies still has been problematic because of the complexities of real sediments and of natural flows.
Sediment transport models have b een broadly developed and applied to investigate sediment dynamics in lakes including a one-dimensional model [17][18][19], a vertical two-dimensional model [20,21], a horizontal two-dimensional model [22][23][24][25][26], and quasi-three-dimensional [27] and three-dimensional models [28][29][30][31][32].The erosion and resuspension of bottom sediment from the active layers is one of the major sources of the total suspended sediment in the water column [33].To accurately simulate the suspended-sediment concentration in lakes, incorporating resuspension and deposition into the suspended-sediment transport model is essential.However, three-dimensional models are particularly appropriate in cases such as regions with complex bathymetries, the vertical gradients in water temperature and suspended-sediment concentration, and density currents in the lakes.
Recent advances in computational fluid dynamics and sediment transport model formulations enable the computation of flows and sediment fluxes in three-dimensions.However, spatial and temporal scales, and model processes must be carefully selected for a particular study to ensure that the model will solve the problem at hand and that computations can be carried out in the practical sense.In the present study, a three-dimensional, unstructured-grid hydrodynamic and suspended-sediment transport models (i.e., SELFE-SED) were developed and applied to simulate the hydrodynamics and sediment dynamics in the Tsuei-Feng Lake.The model simulation was validated against profiles of water level, and time-variation and spatial-variation suspended solids from 2009 to 2011.The validated model then was applied to explore the most sensitive parameter affecting the suspended-sediment concentration and the influence of wind stress on mean currents and sediment dynamics.

Description of Study Area
Tsuei-Feng Lake (TFL) is located between the Taiping Mountain and Dayuan Mountain in the north-central region of Taiwan (Figure 1a).The indicated altitude is 1,900 m above sea level.The area of the lake is approximately 20 ha.It is one of the largest mountain lakes in Taiwan.It is also the most popular focus in the Taipingshan National Forest Recreation Area.
The regulation of the water flow into TFL primar ily results in extreme water level fluctuation, which severely affects the activity and succession of fauna and flora in the area.The average annual temperature is approximately 11.4°C, and the annual precipitation is more than 3,662 mm.TFL is subject to several typhoons in the summer and autumn each y ear, during which more than 1,300 mm of precipitation may fall on the lake.The water column is stratified from early April to October.The water column is usually completely mixed in the winter and is associated with intensive rainfall during the typhoon season.The water quality variables of pH, dissolved oxygen, total nitrogen, and chlorophyll a are in the range of 5.4~6.5, 5.9~7.7 mg/L, 0.5~1.4mg/L, and 0.4~6 g/L, respectively [34].The spatial distributions of suspended-sediment concentration on June 6 and December 11, 2010 are shown in Figure 2.These indicate that the suspended-sediment concentration is higher on December 11 (during the autumn season), while it is lower on June 6 (during early summer).The suspended-sediment concentration in the lake exhibits temporal and spatial variations.According to a sieve analysis, the bottom sediment type can be described as being between fine sand and silt [34].

Hydrodynamic Model
Multi-scale lake circulation modeling has become a mature field.The models used typically describe the free-surface elevation, velocities, and water temperature of various water bodies by solving the Reynolds-averaged Navier-Stokes equations which represent conservation of mass, momentum, and water temperature subject to the hydrostatic and Boussinesq approximations.The equation of state describes water density as a function of water temperature.In this study, the three-dimensional, semi-implicit Eulerian-Lagrangian finite-element model [35] was implemented to simulate the hydrodynamics and hydrothermal characteristics in the Tsuei-Feng Lake.
The SELFE m odel solves the governing equations using a semi-implicit finite-element/volume scheme with a key step of decoupling continuity and momentum equations via the bottom boundary.In SELFE, the horizontal domain is discretized with unstructured triangular grids.Along the vertical direction, a hybrid coordinate system consisting of the S-coordinate and Z-coordinate in the upper and deeper parts of the water column, respectively, can effectively prevent the so-called hydrostatic inconsistency issue.The model handles the advection terms in the momentum equation with an Eulerian-Lagrangian method, allowing larger time steps without compromising computational stability and accuracy.For the advection terms in the transport equations (i.e., water temperature), the Eulerian-Lagrangian scheme is applied.When the Eulerian-Lagrangian scheme is used, the transport equations can be efficiently solved by a finite-element method.To calculate turbulent mixing processes, SELFE uses the generic length scale (GLS) turbulence closure of Umlauf and Burchard [36], which has the advantage of encompassing most of the 2.5-equation closure

Suspended-Sediment Transport Model
The suspended-sediment transport model is treated in two distinct models: the advection-diffusion components are based on the SELFE transport formulation, taking advantage of its numerical characteristics, such as the option to use different methods to solve the transport equation, while the remaining terms are computed based on the sediment transport formulation of a three-dimensional hydrodynamic-eutrophication model [37].
The suspended-sediment transport model was incorporated using the hydrodynamic model as a passive tracer under the same spatial and temporal resolutions.The governing equation for the total suspended-sediment concentration is as follows: where u , v and w are the velocity components in the x , y and z directions provided by SELFE, respectively, s w is the settling velocity of the suspended sediment,  is the vertical eddy diffusivity and s F represents horizontal diffusion.
At the water surface, there is no net transport across the free surface, and, therefore, diffusion flux always counterbalances the settling flux and the boundary condition where  is the free-surface elevation.
At the sediment bed, the net sediment flux is equal to the summation of the sediment erosion flux and the sediment deposition flux.
The erosion term  

E kgm s
 is expressed according to the Partheniades [8] formulation, using the excess-bottom shear-stress concept:  stratification and gravitationally driven flow, the effects of suspended-sediment concentration on density have to be considered.Water density including the suspended sediment was computed as follows: where w  is the water density, which is determined by the equation of state, and s  is the suspended sediment density.

Model Implementation
In this study, the bottom topography data in the TFL measured in 2008 were obtained from the Academia Sinica, Taiwan.The greatest depth of study area is 4.3 m where is close to the buoy station (Figure 1a).The unstructured model meshes for TFL consist of 11617 meshes/6058 nodes in the horizontal direction (Figure 1b).High-resolution unstructured meshes, which include a mesh size of approximately 5 m, were used for TFL.Model simulations were conducted using ten S-levels in the vertical direction.For this model mesh, a time step ( t  = 20 seconds) was used in simulations to guarantee the numerical stability.

Indices of Model Performance
To evaluate the performance of the hydrodynamic and suspended-sediment concentration, two criteria were adopted to compare the predicted results and the observational data, which are the mean absolute error (MAE) and root mean square error (RMSE).MAE is the average of the absolute values of differences between the observed data and simulated values.RMSE basically specifies the overall difference in the sum of squares normalized to the number of observations.RMSE is similar to a standard error of the mean for the uncertainty of the model.These criteria are defined by the following equations: where N is the total number of data points and

Model Validation
Simulation models are increasingly being used to solve problems and to aid in decision-making.To make sure of the model accuracy for further practical applications, the observational data is used to validate the model and to ascerta in its capability for predicting water level and suspended-sediment concentration [38].
The observational data collected from November 2009 to July 2011 was used for model validation.To warm up the model, the spin-up time is specified as 15 days.Th e model was therefore was run from November 1, 2009 to July 31, 2011.The initial conditions for water level and suspended-sediment concentration are specified as 3.0 m and 13.0 mg/L, respectively, to start running the hydrodynamic and suspended-sediment transport model.

Water Level
The model validation of the water level was executed with daily discharge at the inflow and outflow in the lake in addition to precipitation. Figure 3a

Suspended-Sediment Modeling
Several parameters in the suspended-sediment transport model including the settling velocity ( s w ), critical stress for erosion ( ec  ), deposition ( dc  ), and empirical erosion rate ( M ) are difficult to be determined in the lakes because they have a wide range of value to be adopted.The consolidation of the sediment bed affected by interactions with physical, biological, chemical, and the composition of sediment, is an important factor influencing these parameters [39][40][41][42].In this study, the suspended-sediment transport model is implemented to confirm these parameters.
In order to validate the suspended-sediment concentration in the TFL, the monthly measured data collected from November 2009 to July 2011 was compared with the simulated results.The water sample was taken from different water depths to measure the concentration of suspended sediment.Concentrations of suspended sediment were determined using the drying method after filtering samples through GF/F filters [43].The comparison of the simulated suspended-sediment concentration and the measured concentration taken as a vertical average at the buoy station is shown in Figure 4.It reveals that simulated results fairly match the measured suspended-sediment concentration.The MAE and RMSE values between the computed and measured suspended-sediment concentrations are 1.73 mg/L and 2.23 mg/L, respectively.
Figure 5 presents the vertical profiles of simulated and measured suspended -sediment concentrations at the buoy station in 2010.It indicates that the model simulation results match the measured suspended-sediment concentrations in vertical profiles.6 illustrates the comparison of the spatial distribution of the vertically averaged suspended-sediment concentration between simulated and observed results on April 4, July 29, November 6, 2010 and January 28, 2011.It can be seen that the simulated and observed suspended-sediment concentrations in their spatial distribution are quite similar.Through the model validation, the constant values of 1.0x10 -3 m/s for settling velocity ( s w ), 5.0x10 -3 N/m 2 for ec  , 5.0x10 -4 N/m 2 for dc  , and 5.0x10 -6 kg/m 2 .S for M were chosen for the suspended-sediment transport model.

Model Sensitivity analysis
Sensitivity analysis is the study of how the uncertainty in the output of a mathematical model can be apportioned to different sources of uncertainty in its inputs .Three-dimensional suspended-sediment transport models have been applied to examine the sensitivity of their models to alter the physical forcing [44][45][46].In this study, model sensitivity analysis was implemented to explore the effects of settling velocity ( s w ), the critical stress for erosion ( ec  ), deposition ( dc  ), and the empirical erosion rate ( M ) on the suspended-sediment concentration and to determine the most important factor that affects the suspended-sediment concentration in the TFL.
The original bases depend on the simulation of model validation from April 1, 2010 to April 30, 2010.The effects of these parameters on the suspended-sediment concentration w ere investigated with two alternative cases; one involves the value used in model validation plus 50% and the other involves the value minus 50%. Figure 7 presents the model sensitivity results for the four parameters during April 2010.It indicates that the settling velocity would be the most important parameter affecting the suspended-sediment concentration in the TFL.Table 2 lists the model sensitivity results.It indicates that an increase in settling velocity results in a decrease in the suspended-sediment concentration.The maximum rates for decreasing and increasing suspended-sediment concentration are 43.85% and 179.15%, respectively.The maximum rate is the suspended-sediment concentration for the sensitivity run shown in Figure 7.
The most and least important factors to affect the suspended-sediment concentration are the settling velocity and the erosion rate, respectively, shown in Table 2 based on model sensitivity runs.Lee et al. [29] and Hawley et al. [20] used a one-dimensional (vertical) resuspended bed model and a two-dimensional (vertical and cross-shore) sediment transport model, respectively, to implement model sensitivity to identify the important resuspension parameters in Lake Michigan.They found that the settling velocity was crucial parameter in controlling the prediction of suspended-sediment concentration.Their simulated results of sensitivity analysis were similar with the current study with a three-dimensional suspended-sediment transport model.However, we further found that the erosion rate is the least important parameter in the suspended-sediment modeling.

Influence of wind stress on current and suspended sediment
The flows, surface waves, and sediment resuspension in the shallow lakes are significantly driven by wind forcing [47][48][49].Podsetchine and Huttula [50] addressed the correlation and time lag of wind direction/speed to resuspended sediment in Lake Karhijärvi, Finland.Jin and Ji [51] applied a three-dimensional model to investigate the effects of wind velocity and fetch and wind persistence on suspended-sediment concentrations in Lake Okeechobee.The mean circulation in the TFL would be one of the important factors to affect the transport and distribution of suspended sediment.To investigate how the suspended sediment could be transported and distributed in the YYL, different wind speeds and directions including northeast wind: 2.23 m/s (baseline), northeast wind: 9.39 m/s, northwest wind: 9.39 m/s, and southwest wind: 9.39 m/s were used to drive the model simulation to yield mean circulation and suspended-sediment distribution.The wind speed and direction of 2.23 m/s from the northeast in November 2010 were selected as baseline.The maximum wind speed is 9.39 m/s during the model validation period and prevailing wind directions are northeast, northwest, and southwest adopted for model simulations to compare with the baseline condition.
The mean currents at the top and bottom layers for the baseline condition and for different wind speeds and directions, respectively, are shown in Figure 8 and Figure 9.The obvious anticyclonic circular gyre appears at the top layer for the baseline (Figure 8a), while the cyclonic circular gyre is seen at the top layer for the wind (Figure 8d).Basically, the direction of the mean current at the top layer follows the wind direction.A strong wind results in high mean current at top layer.The bottom currents are smaller than the surface currents.The mean surface currents are between 0.018 and 33.0 cm/s, while the mean bottom currents are between 0.01 and 31.1 cm/s.
Figures 10 and 11 show the distributions of suspended-sediment concentration at the top and bottom layers, respectively, for the baseline condition and for different wind speeds and directions.Both figures indicate that the distributions of suspended-sediment concentration are obviously different depending on the wind speed and direction.A strong wind produces high suspended-sediment concentrations in the surface and bottom layers.The highest concentration appears near the buoy station, while the lowest concentration appears at the outlet.The mean surface concentrations are between 19.0 and 489.7 mg/L, while the mean bottom concentrations are between 19.4 and 491.9 mg/L for a wind speed of 9.39 m/s.We found that the spatial distribution of the suspended-sediment concentration was subject to wind speed.Figure 12 illustrates the simulated results of the vertical suspended-sediment profile at the buoy station under conditions of different wind speeds and directions.It shows that the suspended-sediment concentrations in the vertical direction are homogeneous under conditions of a northeast wind, northwest wind, and southwest wind at 9.39 m/s, while stratification occurred with a northeast wind at 2.23 m/s.The different wind speeds and directions also produced different patterns in the vertical suspended-sediment concentration.The northwest wind at 9.39 m/s produces the highest suspended-sediment concentration at the buoy station (Figure 12c).It implies that the water quality conditions in the lake may be affected by the wind stress.Jin and Ji [51] documented that the effects of wind-induced current and wind-induced resuspension on suspended-sediment transport were important factors affecting the suspended-sediment concentration in Lake Okeechobee.Cozar et al. [52] presented empirical correlations between turbidity and wind speed based on the field measurements.They yielded a formula and used the wind speed and water depth to calculate the suspended-sediment concentration.It meant that strength of wind speed was high correlation with the suspended-sediment concentration.In the present study, we also found that wind stresses have significantly influence on the mean circulation and suspended-sediment transport and distribution in a shallow lake.

Conclusions
A three-dimensional hydrodynamic and suspended-sediment transport model was developed and applied to the Tsuei-Feng Lake in the north-central region of Taiwan.The model was validated with observed water level and suspended-sediment concentrations in 2009, 2010, and 2011.The predicted results quantitatively agreed with measured data in the lake.
The validated model was applied to explore the crucial factor that most affects the suspended-sediment distribution in the TFL.Model sensitivity analysis indicates that the settling velocity is the most important parameter affect ing the suspended-sediment concentration and the erosion rate is the least important parameter in the suspended-sediment transport modeling.Therefore, this parameter (i.e.settling velocity) should be carefully adjusted when the model is used for the model validation procedure.The model was also used to probe the mean current and suspended-sediment distribution in the TFL.Wind speed and direction are the dominant factors driving the current patterns and affecting the suspended-sediment concentration.The obvious anticyclonic and cyclonic circular gyres appear in the lake based on the wind stress.The distributions of suspended-sediment concentration are clearly different depending on the wind speed and direction.Strong winds result in high mean currents and suspended-sediment concentrations in the top and bottom layers due to the effect of resuspension in a shallow lake.

Figure 1 .
Figure 1.(a) Map of Tsue i-Fe ng Lake (TFL) in the north-ce nte r re gion of Taiwan and (b) horizontal grid of TFL for three -dime nsional hydrodynamic and suspe nde d-se dime nt transport mode l.The grid size is approximate ly 5 m.

Figure 2 .
Figure 2. The spatial distribution of suspe nde d-se dime nt concentration on (a) June 6, 2010 and (b) De ce mber 11, 2010.The white color in the figure re pre sents the dry area (no wate r).
where cd  is the critical shear stress   2 Nm  for deposition, and b C is the near-bed suspended-sediment concentration.For simulating the hydrodynamic properties such as density Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 29 May 2017 doi:10.20944/preprints201705.0203.v1 are the values the ith data point, modeled and observed, respectively.
presents the recorded rainfall which reached 48.8 mm/hr on October 21, 2010 when Typhoon Megi hit Taiwan.The observed water levels in the TFL from November 1, 2009 to July 31, 2011 were compared with the model prediction shown in Figure 3b.The MAE and RMSE values between computed and observed water surface elevation are 0.05 m and 0.15 m, respectively.It indicates that the simulated results mimic the observed water levels.Based on the model validation of water level, the roughness height used in the hydrodynamic model was adopted to be 0.01 m for the model simulation.

Figure 3 .
Figure 3. (a) Recorde d rainfall from Novembe r 1, 2009 to July 31, 2011 and (b) comparison of water le ve l be twe en model pre diction and observation.

Figure 4 .Figure 5 .
Figure 4.The comparison of suspe nde d-se dime nt conce ntrations pre dicte d by the mode l and me asured at the buoy station from Nove mber 2009 to July 2011.

Figure 6 .
Figure 6.Comparison of spatial distributions of vertically ave rage d suspe nde d-se dime nt conce ntration be twee n mode l simulation and observation on April 4, July 29, Nove mber 6, 2010 and January 28, 2011.Note that figures (a) and (b) show the observation and simulation on April 4, 2010.Figures (c) and (d) prese nt the observation and simulation on July 29, 2010.Figures (e ) and (f) are the obse rvation and simulation on Nove mbe r 6, 2010.Figures (g) and (h) prese nt the observation and simulation on January 28, 2011.
maximum values were determined by the formula represented by suspended-sediment concentration for the base runs shown in Figure4and sens C

Fig ure 7 .
Fig ure 7. Mode l se nsitivity for the (a) se ttling ve locity ( s w ), (b) critical stress for de position (

Figure 8 .
Figure 8. Mean circulation at the top layer for diffe re nt wind spee ds and directions (a) northeast wind: 2.23 m/s, (b) northeast wind: 9.39 m/s, (c) northwest wind: 9.39 m/s, and (d) southwest wind: 9.39 m/s.

Table 1
presents the statistical errors between simulated and measured suspended-sediment concentrations for each measur ed Preprints (

Table 1 .
Statistical e rror betwee n simulate d and measure d suspe nde d-se dime nt conce ntrations from Nove mbe r 2009 to July 2011.

Table 2 .
Re sults of mode l se nsitivity run for four parameters .