INFLUENCE OF TURBULENCE CLOSURE ON ESTUARINE SEDIMENT DYNAMICS AND MORPHODYNAMICS

Turbulence significantly impacts hydrodynamics, mixing and sediment dynamics in coastal environments. We employ a three-dimensional model, the Proudman Oceanographic Laboratory Coastal Ocean Modeling System (POLCOMS), to investigate the effects of implementing various turbulence closure schemes on sediment dynamics and morphodynamics. This model is applied to an idealized estuary, which is represented by a straight rectangular basin. A simple tidal flow is forced at one end and a constant river flow is imposed at the other. Most of the turbulence closure schemes employed are implemented via coupling to the General Ocean Turbulence Model (GOTM). Their effects are also compared to the impact of different erosion parameterizations on the numerical results and observed for different sediment properties.


INTRODUCTION
Sediment transport in estuaries plays a crucial role for coastal zones around the world.Under normal flow conditions, estuaries tend to import sediment (e.g., Meade, 1969) leading to the widespread need for dredging.The estuary turbidity maximum (ETM), a zone of elevated suspended sediment near the limit of salt intrusion, also has an important ecological role.A number of different models are now able to describe such coastal phenomena.However, not all processes can be fully resolved and some need to be modeled instead.There are usually several possible approaches and it becomes necessary to assess their applicability and skill.
A good description of turbulence is fundamental to models describing sediment transport in coastal environments.For example, Warner et al. (2005) showed that different closures can lead to significant differences in suspended sediment concentration.Puleo et al. (2004) showed that the bed shear stress is also impacted by the turbulence closure employed, which would then lead to discrepancies in sediment transport.Given these previous studies, it can be expected that the effect of turbulence should also extend to the morphological evolution that results from convergence or divergence of sediment fluxes.
However, turbulence is not the only modeled process which may impact numerical simulations.In particular, the near-bed boundary conditions for sediment transport have been recognized as an important potential source of error in morphological models.Given the multiplicity of approaches and the existence of user-defined parameters that need to be identified empirically, numerical results may be subject to significant discrepancies.An important issue then concerns how the impacts of different modeled processes compare.
We aim here to investigate the links between implementation of turbulence closure and the resulting sediment transport and morphological evolution for the case of an idealized estuary.These effects are compared to those arising from the use of different approaches for the near-bed sediment boundary condition.To that end, we employ the Proudman Oceanographic Laboratory Coastal Ocean Modeling System (POLCOMS), which we apply to an idealized rectangular estuary.We first provide a succinct description of the coupled hydrodynamics, turbulence and sediment transport model.We then summarize the present numerical setup.The impact of turbulence closures and erosion parameterization is assessed and compared using tidally-averaged cross-sections of salinity and suspended sediment concentration, and the bed level change calculated during the numerical simulations.

MODEL DESCRIPTION
The Proudman Oceanographic Laboratory Coastal Ocean Modeling System (Holt and James, 2001) consists of a three-dimensional hydrodynamic model, which can be coupled to a sediment transport and morphological evolution model (Amoudry et al. 2009), and to the General Ocean Turbulence Model (GOTM) (Holt and Umlauf 2008).

Hydrodynamic model
Hydrodynamics in the model are governed by three-dimensional, incompressible, hydrostatic, Boussinesq fluid equations of motion, which are summarized in Souza and James (1996) and Holt and James (2001).In the present study, the model employs Cartesian horizontal coordinates and the vertical coordinate is taken to be ( ) ( ) , where z is the Cartesian vertical coordinate, h the reference water depth and ζ the water elevation with respect to the reference level.Since we focus on a narrow and shallow estuary, the Kelvin number will remain large, the Ekman number will remain low, and we neglect the Coriolis term.Flow velocities are divided between depth-varying and depthindependent components, and a time-splitting method is used to solve the resulting barotropic and baroclinic equations of motion.
Turbulent stresses and fluxes are modeled following turbulent viscosity and turbulent gradient diffusion hypotheses.As such the model requires the implementation of a turbulence closure scheme, which we discuss in the following section.
Surface stresses are neglected in the present implementation and the bottom stress is calculated based on the near-bed velocity following a logarithmic profile: where F B and G B are the two components of the bed shear stress, and u B and v B are the two components of the near-bed velocity (velocity of the bottom vertical grid) defined at an elevation δ 0 above the bed.κ = 0.41 is the von Karman constant and z 0 the bed roughness, which is taken to be a user-defined constant here.

