Formulating Fine to Medium Sand Erosion for Suspended Sediment Transport Models

The capacity of an advection/diffusion model to predict sand transport under varying wave and current conditions is evaluated. The horizontal sand transport rate is computed by vertical integration of the suspended sediment flux. A correction procedure for the near-bed concentration is proposed so that model results are independent of the vertical resolution. The method can thus be implemented in regional models with operational applications. Simulating equilibrium sand transport rates, when erosion and deposition are balanced, requires a new empirical erosion law that involves the non-dimensional excess shear stress and a parameter that depends on the size of the sand grain. Comparison with several datasets and sediment transport formulae demonstrated the model’s capacity to simulate sand transport rates for a large range of current and wave conditions and sand diameters in the range 100–500 μm. Measured transport rates were predicted within a factor two in 67% of cases with current only and in 35% of cases with both waves and current. In comparison with the results obtained by Camenen and Larroudé (2003), who provided the same indicators for several practical transport rate formulations (whose means are respectively 72% and 37%), the proposed approach gives reasonable results. Before fitting a new erosion law to our model, classical erosion rate OPEN ACCESS J. Mar. Sci. Eng. 2015, 3 907 formulations were tested but led to poor comparisons with expected sediment transport rates. We suggest that classical erosion laws should be used with care in advection/diffusion models similar to ours, and that at least a full validation procedure for transport rates involving a range of sand diameters and hydrodynamic conditions should be carried out.