Turbulence models
Different closure schemes for the eddy viscosity and turbulent diffusivities are implemented either directly in the hydrodynamic model or via GOTM (Umlauf et al., 2005).While these models differ, they all express the eddy viscosity ν t and scalar turbulent diffusivities where the dimensionless quantities μ c and a c μ are referred to as stability functions.All models implemented relate the velocity scale to the turbulent kinetic energy k, which is calculated by solving its balance equation.They mainly differ on the method used to obtain the length scale and on the method used to compute the stability functions.
The length scale can be obtained via an algebraic expression (one-equation model) or via solving a second balance equation for an appropriate quantity (two-equation model).Several quantities can be used for this second equation and here we test the use of the turbulence dissipation rate ε (e.g., Rodi 1987) and turbulence frequency ω (e.g., Umlauf et al. 2003).
Several approaches can also be used to determine the stability functions.The simplest is to set them as constant values.They can also be taken to be empirical functions, usually of the Richardson number, for which several formulations are available.Finally, they can also be derived from algebraic Reynolds stress models under some simplifying assumptions.
We choose to investigate three length scale closure schemes and five methods to calculate the stability functions.All approaches determine the turbulent kinetic energy by solving a balance equation of the type: where the right hand side terms are respectively the diffusive term, the shear production P, the buoyancy production G, and the turbulence dissipation.The diffusive term reads where with k σ the Prandtl number for k.
The length scale can be calculated from an algebraic expression following the ocean turbulence closure implementation of Holt and James (2001): in which σ is negative.The other two approaches solve a second balance equation for ε or ω and then relate the length scale to these quantities: Both models are fully closed by providing values for all the model constants and for the k and ε or ω Prandtl numbers.All are set as the typical default values.The five stability function approaches include three empirical formulations (Munk and Anderson 1948;Eifler and Schrimp 1992;Schumann and Gerz, 1995) and the expressions derived from the two versions of the second-order model of Canuto et al. (2001).In all cases, we also use the length scale limiter of Galperin et al. (1988).

Sediment transport model
The sediment transport model uses the suspended sediment and the morphodynamic modules of the model described in Amoudry et al. (2009).Even though an unlimited number of different sediment classes can be defined, only one is used here and it is described by values for the sediment particle density ρ s and the sediment particle diameter d.The sediment bed is represented by a constant userdetermined number of layers, each of which is characterized by thickness and porosity.
The suspended sediment concentration follows an advection-diffusion equation, which is an expression of suspended sediment mass conservation where the left hand side term accounts for the time rate of change and the concentration change due to advective terms.W s is the sediment settling velocity and is determined following the van Rijn (1993) formula.S c is a source or sink term and D(c) represents the diffusive term.At the free surface, zero sediment flux is imposed.At the sea bed, the diffusive sediment flux is set to vanish, while erosion and deposition are accounted for as source and sink for the bottom grid.
The effect of the suspended sediment on the fluid density is included as where ρ is the density of the fluid-sediment mixture and ρ w is the density of clear water, which is a function of salinity, temperature and pressure.Deposition is due to gravitational settling and occurs at all bed shear stress values, while different erosion parameterizations are implemented.Two approaches are considered in the present study.The first relates the erosion flux directly to an excess bed shear stress via a linear dependence: where E 0 is a user-defined erodibility constant, φ is the top bed layer porosity, τ b the bed shear stress magnitude and τ ce a user-defined critical erosion stress.The second approach relates the erosion flux to a reference concentration value, which we choose to be given by the van Rijn (2007) formula: where θ is the Shields parameter (non-dimensional bed shear stress) and θ cr is the critical bed shear stress.The reference level is chosen as z ref = 20 d with a minimum value of 1 cm, and the nondimensional particle diameter D * is given by For both approaches, the bed erosion is limited by the amount of sediment available in the top bed layer.
The bed morphology follows a mass conservation principle and the location of the sediment bed is updated at each barotropic time step.At each step, a bed level change is calculated from conservation of sediment mass between bed and transported material.Here, only the vertical exchange (deposition and erosion) is considered and horizontal fluxes due to bed load transport are neglected.Such a bed level change can be scaled up by a morphological factor morph F that allows computations for accelerated bed changes (Lesser et al. 2004).Morphological changes are then taken into account in the hydrodynamic model while solving the depth-averaged governing equations.In addition, the bottom boundary condition of the vertical velocity is set as the rate of change of the sea floor.