Introduction
Classically, modelling strategies are different for sand and mud transport. This is due to the fact that, whereas mud is advected with fluid velocity, sands are not, and this is all the more true the larger their diameter [1]. The reason for this difference is the much higher settling velocity of sand particles. Mud transport is generally modelled using an advection-diffusion equation. Transport rates are computed as the vertical integration of the product of fluid velocity and the concentration of suspended sediments (e.g., [2]). Because equilibrium is reached very rapidly, for instance in flume tests, sand transport dynamics have been generally formulated as vertically-integrated horizontal transport rates, using so-called parametric models empirically deduced from experimental data (e.g., [3][4][5]).
Many coastal environments, including estuaries and continental shelves are composed of mud and fine sand carried from rivers or from the sea. Mud and sand are likely to mix, and the resulting physical properties of the sediment, as well as its dynamics can be highly modified [6][7][8]. In this context, models accounting for both cohesive and non-cohesive sediment dynamics have aroused a growing interest (e.g., [9][10][11][12][13][14][15]). Using the same mathematical formalism for sand and mud transport is particularly appealing in those models, mostly for the following reasons. Firstly, it can be seen as a prerequisite to avoid the temporal and spatial numerical discontinuity of erosion and transport flux at the frontier between cohesive and non-cohesive behaviours [11][12][13]. Moreover, managing variations of surficial sediment composition requires considering thin layers of sediment, whose content directly depends on simultaneous vertical fluxes of the sandy and the muddy fractions [11]. Finally, erosion of mixed sediments can be implemented in models using flume experiments. Since equilibrium is generally not reached in most cases of a sediment mixture, formulating sand transport through erosion/deposition fluxes is needed to be consistent with erosion experiments.
In recent sediment modelling studies dealing with sand/mud mixture in estuaries or over shelves, it has become usual to only consider a suspension transport mode for both sand and mud, resolving a different advection-dispersion equation for each sediment fraction (e.g., [11][12][13][16][17][18]). The bed load transport is usually not included in those models for fine to medium sands, and the total transport is expected to be reproduced by the suspended load only. This may only be valid for very fine sand, but for fine to medium sand, it can be assumed that the saltation transport and possibly the bed load can be included in the suspension transport within the boundary layer very close to the bed, provided the erosion rate can be calibrated for that purpose (i.e., fitting the erosion rate to get the right vertically-integrated horizontal flux at equilibrium). The validity of such a strategy is not well documented. Only a few validations of the approach, which consists of computing sand transport by solving an advection-diffusion equation with a pick-up function as the source term and deposition flux as the sink term, have been reported. The capacity of this method to reproduce experimental total fluxes (bed load + suspension) has often been tested for a narrow range of sand (200-250 μm) [19,20]. However, the possible use of this approach in a wider range of sand grain sizes requires further validation.
Using an advection-diffusion equation for sand transport is not easy to implement for several reasons. Sand transport takes place very close to the bed and is highly vertically variable. Most of the numerical models used at the shelf scale are not sufficiently vertically refined to resolve processes in the bottom boundary layer. In addition, the expression of the bottom boundary condition for sand is a major difficulty. Several formulations based on experiments have been proposed using pick-up functions or reference concentrations [4,[21][22][23][24][25][26]. The distance above the bottom where such a condition has to be specified remains a sensitive parameter in those formulations. Moreover, at equilibrium, the horizontal solid fluxes that result from such computations not only depend on the pick-up function, but also on numerical parameters (e.g., the vertical resolution) or physical processes implementation (e.g., model turbulence closure or deposition flux formulation).
The aim of this study is to validate an advection/diffusion sand transport model, which could in the future be implemented in a regional sediment transport model, that is applicable with a low vertical resolution and a discrete distribution of fine to medium sands. Most existing 3D sediment transport models take into account the suspension load to compute the transport of non-cohesive sediment (e.g., [27][28][29][30]). Those models use classical erosion flux formulations from the literature with various numerical implementations of the advection-diffusion equation (e.g., [21,25,31]). However, the capacity of those models to reproduce sand transport rates at equilibrium for a wide range of sand grain sizes is generally not documented. In this paper, we point out that care should be taken to validate the transport rate at equilibrium for a wide range of sand size when using the advection-diffusion formulation. Indeed, we show in this study that existing erosion flux formulations are not necessarily consistent with 3D low vertical resolution models.
Here, we develop and validate a 1DV model (one-dimension vertical) of suspended sand transport (i.e., using an advection-dispersion equation with no bed load formulation included) for sand diameters D in the range of 100-500 μm. Modelled transport rates are compared to both data and parametric formulations reported in other studies under both wave and current conditions. The paper is structured as follows: in the Methods Section, the 1DV model used in this study is presented, before introducing the data and engineering formulations that are used for comparison. In the Results Section, the initial simulations are compared with classical equilibrium sand transport rates. Then, a new formulation of the erosion law is proposed and evaluated. A final discussion concludes the paper.

Methods
In this section, we first describe the one-dimension vertical (1DV) model that was specially developed for sand transport studies and later introduce the validation strategy.

1DV Model
This model is used as a research model with the goal of transferring the technique in the future to any 3D model used for coastal applications (including those at the shelf scale). Thus, vertical resolution should not be too high (10 layers in this study), which requires specific calculations (explained below) for application to sand transport.

Advection-Diffusion Equation
In order to compute sand transport, the following advection-diffusion equation is solved by the 1DV version of the Siam (simulation of multivariate advection) model [32,33]: where C is the suspended sand concentration, Kz the vertical turbulent diffusion determined from the Prandtl mixing length formulation and Ws the settling velocity. At the bottom, the boundary condition reads: where E and Fd are respectively the source and sink terms (erosion and deposition) at the bottom boundary. E and Fd are described hereafter.
Hereafter, C is always associated with the 1DV model concentration (averaged over the model layer), whereas c refers to a concentration resulting from an analytical computation.

Deposition Flux
The deposition flux is computed as: where bot c is the near-bed suspended sediment concentration and Ws is the settling velocity deduced from the sand grain characteristics following Soulsby [1]: 2 30 5 [(10 36 1 049 ) 10 36] with: where ν is the water viscosity, s the relative density (s = ρs/ρ), ρs the sediment density, ρ the water density and g the acceleration of gravity. Most often, the near-bed concentration cbot is assumed to be similar to the concentration Ck = 1 modelled in the bottom layer (whose index k is 1) (Figure 1), which is reasonable when the settling velocity is small. When dealing with sand and bottom layer thicknesses in the range 0.1 m to a few metres, this assumption is no longer valid, and cbot has to be extrapolated from Ck = 1. This can be achieved by coupling an analytical model for the concentration profile to the numerical model. Assuming a configuration near equilibrium conditions, an analytical Rouse-like profile can be used in the bottom boundary layer (e.g., [14,17,34]). However, the near-bed concentration generally becomes infinite when the distance from the bed to the location where it is computed tends to zero. Consequently, it is necessary to set a distance aref from the bed where cbot is being calculated (Figure 1). The reference height is often considered to be twice the sediment diameter (e.g., [25,31]). In our case, aref is chosen to be equal to the grain size roughness 2 5 ss k D   [35]. This corresponds to an equivalent Nikuradse roughness length for a flat sandy bed, considering that under this level, any concentration profile has probably no significance for suspended sediment transport. The method to calculate cbot is dependent on the choice of the analytical models chosen for both velocity and concentration profiles in the boundary layer. Those analytical models are presented below.

Bottom Boundary Layer Model and Near-Bed Concentration
We follow the concept of the double bottom boundary layer of Smith and Mc Lean [22], which has been found to accurately reproduce the partition between the skin shear stress and the total shear stress [36]. When bed forms are present, the flow is split into two layers above the bottom, the internal boundary layer (IBL), which is controlled by skin roughness, and the outer boundary layer (OBL), which is controlled by bedforms. The velocity profile u(z) in both layers is assumed to follow a logarithmic profile in agreement with the von Karman-Prandtl theory ( Figure 1): where sf u  is the friction velocity relative to the IBL, u  is the friction velocity relative to the OBL, z0s the skin roughness (z0s = kss/30 deduced from the grain size; see Appendix A), z0f the total roughness (z0f = ksf/30 deduced from bedform geometry; see Appendix B) and κ the Von Karman constant (= 0.4). It can be noticed that this velocity profile is in agreement with the turbulence closure of the 1DV numerical model (mixing length), but would also have been compatible with other turbulence closures, as most of them are similar in the bottom boundary layer, provided the Reynolds number remains high. Although the deposition and erosion fluxes may not be balanced, the concentration profile in the bottom layer is assumed to be an equilibrium Rouse profile both inside and outside the IBL (Figure 1). Inside the IBL, vertical mixing is caused by the skin shear stress ( sf  ), whereas it is caused by the total hydraulic shear stress ( ff  ) inside the OBL, so that concentration profiles differ ( Figure 1). The method for computing hydraulic shear stress and its skin friction component, under wave and current forcing, is given in Appendix A. Skin roughness z0s is assumed to only depend on the sediment size, while the form roughness is likely to change in response to the forcing, as the bedform itself depends on the forcing. A bedform predictor is used (Appendix B). When sheet flow is predicted (cf., Appendix B), IBL and OBL merge. We approximate the IBL thickness (zi) by the form roughness ksf ( Figure 1). When wave action is important for sediment motion, we assume that the wave boundary layer concerns the whole bottom layer in which the extrapolation is applied. Thus, following Soulsby [1], vertical mixing of sand induced by turbulent wave motion is accounted for by considering the maximum wave + current shear stress during a wave period in the Rouse number expression in both the IBL and the OBL.
The corresponding concentration profile can be written: and the Rouse number within the IBL defined as [1]: with c(zi) computed with Equation (8) for the sake of profile continuity and the Rouse number within the OBL defined as [1]: Finally, to determine the near-bottom concentration cbot from the concentration Ck = 1 computed by the numerical model (Equation (1)) in its bottom layer (k = 1), it is assumed that this concentration is representative of the average concentration within the bottom layer of thickness δzk = 1. The average concentration obtained from the integration of the analytical concentration profile (Equations (8) and (10)) over the bottom layer should then match this representative concentration Ck = 1 (i.e., . This leads to: with: and: 2

Suspended Sediment Horizontal Flux
The transport rate Qk is calculated by the model in each layer following: where ρs, Uk, Ck and δzk are respectively the sediment density, the velocity, the sediment concentration and the thickness of the layer k ( Figure 1). Note that since the eddy diffusivity is based on the Prandtl mixing length formulation, velocity profile is logarithmic. Inside the bottom layer, both the velocity profile and suspended concentration profile are subject to large gradients. To account for such gradients, especially when the vertical resolution is low, it is proposed to correct Q1 in Equation (15). A similar correction procedure has been previously proposed by Waeles et al. [34]. The transport rate thus needs to be corrected to compensate for the low vertical resolution.
The corrected flux can be expressed as where the corrective factor f for the transport rate is then computed as: with c(z) defined by Equations (8) and (10) and u(z) defined by Equations (6) and (7) (Figure 1). In Equations (6) and (7), sf u  and u  are computed from the velocity Uk = 1 in the bottom layer (i.e., at z = δz1/2), ensuring the adequacy of the coupling between the numerical model and the analytical model. The definite integrals in Equations (12) and (16) are computed using a numerical integration at each time step. For this purpose, the bottom layer is discretised in thin sub-layers (100 sub-layers in this study) with a thickness decreasing toward the bottom to account for the increasing concentration gradient toward the bottom. These numerical integrations require only a little computational work in comparison with the whole model.

Erosion Fluxes
The closure of the model requires that the erosion flux be specified. Several formulations based on experiments are proposed in the literature. However, the horizontal transport rates resulting in fine from those formulations of the erosion flux seem not to have been validated. The validation of sand transport rates using erosion flux formulations will therefore constitute the core of our study.
Non-cohesive sand transport is characterized by a rapidly-reached equilibrium state where the deposition flux balances the erosion flux. Based on Equation (3), when equilibrium is reached, the deposition flux reads: where ca represents the so-called reference concentration. At equilibrium, E = Fd, and the erosion flux can then be expressed as s a E W c  . More generally, some authors propose to express the erosion flux following this latter equation even when equilibrium is not reached (e.g., [24,25,31]). When the erosion boundary condition is expressed with a reference concentration, attention has to be paid to which reference height it refers: naturally, it has to be the same as the one selected for the deposition term. In this study, four formulations, detailed in Appendix C, were tested. The formulations of Le Hir et al. [23], van Rijn [21], Engelund and Fredsøe [25] and Zyserman and Fredsøe [31] were implemented in the 1DV model and respectively denoted Siam-ERODI, Siam-VR84, Siam-EF76 and Siam-ZF94. In those model configurations, the default reference height (aref = 2.5D) was changed to match the reference height corresponding to the erosion boundary condition (if prescribed).

Transport Rate Validation Strategy
In this section, we present the various data and sand transport rate formulations used to validate the model for different hydrodynamic conditions. All model and data abbreviations are summarized in Table 1. Many practical models that provide the transport rate at equilibrium are proposed in the literature. To validate our model in the case of current only, three well-known formulations of total transport were selected. The expressions of Engelund and Hansen [37], Yang [38] and van Rijn [39], respectively denoted EH67, Y73 and VR84, are detailed in Appendix D. These engineering or practical models (in the sense that they are straightforward computations for sediment transport) were tested by van Rijn [60] for a set of 1260 data acquired in rivers, estuaries and flumes. These experiments were conducted using sand diameters D ranging from 100 to 400 μm and for fluid velocities ranging from 0.4 to 2.4 m/s. According to van Rijn [60], formulations EH67 and VR84 agree relatively well with the data, whereas formulation Y73 underestimates the transport rate in the case of large water height. Moreover, formulations EH67 and Y73 need to determine hydraulic roughness in order to calculate the friction velocity u  , which may considerably modify total transport rates. Calculations using these engineering models (only) are made under flat bed conditions (i.e., the roughness height 2 5 s k D   [1]), assuming that no ripples are formed or at least that the required friction velocity is related to skin friction. Formulation VR84 is a simplified method derived from the one proposed by van Rijn [3] in which the hydraulic roughness is not required.

Wave and Current
In order to validate our model in the case of combined wave and current forcing, the results of Davies et al. [19] and Camenen and Larroudé [51] are used. These studies compared transport rates deduced from engineering and research models using a wide range of varying parameters: ). In our study, we use the results of Davies et al. [19], which include two 1DV models (with high vertical resolution), the STP model [40] and a TKE model [41,42], as well as engineering formulations called BIJKER A [43,44], SEDFLUX [1,45,46], Dibajnia and Watanabe [47], TRANSPORT [48] and Bagnold-Bailard [49,50]. Like in Davies et al. (2002) [19], our model is tested and compared to other models for a set of five hydrodynamic configurations listed in Table 2. were obtained during strong currents conditions without waves [55]. Laboratory data acquired in oscillating water tunnels with both waves and current (about 260 data) come from Al Salem [56] or from Ribberink and Al Salem [57] (these are called alsalem), from Dohmen-Janssen [58] (called janssen) and from Dibajnia and Watanabee [47] or Dibajnia [59] (called dibajnia). In the present study, the same data are used to validate our model, and the model performance is evaluated using the approach of Camenen and Larroudé [51]. The modelled sediment flux (qs,num) is considered to be acceptable when it approaches the experimental flux (qs,data) within a factor 2 (here, "a factor n" means between 1/n times and n times the measured value). For each model, the percentage of values predicted within a factor 2 (i.e., error lower than 50%) and within a factor 1.25 (i.e., error lower than 20%) is computed for "current only" data and "wave-current" data (indicated as Cc50, Cc80, Cw50 and Cw80, respectively). When the experimental flux is negative (current and wave are in opposite directions), the absolute value is shown on the graph if the predicted flux direction is correct. If not, the error is considered to be infinite.

Results
The Siam 1DV model described previously was used to simulate equilibrium horizontal transport rates for various erosion flux formulations. In the following section, we evaluate the model results in a simple no-wave configuration.

Evaluation of the 1DV Model Transport Rates in the "Current Only" Case
Results were compared to total transport rates computed with engineering formulations (EH67, Y73 and VR84) and a few data mentioned in the Methods Section. Tests were run for currents varying from 0 to 2 m/s (without waves) and various sand diameters (from 100 to 500 μm) and with a water depth of 7.2 m. The comparison in the case of a current velocity of 1.8 m/s is presented in Figure 2. In agreement with numerous studies, transport rates provided by engineering models differ considerably from one another, as they fluctuate by about one order of magnitude ( Figure 2). The few data corresponding to this "current only" experiment fall within the range of values proposed by the three formulations. In agreement with van Rijn [60], formulation Y73 seems to underestimate the transport rate, at least for this particular experiment. Besides, it appears that, whatever erosion flux formulation is used, the transport rates computed by Siam 1DV (hereafter, "1DV transport rates") do not match the range of engineering and experimental fluxes ( Figure 2). Comparisons made for other current velocities (results not shown) confirmed this result. Finer sand fluxes are consistently overestimated, whereas coarser sand fluxes are consistently underestimated.

New Empirical Formulation of Erosion Flux
The previous comparison showed that the erosion fluxes that are given in the literature, after integration by the 1DV model and once equilibrium is reached, lead to largely overestimated horizontal fluxes of fine sand. On the contrary, modelled fluxes are strongly underestimated for coarser sand.
In order to figure out the reason for those discrepancies, we computed the analytical transport rate obtained from the product of a simple Rouse profile (Equation (8)) and a simple logarithmic velocity profile (Equation (6)) in the water column. This test is performed for a mean current velocity of 1.8 m/s and a water depth of 7.2 m (similar to the test in Figure 2). We consider a flat bottom and aref = ksf = 2.5D and z0 = ksf/30. The current velocity, the sediment concentration and the transport rate are computed for three grain sizes (100, 200 and 500 μm) and reported in Figure 3.   cbot). The Rouse number Rsf is indicated for each relative concentration profile. (c) Relative transport rate (i.e., transport rate normalized by the reference concentration cbot). The depth integrated transport rate is indicated next to each profile.
In this simple configuration, it appears that for the same reference concentration cbot, the total transport rates for 100 μm and 200 μm sands are two orders of magnitude different (Figure 3). Further, there is only one order of magnitude between the transport rate for 200 μm and 500 μm. According to formulations EH67, VR84 and Y73, transport rates from 100 μm to 500 μm should be within the same order of magnitude ( Figure 2). With this simple model, in order to balance the transport rate for the three different grain sizes, a large range of reference concentrations is required, especially between 100 and 200 μm. Since at equilibrium, the erosion flux is proportional to the settling speed, which increases with diameter, and to the reference concentration ( s bot E W c  ), a large dependency of the erosion flux with the diameter is needed. In the erosion flux formulation that we used previously (ERODI and VR84), there is less than one order of magnitude difference from 100 μm to 500 μm ( Figure 4). This is likely to explain the discrepancies observed in the Siam 1DV model while using erosion flux formulation from the literature. Therefore, a new erosion flux has been empirically adjusted. Assuming an erosion law in the form of following Equation (18), we determined the erosion rate in the Siam 1DV that fit the best the range of transport rates calculated with VR84, EH67 and Y73 models for various sand diameters (from 100 to 500 μm) and velocities (from 0 to 2 m/s). More details on the determination of the erosion law are given in Appendix E. The result of the fit leads to the following formulation of the erosion law: with: with a1 = 18,586 kg·m −3 ·s −1 and b1 = −3.38 kg·m −2 ·s −1 , and otherwise: with a2 = 59,658 kg·m −3 ·s −1 and b2 = −9.86 kg·m −2 ·s −1 .
The power law α of the relative excess shear stress is adjusted to 0.5. This empirical formulation involves a coefficient E0 that depends on the sand diameter ( Figure 4). The erosion factor E0 increases from 0.03 to 21 kg·m −2 ·s −1 for fine sands between 100 and 500 μm. The critical shear stress τce is computed from the classical critical Shields' parameter for sand movement.

Validation of the New Erosion Flux Formulation
The result of the fit that led to the formulation of the erosion law is presented in Figure 5 for three different velocities (0.5, 1 and 1.8 m/s). The transport rates obtained with the new erosion flux formulation match well with the engineering formulations VR84, EH67 and Y73 in the case of current only. In the following sections, we further assess the capacity of the model to transport sand in suspension for a range of different configurations.

Overall Performance
For the purpose of comparison, Figure 6a shows the results obtained by Camenen and Larroudé [51] with the van Rijn [53] formulation, which is considered to yield reasonable results (Camenen and Larroudé [51]). According to Figure 6b, the transport rates computed with the advection/diffusion model and our erosion law (Equations (18)- (20)) compare reasonably well with the empirical laws used to directly assess equilibrium solid fluxes. Measured transport rates were predicted within a factor 2 in 67% of cases with current only and in 35% of cases with both waves and current (Table 3). However, an overestimation is sometimes observed in the case of the scheldt data (concerning strong currents and fine sands). Under wave-current interaction, the model also overestimates data sometimes. A considerable error is likely to appear when waves become strongly asymmetric, as the model was not designed to account for the aforementioned process. Actually, at least one third of the data corresponding to wave forcing was acquired under asymmetric waves, inducing significant transport rates without current or transport opposite the current direction. Table 3. Comparison between various models and data described by Camenen and Larroudé [51] and in this study. See the explanations of symbols Cc and Cw in the Methods Section. Adapted from Camenen and Larroudé [51].  Figure 6. (a) Comparison of total sand transport rates between the Van Rijn [53] formula (qs,num) and experimental data (qs,data). From Camenen and Larroudé [51]. (b) Comparison of total sand transport rates between Siam 1DV (qs,num) and experimental data (qs,data). See the explanations of the legend and symbols Cc and Cw in the Methods Section. Values predicted within a factor 2 sit in between the two dashed lines (qs,num = 0.5qs,data and qs,num = 2qs,data).