NUMERICAL SETUP
The ability of this sediment model to accurately reproduce suspended sediment concentrations (SSC) and morphological evolution under steady currents has been shown in Amoudry et al. (2009).
Here, we apply it here to the case of a rectangular idealized estuary of length 100 km, width 3 km, and initial uniform depth 15 m.The numerical domain is such that north and south boundaries are land boundaries.A simple tidal flow is prescribed at the open west boundary.Finally, the east boundary is a weir at which a constant river flow is imposed.This domain is discretized with 200 m horizontal resolution in both directions and twenty vertical levels.A higher vertical resolution of 40 levels was tested but was not found to significantly modify the numerical solutions.
At the west boundary, a flux-radiation boundary condition (Flather, 1981) is implemented, and the external elevations and depth-averaged currents are prescribed analytically as where ζ 0 and U 0 are the elevation and current amplitudes, φ a phase lag/lead of the current respect to the elevation and T the tidal period.Such a condition generates tidal waves propagating in the numerical domain and lets waves leave the domain with minimal reflections.The tide is forced with ζ 0 = 1.3 m, U 0 = 1.5 m/s, φ = -0.25 and T = 12 hours.The salinity at the west boundary is relaxed towards a constant value of 30 psu over five grid points, and the temperature is kept constant in the entire domain.At the east boundary, a constant clear freshwater river flow of 150 m 3 /s is prescribed via an increase of the sea surface elevation according to the volume flux and a corresponding adjustment to the salinity and sediment concentration throughout the water column.
Numerical simulations are started from rest with a level free surface and a constant 30 psu salinity in the estuary.Hydrodynamics only are first simulated for 60 tidal cycles, suspended sediment is then simulated for a further 60 cycles, and finally the bed evolution is calculated over 10 cycles with a morphological factor of 10, corresponding to 100 cycles.
We conduct a series of numerical experiments to investigate the impacts of both the erosion parameterization and the turbulence closure on sediment transport for 160 µm diameter sand (particle density of 2650 kg/m 3 ).To that end, eight cases were considered (table 1).Two cases use a k-ω model with the same stability function determined from the second-order model of Canuto et al. ( 2001) but each with a different erosion formulation approach.In all cases, the critical stress is taken to correspond to the Shields curve value (τ ce = 0.17 N/m 2 ), and when equation 10 is employed, the erodibility constant (E 0 = 0.082 kg/m 2 /s) has been determined from model-data comparison of suspended sediment concentration profiles (Amoudry et al., 2009).Three cases investigate the influence of the turbulence length scale closure by using either a k-ω model, a k-ε model, or a oneequation k model.Finally, five cases investigate the influence of the stability function by implementing the different formulations introduced previously.

RESULTS
Numerical results for the tidally-averaged suspended sediment concentration and the tidallyaveraged isohalines are presented in figures 1 and 2. The total bed changes calculated by the morphological model are presented in figures 2 and 4. The case using the reference concentration (case 2) is included in all four figures for comparative purposes.In figures 1 and 2, all cases but one clearly exhibit the presence of a maximum in suspended sediment concentration around 40 km from the mouth.This maximum occurs near the salt intrusion length, which is consistent with observations of the ETM in partially mixed estuaries (e.g., Postma, 1967;Schubel, 1968;Allen et al., 1980;Sanford et al., 2001).While the present model application can not be compared to experimental data, this still highlights the ability of the model to reproduce known observed physical phenomenon at least qualitatively.Figures 3 and 4 display significant bed evolution for this idealized estuary scenario.In all cases, erosion occurred near the mouth and no change is observed at the upstream end.Such a behavior at the upstream end is expected and consistent with numerical suspended sediment concentrations (figures 1 and 2), as bed shear stresses are not large enough to suspend particles and result in transport.In the rest of the estuary, two main areas of deposition can be seen.One is located around 40 km from the mouth and relates to both the ETM and the limit of salt intrusion observed in figures 1 and 2. A possible simple explanation is that the convergence of the estuarine circulation at the limit of salt intrusion results in a convergence of sediment fluxes, and thus deposition.The second deposition area occurs 15 to 20 km from the mouth.In this case, a localized net deposition is explained by smaller near-bed velocities during both ebb and flood, which are related to a stronger horizontal salinity gradient.