Figure 7.
Comparison of sand transport rates as a function of sand diameter in the case of prevailing waves (a); the case of wave and current (b) and the case of current only (c). Bold black curves, corresponding to Siam 1DV transport rates using the new erosion law (Equations (18)- (20)), are added to other models' results (copied from Camenen and Larroudé [51]). Uc is the mean current velocity; h the water depth; Tw is the wave period; and Uw is the wave orbital velocity at the bottom.

Validation for Varying Grain Size
The effect of grain size variation was investigated with the Siam 1DV approach, and results were added to the comparison of Camenen and Larroudé [51] (Figure 7). Results are coherent with both data and other transport rate formulations in the case of prevailing waves, waves and current or prevailing current (Figure 7). The large disparities in transport rates for finer sands highlight the difficulty in validating the model for fine sands. This is due to the fact that most studies focus on sand around  Table 2) (sand diameter 250 μm). Bold black curves, corresponding to Siam 1DV transport rates using the new erosion law (Equations (18)- (20)), are added to other models' results (copied from Davies et al. [19]).

Validation for Varying Hydrodynamic Conditions
The robustness of the model for a wide range of wave and current conditions was finally checked by adding Siam 1DV results to the comparison study presented by Davies et al. [19] (Figure 8). These comparisons confirmed the good behaviour of the model for different hydrodynamic conditions. However, in the case of strong waves, the Siam model provides among the highest transport rates of the comparison. This matches with the previous comparison with data that showed a slight overestimation in the case of waves (Figure 6b).