Influence of the erosion parameterization
The influence of the erosion parameterization is easily assessed by comparing the top two panels of figures 1 to 4. Overall, there is relatively little change.In particular, the qualitative observations about ETM, deposition areas, and their locations are not modified.However, the near bed concentrations using the van Rijn (2007) reference concentration are significantly smaller than those using the linear erosion formula for small to moderate stresses (around 70 to 80 km in figure 1).This is not the case for higher stresses (closer to the mouth in figure 1).Such discrepancy is probably due to the different power relationships in equations 9 and 10.

Influence of the turbulence length scale closure
As displayed in figures 1 and 3, different turbulence length scale closures can lead to dramatic discrepancies in the numerical results.For example, no clear ETM can be distinguished in the bottom panel of figure 1 when using the one-equation k model.Both the isohalines and the salt intrusion length are also drastically different.The absence of ETM indicates an important failure of the model to reproduce a known physical phenomenon.While this could be caused by both the lack of a second balance equation and a different stability function, the impact of stability functions (see next section) is not sufficient to justify such discrepancies, and we therefore attribute the observed failure to the use of an algebraic expression for the length scale.Comparatively, both two-equation models (k-ω and k-ε) lead to very similar salinities and suspended sediment concentrations, but slight changes in position the salt intrusion and the ETM can be observed.The amount of suspended sediment also differs slightly between the two closures mainly due to differences in calculated bed shear stress.Significant changes can also be observed for the bed response (figure 3).Once again, the oneequation k model displays very different behavior.Unlike the two-equation closures, there is no clear large deposition area 20 and 40 km from the mouth.The morphological changes also extend further upstream to 70 km, which exceeds the increase in salt intrusion length of 60 km observed in figure 1.Both two-equation models show similar qualitative behavior over the estuary, in that deposition occurs at the same locations.However, there is a significant difference (up to a factor of 2) in the amount of sediment deposited, which is much larger than the impact due to the erosion parameterization.

Influence of the stability functions
Implementing different stability functions can also have a significant impact on the sediment dynamics and morphodynamics (figures 2 and 4).While the suspended sediment concentration and isohalines are qualitatively somewhat similar, there can be substantial discrepancies in the salt intrusion length.This, in turn, results in a shift of the ETM (figure 2) and of one deposition area (figure 4).The amount of suspended sediment also differs depending on the stability function employed (figure 2), and the empirical stability function seems to induce more suspended sediment in the water column (especially between 20 and 30 km from the mouth).While one deposition area is significantly shifted due to differences in salt intrusion length, the location of the other is not impacted and stays around 20 km from the mouth.The amplitudes of observed numerical bed evolution can also vary substantially (again by a factor of 2) and are more important than the impact of the erosion parameterization.
However, the effects of the stability functions mainly occur between implementations relying on a second-order model and empirical expressions.Comparatively, there are only moderate differences between the results obtained with the two implementations derived from second-order models, or between the results obtained with the three empirical formulations.

Persistence of observations for different conditions
Similar numerical simulations to those summarized in table 1 (with the exception of case 2) have been performed for other sediment characteristics and for other tidal flow inputs.As an example, figures 5 and 6 present the tide-averaged suspended sediment concentration, salinity contours and the bed evolution numerically observed for different turbulence length scale closures for smaller particles.The particle diameter for these cases is d = 40 μm and the linear erosion is employed (i.e., equation 10) with E 0 = 10 -4 kg/m 2 /s and τ ce = 0.20 N/m 2 .Similar characteristics can be observed in figures 1 and 5.There is an ETM for the two-equation models at the salt intrusion length, which is the same in both cases.This is not surprising given that the salt intrusion length is determined by the hydrodynamics.There is, however, significantly less sediment suspended for the fine particles upstream of the salt intrusion, which is due to a relatively much higher critical erosion stress.For the bed evolution, only the downstream deposition area is observed and the turbulence length scale closure has an effect similar to that observed in figure 3. The deposition area related to the ETM is not observed in figure 6 because the flow convergence does not result in a sediment flux convergence with such little suspension upstream of the salt intrusion length.
In summary, similar observations can be made concerning the effect of turbulence closures for smaller particles (non-cohesive silt of diameter 40 μm here).Even though they are not shown here, similar observations can also be made for different tidal inputs.The one-equation k model leads to numerical results that show very significant differences to all other cases and that are not able to predict the presence of an estuary turbidity maximum.Comparatively, the two two-equation models (kω and k-ε) give qualitatively very similar results, but with significant quantitative differences in deposition.The effect of the stability functions can be summarized as occurring mainly between implementations derived from second-order models and empirical formulations, and it focuses on the salt intrusion length, and quantities that depend on it.