Relevance of Analytical Refinements in the Bottom Layer
We tested the impact of the vertical resolution on the transport rate. First, we used the model without either extrapolating the near bed concentration from the model concentration in the bottom layer (for deposition) or correcting the expression of horizontal flux in the bottom layer (i.e., f = 1 and Fd = WsCk = 1). In this configuration, the sand fluxes in the case of current only are highly dependent on the vertical resolution ( Figure 9). For two different resolutions of the bottom layer (20 cm and 5 m), the transport rates (grey lines) vary by about one order of magnitude. On the contrary, transport rates calculated with the model presented in this study are nearly independent of the resolution of the model (Figure 9). We can however note that some discontinuities, related to the appearance of the sheet-flow regime, are encountered. The change of roughness under sheet-flow and its impact on the suspended flux depend on the model resolution in our model.

Discussion
We tested the model, using several erosion flux formulations from the literature, against data and equilibrium transport rate formulations. None of these erosion flux formulations tested in the model led to acceptable simulations of sand transport rates. Three hypotheses could be suggested to explain this inconsistency:  Bed load has not been accounted for in our model. For the experiment leading to Figure 2, adding the bed load in the model leads to the same conclusion. For instance, using VR84, which gives both suspended load and bed load, it appears that, whatever the velocity, the bed load only accounts for 6%-22% of the total transport from the finer to the coarser sand (following VR84, the ratio of the suspended load to the total load is independent of the velocity). Moreover, for the "current only" experiment ( Figure 2), the ratio / s u W  is always greater than one, which means that the suspended load is the predominant mode of transport [61]. The bed load is therefore not expected to explain the discrepancies observed when using classical erosion flux formulations, at least for this experiment.  The 1DV model could suffer from weaknesses, mostly in the way the boundary layer is formulated. It is for instance not fully clear if the Rouse profile holds very close to the bottom. Different tests were conducted with varying parameterization (reference height, hindered settling velocity or the parameterization of vertical mixing were modified) and led to the same conclusion. A simple Rouse concentration profile (assuming a single bottom boundary layer) has also been tested, but did not enable us to match the total flux range given by engineering models. Further, we cannot take into account the impact of sediment concentration on settling velocity near the bottom in the model. This is due to the fact that we are using an analytical solution for the bottom layer of the model (a Rouse profile), which is only valid for constant settling velocity (with depth). We acknowledge that this could induce biases in the transport rates computed by the model.  The erosion fluxes proposed in the literature could be inconsistent with our modelling strategy.
We suggest here that this hypothesis is the most relevant one. Our approach was thus to search for an empirical erosion law able to produce results that match the horizontal flux range given by empirical engineering models.
The ability to transport sand with an advection-diffusion equation was then achieved by introducing a new empirical erosion law. It should be noted that this erosion law is consistent with a reference height (elevation where the deposition flux is evaluated) of 2.5D. The erosion flux involves a power law of the excess shear stress multiplied by a nonlinear function of the sand diameter. The power α of 0.5 is in agreement with the formulation proposed by Waeles et al. [13], but is lower than in other classical formulations [21,23]. We tested different power α, in agreement with other formulations [21,23], but this did not lead to acceptable results. Indeed, in Figure 8a, for instance, the shape of the curve was not matching the expected shape, with a much steeper increase of the transport rate with current velocity than expected. The variation of E0 with D is much more amplified than in other formulations ( Figure 4). We acknowledge that this erosion flux formulation may be dependent on the conceptual model of the boundary layer that we adopted in this study. Transposing this erosion flux formulation to other studies would therefore require using the same boundary layer conceptual model. However, we suggest that the strong variation of E0 with D is not specific to our model and that the same behaviour would be observed for any other comparable model. This aspect would require further investigations. We also suggest that the same validation approach should be carried out when using multi-class sediment transport models using advection-diffusion equations for sand. Further, due to the strong variation of E0 with D and the discontinuity in the formulation around 180 μm, we acknowledge that some discontinuities in the transport rates can be produced (e.g., Figure 7c). This, however, should not be an issue for implementation in 3D models using a discrete form of E0 (i.e., using only a few sand classes).
In this study, we looked at transport rates at equilibrium, so that the coupling between an unsteady numerical model and a stationary analytical model remains coherent. Further investigations should, however, be carried out to evaluate the impact of this coupling while studying unsteady flows. We also acknowledge that, since the model does not account for bed load, it has been validated for medium sands, probably thanks to some compensation in the fitted erosion law. If our approach seems reasonable in the case of fine sand particles, it would probably reach its own limit when the bed load is prevailing.
The strength of the 1DV approach compared to engineering or practical models is that it computes the sediment concentration profile in the whole water column. In the present study, we only validated the total sediment transport at equilibrium. In the future, it would also be interesting to compare the sediment concentration profile with field data, especially close to the sediment bed, where the analytical model is coupled, and in non-stationary conditions.
In parallel with the validation of the suspended sediment transport model at equilibrium (this study), the model was implemented in a 3D multi-class sediment transport model [12]. It was tested at shelf scale using one class of mud and two classes of sand (160 and 380 μm). For the sand fractions, the erosion flux formulation and the conceptual model of the boundary layer presented in this study were adopted, and the sediment dynamics in that modelling study proved to be consistent with in situ measurements [12].

Conclusions
Transporting sand in suspension is necessary in sediment transport models dealing with cohesive sediment, in order to allow the management of mixtures of sand and mud. The initial objective of this work was to evaluate the ability to predict fine to medium sand horizontal fluxes under varying wave and current conditions by using an advection-diffusion model. It has been shown that using classical erosion rate formulations from the literature in our model did not lead to consistent results. Those classical erosion laws should not be used in advection-diffusion models without a full validation procedure involving a range of sand diameters and hydrodynamic conditions. The ability to transport fine sands in our modelling framework was only made possible by introducing a new empirical erosion law. Further, the transport rate computed from the advection-diffusion model had to be independent of the vertical resolution, so that it can be implemented in regional models with operational applications. This was successfully achieved by using a correction procedure for extrapolating the near-bed concentration (in the deposition flux) from the concentration computed in the bottom layer. In addition, the actual horizontal flux in this bottom layer deserved a correction to account for very strong inverse gradients of velocity and concentration in the boundary layer. It is suggested that in a three-dimensional model, such a flux correction should be applied when solving the advection-diffusion equation for suspensions.
Using the new empirical erosion flux formulation, the model demonstrated its capacity to simulate transport rates for sand diameter in the range 100-500 μm. The proposed erosion law was evaluated as best as possible for a wide range of different hydrodynamic conditions and grain sizes. This has been done despite the difficulty involved in qualifying sand transport models, illustrated by the very scattered results of classical sand transport formulations [19,41,51,60,62]. This evaluation showed that our model compares well with various other sand transport models. While comparing the model with data, we showed that measured transport rates were predicted within a factor 2 in 67% of cases with current only and in 35% of cases with both waves and current. In comparison with the results obtained by Camenen and Larroudé [51], who provided the same indicators for several practical transport rate formulations (whose means are respectively 72% and 37%), our model gives reasonable results. Finally, the model presented in this study has been formulated for well-sorted, non-cohesive sediments. In the future, it would have to be adapted in the case of heterogeneous sediment composition to reach our longer term goal. -Ub and A are respectively the orbital velocity and half excursion above the wave boundary layer, , with T the wave period).
The formulation of Soulsby [1] is used to account for non-linear wave-current interactions: where τm represents the mean (wave-averaged) shear stress in the direction of the current and φ is the angle between current and wave directions. τ is the maximum shear stress generated during a wave period and is either used for the expression of the skin sf  or the hydraulic shear stress ff  , depending on the roughness used (ks = kss or ksf).

B. Bedform Predictor
The form roughness ksf is predicted following Yalin [  is the shear stress for which a "sheet-flow" regime appears and the bed is flat. We finally needed to determine E0. We found out that E0 depend strongly on D. Therefore for each different diameter (ranging from 100 to 500 μm) we computed the E0 value minimizing the RMSE between the Siam 1DV transport and the average transport obtained from VR84, EH67 and Y73 for current velocities ranging from 0 to 2 m/s. A curve of best fit has then been computed to present E0 as a continuous function of D.