DISCUSSION
The present numerical results confirm that the effects of turbulence closure significantly impact estuarine dynamics and suspended sediment.They also highlight the impact on potential morphological change.In particular, in the present case, the effects of turbulence are more important than that of a change in erosion parameterization.However, precise information on which turbulence closure approach ought to be preferred is hindered by the lack of data for the present study.
In the case of the one-equation k model, this problem is actually relatively simple.The results display very large discrepancies with all other cases and are unable to reproduce key known physical processes.Such obvious shortcomings are also compounded by important theoretical flaws.Results are inherently dependent on the algebraic expression for the length scale, and specifications are flow specific.In that sense, the model is incomplete, and thus significant errors can be found for other flows.
A clear superiority is less evident between the two two-equation models, as we did not observe a big difference between them.Even though they are mathematically different, the k-ε and k-ω models are physically identical with the exception of the diffusive terms.The models are actually mathematically identical for homogeneous turbulence (Pope, 2000).Given such similarity, the qualitative agreement between the numerical results should not be surprising, and the quantitative differences can probably be attributed to the different diffusive terms.It thus appears that the k-ω model could be a viable alternative to the k-ε model.It is also worth noting that the k-ε model is known to fail in some situations for which the k-ω model does not (see, for example, Pope, 2000), and the latter approach could thus provide a more generally robust solution.
Similarly, the superiority of a given stability function is not straightforward.However, the stability functions derived from second-order models are actually a consequence of the model and can be mathematically derived.They thus offer an important theoretical advantage compared to empirical functions of the Richardson number.

CONCLUSION
We used a coupled turbulence, sediment transport, and three-dimensional hydrodynamic model to investigate the impact of several turbulence closures on the numerical results in an idealized estuarine scenario.This model has been able to appropriately reproduce the well-known physical phenomenon that is the estuarine turbidity maximum.
The effect of turbulence closures has been found to be significant.In particular, it is found here to be more important than the change in erosion parameterization implemented.While most turbulence closures lead to qualitatively similar results, the one-equation k model produces results which are qualitatively erroneous.As such, we advise against its use in estuarine applications.While some significant quantitative differences were observed for the remaining two-equation closures, the absence of data makes detailed assessment of the quality of a given turbulence closure difficult.

Figure 1 .
Figure 1.Sediment suspended concentration (SSC) and salinity contours for the two erosion parameterizations and three different turbulence length scale closures.From top to bottom are cases 2, 1, 3 and 4. The salinity contours are plotted for 1 psu, 2 psu, and every 2 psu after from right to left, the bold contours correspond to 10 psu and 20 psu.

Figure 2 .
Figure 2. Sediment suspended concentration (SSC) and salinity contours for the two erosion parameterizations and the five different stability functions.From top to bottom: cases 2, 1, 5, 6, 7 and 8.The salinity contours are plotted for 1 psu, 2 psu, and every 2 psu after from right to left, the bold contours correspond to 10 psu and 20 psu.

Figure 3 .
Figure 3. Bed level change (m) for different erosion parameterizations and turbulence length scale closures.From top to bottom: cases 2, 1, 3 and 4.

Figure 5 .
Figure 5. Suspended sediment concentration (SSC, g/m3) and salinity contours for 40 μm diameter particles under case 1 (top), case 3 (middle) and case 4 (bottom).The salinity contours are plotted for 1 psu, 2 psu, and every 2 psu after from right to left, the bold contours correspond to 10 psu and 20 psu.