Modelling Fine Sediment Dynamics: Towards a Common Erosion Law for Fine Sand, Mud and Mixtures

Abstract

This study describes the building of a common erosion law for fine sand and mud, mixed or not, in the case of a typical continental shelf environment, the Bay of Biscay shelf, characterized by slightly energetic conditions and a seabed mainly composed of fine sand and muddy sediments. A 3D realistic hydro-sedimentary model was used to assess the influence of the erosion law setting on sediment dynamics (turbidity, seabed evolution). A pure sand erosion law was applied when the mud fraction in the surficial sediment was lower than a first critical value, and a pure mud erosion law above a second critical value. Both sand and mud erosion laws are formulated similarly, with different parameters (erodibility parameter, critical shear stress and power of the excess shear stress). Several transition trends (linear or exponential) describing variations in these erosion-related parameters between the two critical mud fractions were tested. Suspended sediment concentrations obtained from simulations were compared to measurements taken on the Bay of Biscay shelf with an acoustic profiler over the entire water column. On the one hand, results show that defining an abrupt exponential transition improves model results regarding measurements. On the other hand, they underline the need to define a first critical mud fraction of 10 to 20%, corresponding to a critical clay content of 3–6%, below which pure sand erosion should be prescribed. Both conclusions agree with results of experimental studies reported in the literature mentioning a drastic change in erosion mode above a critical clay content of 2–10% in the mixture. Results also provide evidence for the importance of considering advection in this kind of validation with in situ observations, which is likely to considerably influence both water column and seabed sediment dynamics.

1. Introduction

The transport of fine sediments can be assumed to mainly occur in suspension. Suspended transport is generally simulated by solving an advection/diffusion equation, assuming that sediment particles have the same velocity as water masses, except the vertical settling component (e.g., [1]). Such an equation involves sink and source terms at the bed boundary, which are deposition and erosion fluxes under conditions defined by the hydraulic forcing and the behaviour and composition of both the suspended and deposited sediments. This means that dealing with fine sediment dynamics, whether or not cohesive, requires the formulation of an erosion law. Such an erosion law should be applicable for fine sands (and even medium sands under strong shear stresses) as well as for mud, and naturally for mixtures of sand and mud. Despite these similarities between the suspended transport of fine sand and mud, their erosion processes have generally been investigated separately due to their contrasting behaviours.

An abundant literature is available on cohesive sediments, based on the existence of a critical shear stress for erosion (τ_e in N·m⁻²) and several empirical relationships between the erosion flux E (in kg·m⁻²·s⁻¹) and excess shear stress (i.e., the difference between the actual shear stress (τ in N·m⁻²) and the critical value τ_e, either normalized by the latter or not). The parameters of such an erosion law involve bed characteristics, which may concern electrochemical forces, mineral composition, and organic matter content [2], or pore water characteristics [3]. They also depend on the consolidation state [4], and may be altered by biota effects [5]. Bulk density has often been proposed as a proxy for characterizing the bed, but the plasticity index has also been suggested [6], as well as the undrained cohesion and the sodium adsorption ratio [3].

Literature on sand erosion laws is scarce, partly because of experimental difficulties linked to simultaneous settling and resuspension processes, and partly because sand transport has mainly been considered through the formulation of transport capacity, even in the case of suspension. The need to simultaneously simulate transport of mud and (fine) sand, and their mixture, updates the need for an erosion law for fine sand, often named pick up functions [7]. Pick up functions link the erosion rate with the particle characteristics (size and density) and the shear stress (or, equivalently to the Shields parameter) [8,9]. As for cohesive sediment erosion laws, the dependence on the forcing stress is either absolute (~τ) or relative (~τ/τ_e), involving (or not) a threshold value τ_e.

The analogy between erosion formulations for sand and mud can be highlighted, especially considering the Partheniades-Ariathurai law [10,11] for cohesive sediments, the more usual in numerical models. The erosion law can then be written as

$$E=E_0·{\left(\tau /{\tau }_e-1\right)}^n,if \tau \ge {\tau }_e$$

(1)

where E₀ is an erodibility parameter in kg·m⁻²·s⁻¹ and n a power function of the sediment composition. In Equation (1), E₀ and τ_e are functions of the sediment composition and its consolidation state in the case of cohesive sediments, or functions of the particle diameter and density in the case of non-cohesive sediments. Assuming a similar erosion law for the whole range of mixed sediments, the problem becomes the assessment of the critical shear stress for erosion on the one hand, and the erosion factor E₀ on the other hand, in the full transition range between cohesive and non-cohesive materials. Following Van Ledden et al. [12], the proxy for such a transition could be clay content. An alternative proxy could be mud content (hereafter referred to as f_m), considering that the clay to silt ratio is often uniform in a given study area [13].

Experiments have provided evidence that the resistance to erosion (i.e., critical shear stress τ_e) increases when mud is added to sand, either because of electrochemical bonds which take effect in binding the sand grains or because a cohesive matrix takes place between and around sand grains (e.g., [6,7,14,15,16,17]). Literature on erosion rate for mixed sediments is much less abundant, but a significant decrease (several orders of magnitude) in the erosion rate with mud content is most often reported (e.g., [14,17,18]). For instance, Smith et al. [17] presented laboratory measurements showing a decrease of about two orders of magnitude when the clay fraction in the mixture increased from 0% to 5–10%, and up to one order of magnitude more when the clay fraction increased from 10% to 30%.

From the point of view of modelling, Le Hir et al. [1] reviewed the different approaches developed in the last 20 years to manage sand/mud erosion in numerical process-based models. A non-cohesive and a cohesive regime separated by one or two critical mud fractions were commonly introduced and simulated by specific erosion laws. Le Hir et al. [1] and Waeles et al. [19], partly followed by Bi and Toorman [20], used a three-stage erosion law. Below a first critical mud fraction f_mcr1, they considered a non-cohesive regime where the erosion flux of any class of sandy and muddy sediments remains proportional to its respective concentration in the mixture, but is computed according to a pure sand erosion law (with a potential modulation of the power applied to the excess shear stress). Starting from a value characteristic of a pure sand bed, critical shear stress in this first regime is either kept constant [1], or linearly increases with the mud fraction f_m [19,20]. Above a second critical mud fraction f_mcr2, these authors defined a cohesive regime: Waeles et al. [19] and Le Hir et al. [1] formulate the erosion law using the relative mud concentration (the concentration of mud in the space between sand grains), considered as more relevant than the mud density in the case of sand/mud mixtures, which was in agreement with observations from Migniot [21] or Dickhudt et al. [22]. Between f_mcr1 and f_mcr2, they ensured the continuity by prescribing a linear variation of E₀ and τ_e between non-cohesive and cohesive erosion settings. Carniello et al. [23] used a two-stage erosion law built by Van Ledden [24]: below a critical mud fraction f_mcr, the erosion factor of the sand fraction is steady and the one for the mud fraction varies slightly according to the factor 1/(1 − f_m), while above f_mcr the erosion rate is the same for sand and mud fractions and logarithmically decreases according to a power law. Regarding the critical shear stress for erosion, it first slightly increases with f_m and then varies linearly to reach the mud shear strength. Dealing only with the critical shear stress, Ahmad et al. [25] proposed an alternative to the Van Ledden [24] expression: without any critical mud fraction, τ_e varies linearly with f_m for low values of f_m and more strongly for high f_m values, using a parameter representing the packing of the sand sediment in the mixture. Generally speaking, the transitional erosion rate between the two regimes is poorly documented.

The aim of this paper is to fit an erosion law for mixed sediments to be applied in environments dominated by fine sediment, such as continental shelves, including our area of interest, the Bay of Biscay continental shelf (hereafter referred to as BoBCS). In these environments that develop outside the coastal fringe where wave impact is higher and where stronger tidal currents may take place, a common sedimentological feature is a mixture of fine sands and mud. In addition, bathymetric gradients are often gentle so that the shear stress gradients are likely to be small: as a consequence, the critical shear stress for erosion of surficial material in equilibrium with such environments is also likely to be low contrasted. Last, suspended sediment concentrations (hereafter referred to as SSC) are rather low in such “deep” coastal waters (at least they are much lower than in shallow waters), and sediment exchanges remain very small: as a consequence, freshly deposited sediment does not consolidate under its own weight, and erosion is so small that the surficial sediment can hardly become over-consolidated. This means that in this context, consolidation appears as a second order process, and that when the surficial sediment is muddy, its shear strength never reaches high values. On the contrary, the surface shear strength is likely to remain close to typical values for fine sand in the order of 0.1–0.2 Pa (for sands, this range corresponds to the decreasing trend of the Shields parameter with diameter, when the latter is below 0.6 mm).

Due to difficulties in measuring erosion fluxes (and not only shear strengths) for sand and mud mixtures, the assessment of the erosion law for fine sediment of the BoBCS is achieved by means of continuous measurements in the field, completed by a realistic hydro-sedimentary model of the specific area for testing the erosion law and its parameterization. In this way, the paper also proposes a methodology for fitting the erosion law in a poorly known natural environment. The paper is organized as follows. Section 2 gives the strategy for assessing the erosion law and describes measurements acquisitions on the one hand, and the modelling framework on the other. The description of the model includes the main features of hydrodynamic and sediment models, and the validation of hydrodynamics computation. Section 3 details the successive steps which led to the building of a new erosion law for sand/mud mixtures. Section 4 presents model results obtained in different erosion settings along with an assessment of their relevance with respect to observations. Section 5 discusses the main results in light of previous studies, followed by a short conclusion in Section 6.

2. Strategy and Modelling Background

2.1. Strategy for Assessing an Erosion Law, and Its Application to the BoBCS

Most often erosion laws have been investigated by comparing the critical erosion shear stress and erosion fluxes deduced from the tested law with measurements in a laboratory flume. However, this process has some difficulties. When the sediment is placed or settled in a flume, its behaviour is likely to differ from natural sediment one, especially when cohesive properties begin to develop. In addition, while critical shear stresses can be determined straightforwardly, the erosion flux of the sandy fraction is questionable, mainly because computations generally apply to rough erosion, while measurements often concern a net erosion flux, resulting from this rough erosion and possible simultaneous deposition, which are not easy to characterize. This concern does not apply to the mud fraction.

Here, an alternative methodology is proposed. Selecting an environment where local sediment is representative of the study area, a continuous measurement of suspended sediment concentration and local forcing is used to test the erosion law by simulating observed features by means of modelling. In the likely case of non-negligible horizontal gradients, a full 3D modelling is preferred, as it accounts for all processes, including erosion, deposition and horizontal advection. Deploying field measurements over several months increases the probability of investigating very different forcing conditions, especially if a winter period is selected. The use of an acoustic current profiler enables a simultaneous measurement of local forcing (both current patterns and waves) and resulting resuspension over the whole water column, by analysing the acoustic backscatter. The former observation is used for validating hydrodynamics simulations, while the latter can be compared to predicted SSC according to tested erosion parameters.

2.2. Measurements Used for Erosion Law Assessment and Model Validation

As in many continental shelves (e.g., [26,27]), the surficial sediment of the BoBCS is mainly composed of fine sand (~200 µm) and mud, since water depths exceed the depth of wave-induced frequent reworking, typically about 20–30 m, for the wave regime of the Bay of Biscay. The erosion law for fine sediment mixtures has been fitted by selecting measurements in a rather representative environment of the shelf, both in terms of sediment cover, hydrodynamic forcing and near bottom turbidity. A station located near the coast of southern Brittany, close to Le Croisic (hereafter referred to as LC station, Figure 1) meets these criteria: the local water depth is 23 m on average and the exposure to waves is attenuated by a rocky bank localised in the southwest nearby. Despite a tidal range of 4–5 m, tidal currents remain low, and flow is likely to be controlled by wind-induced currents [28]. The seabed is constituted by muddy sands (d₁₀ of 7.5 µm, d₅₀ of 163 µm, 5.1% clay (% < 4 µm) and 25% mud (% < 63 µm) contents), and exhibits some gradients around the station, with muddy facies to the north, and more sandy ones to the south (Figure 1).

For measurements, a mooring line was deployed over the whole water column for two months between 25 November 2007 and 31 January 2008. This period corresponded to typical winter hydrodynamic and meteorological conditions: mainly south-westerly winds, high rates of river discharge (e.g., Loire river: average flow rate of 1161 m³·s⁻¹, and maximum flow rate of 2240 m³·s⁻¹), and rather large swells (H_s (significant wave height) peaks > 3 m, T_p (peak period) 10–18 s).

Several instruments were placed along the mooring line. A MS5-Hydrolab sensor (OTT Hydromet, Kempten, Germany) providing temperature and salinity measurements (at 60 min intervals) was fixed 1.50 m below the sea surface. An upward looking acoustic Doppler profiler (Acoustic Wave And Current profiler, 1000 KHz, Nortek AS, Vangkroken 2, Norway; hereafter referred to as AWAC) was fixed on a structure placed near the seabed to record the backscattered acoustic intensity (at 30 min intervals), the intensity and direction of the currents (1 min of time-integrated data at 30 min intervals), water elevation (at 30 min intervals), as well as wave height, period, and direction (at 60 min intervals). The return echo was sampled over the entire water column with 25 cm-thick cells. In this configuration, the first sampling cell of the AWAC profiler was located 1.67 m above the seabed (including the frame height and the AWAC blank distance). At the same elevation, a turbidity sensor (WETlabs; Sea-Bird Scientific, Bellevue, Washington, DC, USA) provided hourly measurements for calibration. This turbidity signal was further transformed to SSC in mg·L⁻¹ using water samples. The backscatter index (BI) from the AWAC profiler was evaluated from the sonar equation [29], following the procedure described by Tessier et al. [30], in particular by considering the geometrical attenuation for spherical spreading, the signal attenuation induced by the water, and the geometric correction linked to the expansion of the backscattering volume with increasing distance from the source. Given that SSC derived from the turbidity sensor did not exceed 100 mg·L⁻¹, the signal attenuation caused by the particles was disregarded when estimating backscatter [30]. Then, an empirical relationship was established between the BI of first AWAC cell and SSC measurements of turbidity sensor, following Tessier et al. [30]:

$${10log}_{10}\left(SSC\right)=c_1·BI+c_2$$

(2)

We finally obtained a determination coefficient R² of 0.78 with c₁ = 0.42794 and c₂ = 32.8907. Changes in SSC concentrations in the water column could thus be quantified (Section 4).

2.3. Hydrodynamics Models (Waves, Currents)

2.3.1. Brief Description

Current patterns and advection of suspended sediments were computed with the three-dimensional MARS3D model (3D hydrodynamic Model for Applications at Regional Scale, IFREMER), described in detail by Lazure and Dumas [31]. This code solves the primitive equations under classical Boussinesq and hydrostatic pressure assumptions. The model configuration spreads from 40°58′39′′ N to 55° N in latitude, and from 18°2′1′′ W to 9°30′ E in longitude (Figure 1) with a uniform horizontal resolution of 2500 m (822 × 624 horizontal grid points). The vertical discretization is based on 40 generalized sigma layers. For water depths lower than 15 m, sigma levels are uniformly distributed through the water column. Above this depth, a resolution increase is prescribed in the lower and upper parts of the water column, following the formulation of Song and Haidvogel [32]. For instance, this setting leads to a bottom cell of 40 cm for a total water depth of 25 m. The numerical scheme uses a mode splitting technique with an iterative and semi-implicit method which allows simultaneous calculation of internal and external modes at the same time step. In our case, the time step is fixed at 150 s. Vertical viscosity and diffusivity for temperature, salinity and momentum were obtained using the k-ε turbulence closure of Rodi [33]. Simulations were performed with the realistic forcings detailed in

Table 1. Realistic forcings used in 3D simulations.

Forcing	Source
Initial & boundary conditions (3D velocities, temperature, salinity)	GLORYS global ocean reanalysis [34]
Wave (Significant height, peak period, bottom excursion and orbital velocities)	WaveWatch III hindcast [35]
Meteorological conditions (Atmospheric pressure, wind, temperature, relative humidity, cloud cover)	ARPEGE model [36]
Tide (14 components)	FES2004 solution [37]
River discharge (flow and SSC)	Daily runoff data (French freshwater office)

.

Regarding wave forcing, a wave hindcast database built with the WaveWatchIII (WWIII) model (realistic and validated configuration of Boudière et al. [35]) enabled the computation of bottom wave-induced shear stresses (τ_w in N·m⁻²). The wave-induced shear stress was computed according to the formulation of Jonsson [38], with a wave-induced friction factor determined following Soulsby et al. [39]. Then, the total bottom shear stress τ was computed from the estimated τ_w and from the current-induced shear stress (τ_c in N·m⁻²) provided by the hydrodynamic model, according to the formulation of Soulsby [40], i.e., accounting for a non-linear interaction between waves and currents. Both wave and current shear stresses were computed by considering a skin roughness length z₀ linked to a 200 µm sand, representative of the sandy facies widely encountered on the BoBCS (z₀ = k_s/30 = 2 × 10⁻⁵ m with k_s the Nikuradse roughness coefficient).

2.3.2. Hydrodynamic Validation of the Model

Figure 2 shows model validation regarding hydrodynamics at the LC station from 1 December 2007 to 18 January 2008 in terms of significant wave height (H_S), surface temperature and salinity (T_S and S_S, respectively), and current intensity and direction (Vel_INT and Vel_DIR, respectively).

Measurements highlight typical late autumn/winter energetic conditions. Average H_s were around 1.5 m and occasionally exceeded 4 m in stormy conditions. Peak periods (not illustrated here) ranged from 8 to 18 s (around 10 s on average). Even during “calm” periods, H_S values are generally no lower than 0.8 m. Figure 2a illustrates the ability of the model of Boudière et al. [35] to describe H_S over the period, with a root mean square error (RMSE) of 0.24 cm and a R² of 0.95. Measured and modelled T_S and S_S are illustrated in Figure 2b,c respectively, and demonstrate the correct response of the model with respect to observations with RMSE/R² values of 0.5 °C/0.86 for T_S and 2 PSU/0.7 for S_S. The weaker correlation obtained for S_S is mainly due to model underestimation around 2 January 2008 (i.e., about 5 PSU). It should be underlined that the model accurately reproduces the abrupt change (decrease) in surface temperature and salinity on 11 December 2007, linked to the veering of the Loire river plume caused by easterly winds. This plume advection led to stratification which in turn influenced the vertical profiles of currents in terms of direction and intensity (Figure 2e,g, respectively), which are well represented by the model. More generally, the model provides an appropriate response regarding the direction and intensity of the current (Figure 2d,f) over the entire water column throughout the period. For instance, measured bottom current velocities (0.11 m·s⁻¹ on average, max of 0.44 m·s⁻¹) are correctly reproduced by the model with a RMSE of 0.05 m·s⁻¹.

2.4. Sediment Transport Model

The sediment model used in the present 3D modelling system is described in Le Hir et al. [1]. The model makes it possible to simulate the transport and changes in all kinds of sediment mixtures under the action of hydrodynamic forcing, including that of waves and currents. According to the model concept, the consolidation state of the bottom sediment is linked to the relative mud concentration, i.e., the concentration of mud in the space between the sand grains (C_{rel mud}). Given that the surficial sediment in our study area is weakly consolidated (erosion actually occurs in surficial layers only, typically a few mm or cm thick) and mainly composed of fine sand, consolidation was disregarded and C_{rel mud} was set at a constant value of 550 kg·m⁻³ (representative of pre-consolidated sediment according to Grasso et al. [41]). In addition, bed load was not taken into account in the present application, assuming that surficial sediment in our study area is mainly composed of mixtures of mud and fine sand. In our case, three sediment classes (for which the mass concentrations are the model state variables) are considered: a fine sand (S1), a non-cohesive material with a representative size of 200 µm, and two muddy classes (M1 and M2), which can be distinguished by their settling velocity, in order to be able to schematically represent the vertical dispersion of cohesive material over the shelf. The sediment dynamics was computed with an advection/dispersion equation for each sediment class, representing transport in the water column, as well as exchanges at the water/sediment interface linked to erosion and deposition processes. Consequently, the concentrations of suspended sediment, the related horizontal and vertical fluxes, and the corresponding changes in the seabed (composition and thickness) are simulated. This section details the way of managing sediment deposition, seabed initialization, and technical aspects linked to vertical discretization within the seabed compartment. The erosion law establishment and the numerical modelling experiment aiming to fit an optimal setting will be addressed independently in Section 3.

2.4.1. Managing Sediment Deposition

The deposition flux D_i for each sediment class i is computed according to the Krone law:

$$D_i=W_{s,i}·SS{C}_i·(1-\frac{\tau }{{\tau }_{d,i}})$$

(3)

In Equation (3), τ is the total bottom shear stress (waves and currents; see Section 2.3.1), W_s,i is the settling velocity, SSC_i is the suspended sediment concentration, and τ_d,i is the critical shear stress for deposition. In the present study, the latter was set to a very high value (1000 N·m⁻²) so it is ineffective: considering that consolidation processes are negligible near the interface and that, as a result, the critical shear stress for erosion remains low (Section 3), deposited sediments can be quickly resuspended in the water column if the hydrodynamic conditions are sufficiently intense, which replaces the role played by the term between parentheses in Equation (3) [1].

The settling velocity of non-cohesive sediment (W_s,S1) is computed according to the formulation of Soulsby [40]. The one related to the mud M1 (W_s,M1) is assumed to vary as a function of its concentration in the water column (SSC_M1) and the ambient turbulence according to the formulation of Van Leussen [42]:

$$W_{s,M1}=max\left\{W_{s,min }; min\left\{W_{s,max} ; k·SS{C}_{M1}{}^m·\left(\frac{1+a·G}{1+b·G^2}\right)\right\}\right\},with G=\sqrt{\frac{\epsilon }{\nu }}$$

(4)

where k, m, a and b are empirical constants, G is the absolute velocity gradient, ε is the turbulent dissipation rate, and ν is the water viscosity. For a and b, respective values of 0.3 and 0.09, set by Van Leussen [42] from experiments, were used. In agreement with the settling velocity setting used by Tessier et al. [43] in a modelling study of turbidity over the southern Brittany continental shelf and recent experiments by Verney et al. [44] in similar environments in the Gulf of Lions, k and m were respectively set to 0.005 and 0.7. Following Tessier et al. [43], W_s,M1 is limited by minimum (W_s,min) and maximum (W_s,max) values, respectively set to 0.1 mm·s⁻¹ and 4 mm·s⁻¹, the latter being reached for SSC ≥ 700 mg·L⁻¹ (thus ignoring the hindered settling process which actually does not occur in the range of SSC over the shelf). Lastly, the mud M2 linked to ambient turbidity in the water column is assumed to have a constant very low settling velocity set to 2.5 × 10⁻⁶ m·s⁻¹.

After resolving the sediment transport equation, the actual deposition on the bottom can be computed. In agreement with Le Hir et al. [1], sand and mud particles are deposited successively, by filling the pores or by creating new layers. A volume concentration of well-sorted sediment (C_{vol sort} = 0.58) was attributed to sediment when only one class of non-cohesive sediment is present. However, the bed concentration can reach a higher value if several classes are mixed (C_{vol mix} = 0.67). These typical volume concentrations [45] led to mass concentrations assuming a fixed sediment density ρs (2600 kg·m⁻³) for all classes. In the case of simultaneous deposition of sand and mud, the sand is first deposited at a concentration which depends on surficial sediment composition: in the case of mixed sediment, the sand is first mixed with the initial mixture until C_mix (= C_{vol mix}.ρ_s = 1742 kg·m⁻³) is reached. The remaining sand is deposited with the concentration C_sort (= C_{vol sort}.ρ_s = 1508 kg·m⁻³), from which an increase in the thickness of the layer can be deduced. The same kind of deposition occurs when the surficial sediment is only comprised of sand. In addition, the thickness of any layer is limited to dz_{sed, max}, a numerical parameter of the model. Any deposition of excess sand leads to the creation of a new layer (see Section 2.4.2). Next, mud is deposited: it progressively fills up pores between the sand grains until either C_mix or C_{rel mud} is reached. Considering these criteria, mud is mixed within the initial and new deposits starting from the water/sediment interface. If any excess mud remains after the mixing step, it is added to the upper sediment layer, contributing to its thickening.

2.4.2. Sediment Discretization within the Seabed

In the model, an initial sediment height (h_{sed, ini}) is introduced, and the seabed is vertically discretized in a given number of layers of equivalent thicknesses (dz_{sed, ini}). An optimal vertical discretization of the sediment was assessed by Mengual [46]. By means of sensitivity analyses, it was shown that beyond a 1/3 mm resolution within the seabed compartment (i.e., thickness of each layer dz_{sed, ini}), the SSC response of the model in the water column did not change anymore. According to the conclusions drawn from this sensitivity analysis, the initial sediment thickness h_{sed, ini} was set at 0.03 m, corresponding to 90 sediment layers of thickness dz_{sed, ini} (1/3 mm).

As previously mentioned in Section 2.4.1, a new layer is created when the actual deposition leads to a thickness of the surficial sediment layer higher than a maximum value dz_{sed, max}, corresponding to a parameter of the model. Nevertheless, a maximum number of layers in the seabed compartment (N_max) needs to be defined in order to make computational costs acceptable. While N_max is reached, a fusion of the two sediment layers located at the base of the sedimentary column occurs. By this way, the creation of new deposited layers becomes once again possible. The parameters dz_{sed, max} and N_max control changes of the “sediment vertical discretization” during the simulation, and are likely to influence the sediment dynamics. To prevent any variations in the SSC response of the model linked to changes in the seabed vertical resolution, the maximum thickness of sediment layers, dz_{sed, max}, was set at the same value than the initial one dz_{sed, ini}.

This seabed management constitutes a compromise between the likely maximum erosion depth in most of the shelf and the possibility of new layers being deposited in other places.

2.4.3. Sediment Facies Initialization for the Application to the BoBCS

The seabed was initialized using the distribution of the three sediment classes according to existing surficial sediment maps (e.g., [47]) and for the sand and mud fractions, new sediment samples taken in the BoBCS (Figure 1). No sediment initialization was undertaken in areas indicated as rocky on sediment maps. The mud M2, characterized by a very low settling velocity, is linked to river inputs and was also uniformly initialized (20 kg·m⁻³) on the seabed over the entire shelf. This particle class enables the representation of a non-negligible part of the ambient turbidity of a few mg·L⁻¹ near the coast.

Given the sand fraction f_s of the mixed sediment, a relationship could be derived between C_{rel mud} (relative mud concentration set at a constant value of 550 kg·m⁻³) and the bulk sediment density (C_bulk). The latter can be written as:

$$C_{bulk}=\frac{C_{rel mud}}{1+f_s·\left(\frac{C_{rel mud}}{{\rho }_s}-1\right)}$$

(5)

For large sand fraction (above 90%), the given relative concentration of 550 kg·m⁻³ becomes unlikely, so that $C_{bulk}$ is maximised by $C_{max}$, which corresponds to a dense sediment (C_sort or C_mix depending on the properties of the sediment mixture, i.e., well-sorted or mixed sediments). This formulation could be validated by comparing with recent data in two sectors of the BoBCS [48]: bulk densities of 1540 kg·m⁻³ and 1380 kg·m⁻³ were obtained for 85% and 75% of sand, respectively, while the application of Equation (5) with sand fractions f_s of 0.85 and 0.75 led to C_bulk values of 1668 kg·m⁻³ and 1346 kg·m⁻³ in fair agreement. Sediment concentrations linked to each particle class were then deduced according to their respective fraction.

Lastly, the sediment distribution has been prescribed uniformly along the vertical dimension (i.e., in each layer of the seabed compartment) in the absence of 3D data on grain size distributions.

3. Erosion Law Setting: Building and Numerical Experiment

3.1. General Formulation

As previously evoked in the introduction, the Partheniades-Ariathurai law [10,11], the most commonly used for cohesive sediments, can also represent the pick up process of fine sands. Such an erosion law is then assumed whatever the sediment composition. In the present study, the erosion flux E is thus expressed following Equation (1). A simultaneous erosion of sandy and muddy fractions is assumed, according to their respective concentration in the mixture.

This erosion law involves three parameters, E₀ (erodibility parameter, kg·m⁻²·s⁻¹), τ_e (critical shear stress, N·m⁻²), and n (hereafter referred to as “erosion-related parameters”), which are set at different values depending on the mud content f_m (<63 µm) of the surficial sediment. The concept of critical mud fraction is retained, with the definition of a first critical fraction f_mcr1 below which a non-cohesive behaviour is prescribed, and of a second one, f_mcr2, above which the sediment is assumed to behave as pure mud. Two sets of erosion-related parameters will be defined below and above these mud fractions (see Section 3.2 and Section 3.3 for sand and mud, respectively). As already mentioned in Section 1, a transition in erosion-related parameters has to be prescribed between f_mcr1 and f_mcr2 (between pure sand and pure mud parameters) to manage the erosion of “transitional” sand/mud mixtures. Such a transition is investigated in Section 3.4.

3.2. Pure Sand Erosion

Below the first critical mud fraction f_mcr1, the non-cohesive erosion regime is prescribed by defining a first set of erosion-related parameters linked to pure sand erosion in Equation (1): E_0,sand, τ_e,sand, and n_sand. Le Hir et al. [1] suggested to compute f_mcr1 as a function of the sand mean diameter D (f_mcr1 = α₀ × D with α₀ = 10³ m⁻¹), leading to a value of 20% considering the fine sand S1 of 200 µm (considered as a reference value hereafter). However, the model erosion dynamics is likely to significantly vary while the surficial mud fraction f_m is close to f_mcr1. Such a sensitivity is addressed in Section 4.2.

In Equation (1), parameters linked to pure sand erosion E_0,sand, τ_e,sand, and n_sand were deduced from numerical simulations and empirical formulations. According to Van Rijn [49] and many other numerical models (e.g., [23,24]), the best fit for n_sand is 1.5. The critical shear stress τ_e,sand was determined from the Shields critical mobility parameter computed according to the formulation of Soulsby [40], leading to a value of 0.15 N·m⁻² for a sand of 200 µm.

Considering that equilibrium conditions are usually met for non-cohesive sediment transport, and that such an equilibrium requires compensation between the erosion rate and the deposition flux, many authors formulate the erosion rate (or pick up function) as:

$$E=W_s·C_{ref}$$

(6)

where W_s is the settling velocity and C_ref is a reference concentration which characterizes the equilibrium. The description of C_ref is generally associated with the reference height, h_ref, the distance from the bed where the concentration is considered. In point of fact, this location has to be the one where the equilibrium between deposition and erosion is considered, with respective values that largely depend on the reference height, because of large concentration gradients near the bed. From the point of view of sediment modelling, this means that the deposition flux at the base of the water column has to be expressed at the exact location where the erosion flux is considered, that is, at the reference height where the equilibrium concentration is given. Van Rijn [49] used the concept of equilibrium concentration as a boundary condition of the computation of suspended sediment profile and fitted the expression:

$$C_{ref}=0.015{\rho }_{s}D/\left(h_{ref}·D_{*}{}^{0.3}\right)·{\left(\tau /{\tau }_e-1\right)}^{1.5}$$

(7)

where $D_{*}$ is the non-dimensional median sand diameter ($D_{*}=D{\left(\sigma g/{\nu }^2\right)}^{1/3}$ with $\sigma =\left({\rho }_s/\rho -1\right)$, ρ being the water density, g the gravitational acceleration, and ν the water viscosity). Using the C_ref formulation of Van Rijn [49] in Equation (6) enables us to express the E_0,sand constant in the “Partheniades” form of the erosion law as:

$$E_{0,sand}=\frac{0.015 {\rho }_s D W_{s,S1}}{h_{ref}{D}_{*}{}^{0.3}}$$

(8)

For D = 200 µm (S1), h_ref = 0.02 m, and W_s,S1 = 2.5 cm·s⁻¹, Equation (8) leads to E_0,sand = 5.94 × 10⁻³ kg·m⁻²·s⁻¹. The relevance of this E_0,sand value was assessed by simulating an equilibrium state under a steady current and by comparing the depth-integrated horizontal sediment flux with some standard transport capacity formulations. Using a 1DV version of the code, several computations of fine sand resuspension were performed under different flow intensities ($Ve{l}_{INT}$), and once the equilibrium was reached (deposition = resuspension) the total transport (Q_sand) was computed as:

$$Q_{sand}=\frac{1}{{\rho }_s}{\int }_{k=1}^{n}Ve{l}_{INT}\left(k\right)·SS{C}_{S1}\left(k\right)·dz\left(k\right)$$

(9)

where dz(k) refers to the thickness of the cell k (in m) in the water column (discretized in n layers along the vertical dimension) and SSC_S1(k) to the suspended sediment concentration of sand (S1) in the cell k. Sand transport rates Q_sand were then compared to the rates deduced from the formulations of Van Rijn [50], Engelund and Hansen [51], and Yang [52] for similar flow velocities (hereafter Q_sand,VR1984, Q_sand,EH1967, and Q_sand,Y1973 respectively). Results are illustrated in Figure 3. The results obtained by Dufois and Le Hir [53], who also used an advection/diffusion model to predict sand transport rates for a wide range of current conditions and numerous sand diameters, have been added in Figure 3 (Q_sand,DLH2015). Figure 3 shows that the sand transport rates obtained from our computations are in a consistent range regarding those obtained with other formulations or studies cited in the literature, demonstrating the suitability of our E_0,sand parameter.

The erosion-related parameters E_0,sand, τ_e,sand, and n_sand for fine sand are summarized in

Table 2. Erosion parameters for pure sand and pure mud sediment.

Erosion Regime	E0 (kg·m−2·s−1)	τe (N·m−2)	n
Non-cohesive (pure sand)	E0,sand = 5.94 × 10−3	τe,sand = 0.15	nsand = 1.5
Cohesive (pure mud)	E0,mud = 10−5	τe,mud = 0.1	nmud = 1

and constitute reference parameters characterizing pure sand erosion.

3.3. Pure Mud Erosion

Above the second critical mud fraction f_mcr2 (reference value of 70% according to the default value used by Le Hir et al. [1]), the cohesive erosion regime is prescribed by defining a second set of erosion-related parameters linked to pure mud erosion in Equation (1): E_0,mud, τ_e,mud, and n_mud.

As frequently specified in the Partheniades-Ariathurai formulation, the n_mud exponent was set to 1. Given the lack of established formulation of the erosion factor for pure mud, experimental approaches are often used to calibrate it for specific materials, preferably in situ when possible. For this purpose, a specific device had been designed: the “erodimeter” is described in Le Hir et al. [7]. It consists of a small recirculating flume where a unidirectional flow is increased step by step and interacts with a sediment sample carefully placed at the bottom after transfer from a cylindrical core. When measurements are made on board an oceanographic vessel, the test can be considered as quasi in situ. On the BoBCS, erosion tests had been conducted on board the Thalassa N/O: a few of them were performed on muddy sediment samples (mud content higher than 80%). Figure 4 illustrates the critical shear stress for erosion τ_e,mud, estimated to be 0.1 N·m⁻², suggesting a barely-consolidated easily erodible sediment.

Concerning the erosion coefficient E_0,mud, the range of values cited in the literature extends from 10⁻³ to 10⁻⁵ kg·m⁻²·s⁻¹ for natural mud beds in open water (e.g., [54]). Simulating fine sediment transport along the BoBCS, Tessier et al. [43] applied the Partheniades erosion law but used an even lower erosion constant (E₀ = 1.3 × 10⁻⁶ kg·m⁻²·s⁻¹). As a first attempt in the present study, a low value E_0,mud = 10⁻⁵ kg·m⁻²·s⁻¹ was used, and its appropriateness was demonstrated by comparing the computed erosion fluxes (E_modelled) from Equation (1) (with τ_e = τ_e,mud and n = n_mud otherwise) with measurements from erodimetry experiments (E_measured) conducted on three muddy sediment samples from the BoBCS (mud fraction higher than 70%) (Figure 5).

The above-mentioned E_0,mud, τ_e,mud, and n_mud values are summarized in Table 2 and constitute reference parameters characterizing pure mud erosion.

3.4. Erosion of Transitional Sand/Mud Mixtures: Selection of Transition Formulations to be Tested

Between critical mud fractions (f_mcr1 and f_mcr2), E₀, τ_e, and n ranged between pure sand (E_0,sand, τ_e,sand, and n_sand) and pure mud (E_0,mud, τ_e,mud, and n_mud) parameters, following a transition trend which had to be specified. We defined several expressions of the erosion law for the transition between non-cohesive and cohesive behaviours as a function of the surficial sediment mud fraction (f_m).

First, we considered a linear transition type in which erosion-related parameters are linearly interpolated from the respective sand and mud parameters according to their respective concentrations in the mixture (default solution in the original paper by Le Hir et al. [1]).

Several experimental studies in the literature revealed that the transition from non-cohesive to cohesive could be more abrupt. For instance, Smith et al. [17] performed erodibility experiments on natural and artificial sand/mud mixtures, and showed a rapid decrease in erosion rates with increasing mud/clay contents. Nevertheless, the transition between the two regimes remains poorly documented. Here, we propose an exponential formulation that specifically enables adjustment of the sharpness of the transition, which is not possible with the few existing expressions (e.g., [1,19,20,23,24]). The exponential transition trend (Equation (10)) was applied to all erosion-related parameters (X_exp) as a function of mud content, with a coefficient C_exp allowing the adjustment of the sharpness of the transition, which becomes more abrupt with an increase in C_exp:

$$X_{exp}=\left(X_{sand}-X_{mud}\right)e^{C_{exp}·P_{exp}}+X_{mud}$$

(10)

where P_exp = (f_mcr1 − f_m)/(f_mcr2 − f_mcr1); X_sand = {E_0,sand; τ_e,sand; n_sand}; and X_mud = {E_0,mud; τ_e,mud; n_mud}.

Different settings of the erosion law were tested in 3D simulations to assess the correctness of the model response in terms of SSC and changes in the seabed. All settings are illustrated in Figure 6. Three kinds of transition, one linear, and two exponentials (with C_exp = 10 and 40 in Equation (10)), were defined to evaluate the effect of the transition trend only, using the reference critical mud fractions (f_mcr1 = 20% and f_mcr2 = 70%). The corresponding settings are named S1_LIN, S1_EXP1, and S1_EXP2, respectively.

A second series of simulations was run to evaluate the sensitivity of the results to critical mud fraction values. Thus, the f_mcr1 value was successively reduced to 10%, 5%, and ~0% (with the corresponding f_mcr2 = 60%, 55%, and 50%), but only the exponential transition regime (with C_exp = 40) was considered (simulations S2_EXP2, S3_EXP2, and S4_EXP2 respectively), as it produced better results than the other transitions (see results Section 4.1).

4. Results

4.1. Influence of the Transition Trend between Non-Cohesive and Cohesive Erosion Modes in the Erosion Law

The first step consisted in assessing the 3D model response in the model cell located closest to the LC station in terms of SSC and changes in the seabed (i.e., composition, thickness) by considering two transitions of the erosion law, one linear and one exponential (with C_exp = 40 in Equation (10)), between the non-cohesive and cohesive regimes (Figure 7). This first comparison was performed using the “reference” f_mcr1 and f_mcr2 values of 20% and 70%, respectively. The two erosion settings, S1_LIN and S1_EXP2, are illustrated in Figure 6 (Section 3.4).

Total bottom shear stresses (i.e., those caused by waves and currents) and barotropic currents in Figure 7a illustrate changes in forcing throughout the period. In the water column, we successively represented SSC over the entire water column (Figure 7b–d) and at 1.67 m above the seabed (i.e., at the level of the first AWAC cell; Figure 7e) for the two simulations and for the AWAC measurements. Changes in the seabed in the two simulations are presented in Figure 7f,g. Lastly, a global sediment budget (in kg) was applied to the model cell used for the comparison, and the contribution of advection (hereafter referred to as F_OBC with F_{OBC, mud} for mud and F_{OBC, sand} for sand, representing the total amount (integrated) of sediment that crosses the borders of the cells along the water column, as net inflow if F_OBC increases or as net outflow if F_OBC decreases) is illustrated in Figure 7h.

Regarding the dynamics of sediments in suspension, SSC from the AWAC profiler ranged between 10 and 80 mg·L⁻¹ over the study period. Four main resuspension events can be identified: the first event lasted from 1 to 5 December 2007 (Event 1), the second from 8 to 11 December 2007 (Event 2), the third from 11 to 12 January 2008 (Event 3), and the last from 15 to 17 January 2008 (Event 4). During these events, SSC values ranged between 20 and 80 mg·L⁻¹ near the seabed, and did not exceed 40 mg·L⁻¹ close to the surface (Figure 7). The rest of the time, a higher frequency turbidity signal linked to the semi-diurnal tide resuspension was recorded near the seabed with SSC in the range of 10–15 mg·L⁻¹. Turbidity peaks were regularly detected in the surface signal whereas there was no significant increase near the seabed: these peaks are probably due to the signal diffraction caused by wave-induced air bubbles (wave mixing; see Tessier [55]). Such a phenomenon also occurs during energetic events (H_s > 2 m) with a SSC signal in the surface higher than in the rest of the water column.

Despite the fact that the model response generally agreed with observations, computing SSC with the S1_LIN erosion setting highlighted some periods during which turbidity was overestimated, for instance in the upper half of the water column during Event 2, and several times between 30 December and 9 January (Figure 7c,e). Overestimations were particularly noticeable in the SSC series at 1.67 m above the seabed (Figure 7e): modelled SSC regularly exceeded observed SSC by 20 to 40 mg·L⁻¹ during Events 2 and 4, and even during calmer periods (e.g., between 30 December and 9 January). In contrast, modelled SSC were underestimated by a factor of 2 during Event 3. The S1_LIN erosion setting led to a representation of observed SSC with a RMSE of 14 mg·L⁻¹ over the study period.

Using the S1_EXP2 erosion setting enabled a general improvement in modelled SSC with a RMSE of 10.5 mg·L⁻¹ over the period, and a correct response during the four energetic periods (Figure 7d,e). Differences in SSC between simulations during Events 1, 2, and 4 were mainly due to different erosion rates (especially due to E₀ in Equation (1)) prescribed for the same intermediate mud fraction in the surficial layer (in Figure 6a, this rate is clearly higher in the S1_LIN setting). However, other differences, especially between 30 of December and 9 of January and during Event 3, are linked to the contrasted changes in the seabed between the two simulations.

One major difference in these changes occurs after Event 2 with more significant mud deposition in simulation S1_EXP2 (Figure 7f,g). This difference can also be seen in Figure 7h, which highlights a significant decrease in F_{OBC, mud} (i.e., flow of mud out of the cell) in simulation S1_LIN, but not in S1_EXP2. The decrease in F_{OBC, mud} in simulation S1_LIN is probably due to less mud inputs from adjacent cells, resulting in relatively more mud exported by advection and thus less mud deposition during the decrease in shear stress following Event 2. Note that this difference in seabed changes influences the sediment dynamics in both simulations throughout the period. Following Event 2, a transition in surficial sediment from muddy to sandy occurs in both simulations but at different times. In the S1_LIN simulation, the transition occurs half way through the period, around 30 December, and manifests itself as a SSC peak linked to mud resuspension near the bottom (Figure 7e), and by a decrease in F_{OBC, mud} (relative loss by advection), which does not occur in the S1_EXP2 simulation. Following the transition in the nature of the seabed, the S1_LIN simulation regularly gives incorrect SSC responses (e.g., overestimation around 6 January, underestimation during Event 3). These results underline the potential role of advection processes in the contrasted results of the two simulations, and the need for full 3D modelling to obtain a final fit of the erosion law. In the S1_EXP2 simulation, the transition occurs later in the period, around 11 of January, and enables a correct SSC response regarding Event 3, associated with a decrease in F_{OBC, mud}. It may mean that, on the one hand, setting S1_EXP2 allows a more accurate representation of resuspension dynamics in response to a given forcing, and on the other hand, it induces a more correct change in the nature of the seabed with respect to the variations in forcing over time. Note that despite contrasted sediment dynamics in the different simulations, the mud budget at the scale of the cell summed over the whole period led to similar trends corresponding to a relative loss by advection of around −5 × 10⁶ kg (that is −5/2.5² kg·m⁻²).

Sand contributes to turbidity over shorter periods than mud, mainly during Events 1, 2, and 4. Results provide evidence that sand is not subject to the same dynamics as mud. Until the transition in the nature of the surficial sediments in the middle of the period (30 December), sand dynamics appear to be quite similar in the two simulations (e.g., F_{OBC, sand} in Figure 7h). Beyond this date, the contrasted nature of the seabed results in more regular sand resuspension in S1_LIN, with an advection component leading to a relative local sand loss (in the cell). Starting from Event 3, F_{OBC, sand} increases (i.e., relative sand inflow into the cell) in both simulations while the advection flux of mud decreases (in S1_EXP2) or does not change (in S1_LIN). This highlights the fact that sand and mud dynamics are likely to differ depending on the nature of the seabed in adjacent cells.

The erosion setting S1_EXP1, which is characterized by a less abrupt exponential transition in erosion law (C_exp = 10 in Equation (10)), appears to be less accurate in terms of SSC (not illustrated here) with a RMSE of 12 mg·L⁻¹ over the study period (versus 10.5 mg·L⁻¹ in S1_EXP2). The turbidity response provided by the model would be expected to be degraded while progressively reducing the decreasing trend of the transition (until a linear decrease is reached).

Considering the accurate representation of resuspension events and the lower RMSE obtained with the erosion setting S1_EXP2, the latter was considered as an “optimum” setting, suggesting that the definition of an exponential transition to describe sand/mud mixture erosion between non-cohesive and cohesive erosion modes may be appropriate in hydro-sedimentary numerical models.

4.2. Influence of Critical Mud Fractions

The sensitivity of the SSC model results to critical mud fractions, f_mcr1 and f_mcr2, was assessed, starting from the “optimum” erosion setting deduced in Section 4.1 and characterized by an exponential transition between f_mcr1 = 20% and f_mcr2 = 70% with C_exp = 40 in Equation (10) (i.e., S1_EXP2 setting). Both critical mud fractions were progressively reduced by 10% (f_mcr1 = 10%, f_mcr2 = 60%), 15% (f_mcr1 = 5%, f_mcr2 = 55%), and 20% (f_mcr1 ≈ 0%, f_mcr2 = 50%). The corresponding settings, S2_EXP2, S3_EXP2, and S4_EXP2 are illustrated in Figure 6. Results linked to the application of these different settings are illustrated in Figure 8. Note that the second critical mud fraction f_mcr2 appears in the extension of the exponential decay (Equation (10)), but, due to the shape of the exponential trend, it does not constitute a real critical mud fraction but rather an adjustment parameter for the transition. We can thus consider that this sensitivity analysis mainly deals with the setting of the first critical mud fraction f_mcr1.

First, it can be seen that the SSC linked to the four resuspension events progressively decrease with the decrease in f_mcr1 (Figure 8b), which results in underestimation of turbidity with respect to observed values. While no clear differences in SSC appear between S1_EXP2 and S2_EXP2 (the latter is not illustrated in Figure 8), i.e., with a reduction of f_mcr1 from 20% to 10%, significant SSC underestimations occur for f_mcr1 < 10%. The average turbidity during resuspension events is underestimated by 15–20% (respectively, 30%) in simulation S3_EXP2 (respectively, S4_EXP2). Regarding maximum SSC, underestimations of SSC peaks are around 15–30% (respectively, 40–50%) during Events 1 and 2 in simulation S3_EXP2 (respectively, S4_EXP2). In addition, SSC peaks during Events 3 and 4 are completely absent in these two simulations with an underestimation of about 60%. Other simulations with linear trend but low f_mcr1 were tested and showed no improvement compared with the settings illustrated in Figure 7 and Figure 8.

Changes in the seabed linked to the optimum erosion setting S1_EXP2 and the simulation S3_EXP2 (f_mcr1 = 5%) are illustrated in Figure 8c,d. Following Event 2, contrary to results in S1_EXP2, no drastic change in the nature of the seabed occurs in S3_EXP2 in the rest of the period, with a surficial sediment containing at least 30–40% of mud. This less dynamic change in the seabed is consistent with the lower SSC obtained in the water column. A reduction of f_mcr1 led to the application of a pure mud erosion law starting from a lower mud content in the surficial sediment. This mostly resulted in less erosion (E_0,mud << E_0,sand) with weaker SSC and slower changes in the seabed. This is also visible in the variations in F_{OBC, mud} (Figure 8e) which highlight the fact that the reduced sediment dynamics obtained by reducing f_mcr1 results in weaker gradients (SSC, seabed nature) with adjacent cells, and a less dynamic advection term over the study period.

5. Discussion

5.1. Setting Describing Erosion of a Sand/Mud Mixture

All experimental studies in the literature on the erosion of sand/mud mixtures mentioned a transition in the erosion mode when fine particles progressively fill the spaces between non-cohesive particles. Panagiotopoulos et al. [15] proposed a conceptual model showing the mechanism for the initiation of sediment motion for sand-mud mixtures, based on the forces acting on an individual grain and the associated angle of internal friction. When the mud content increases, clay minerals progressively fill the spaces between the sand particles, which slightly alters the pivoting characteristics and consequently the internal friction angle, and thus a slight change in erosion resistance. As soon as the sand particles are no longer in contact with one another, pivoting stops being the main mechanism behind the initiation of motion, and the resistance of the clay fraction mainly controls erosion. Depending on the authors, this transition is expressed by reasoning in terms of mud or clay content in the mixture. For instance, Mitchener and Torfs [14] proposed a transition between 3% and 15% mud content, and suggested using cohesive-type sediment transport equations above this transition value and sand transport theories below it. Other authors suggested that the transition occurs at much higher mud content, i.e., between 20% and 30% (e.g., [7,15]). More generally, previous investigations emphasized that only 2% to 10% of clay minerals (dry weight) added in a non-cohesive sediment matrix were sufficient to control the soil properties and increase the resistance of the bed to erosion (e.g., [17,24,45,56]).

Modelling results in the present study highlighted the fact that the critical mud fraction (f_mcr1), above which a transition toward a cohesive erosion mode would start, is at least 10% mud content. Grain size analyses of numerous sediment samples (from several locations, and at different depths in the sediment) from the BoBCS revealed that the clay content (per cent < 4 µm) corresponds to 30% (±3%, R² = 0.96) of the mud content (per cent < 63 µm). Such a constant ratio between the clay and mud fractions (or between the clay and silt fractions) in a given area has been observed in many sites worldwide (e.g., [13]). Thus, the critical clay content deduced from our modelling fitting would be around 3%. Therefore, our results are in agreement with experimental results of previous studies regarding the existence and the value of a critical mud/clay fraction indicating a transition in the mode of erosion.

Multiple transition trends of the erosion law (linear, exponentials) were tested to describe the erosion behaviour of transitional sand/mud mixtures, i.e., when the mud fraction exceeds the critical value of f_mcr1 in the mixture. The quality of the model response was evaluated by comparing SSC results with turbidity measurements provided by the AWAC profiler over the entire water column. Based on RMSE and average or maximum SSC values reached during resuspension events, the results provided a more accurate representation of observations while considering an abrupt exponential transition of erosion parameters (i.e., E₀, τ_e, and n in the Partheniades form of the erosion law, see Equation (1)). Actually, changes in SSC produced by this transition formulation mainly hold in the contrasted E₀ values prescribed in erosion law depending on the seabed mud fraction (E_0,mud << E_0,sand; see Table 2). This result agrees with results recently obtained by Smith et al. [17] who performed erosion experiments on mixed sediment beds prepared in the laboratory (250–500 µm sands mixed with different clayey sediments corresponding to kaolinite or kaolinite/bentonite). In particular, they observed a rapid decrease in erosion rates, from 1.5 to 2.5 orders of magnitude, over a range of 2% to 10% clay content. In the present study, exponential transitions prescribed in settings S1_EXP1 and S1_EXP2 (C_exp = 10 and C_exp = 40 in Equation (10), Figure 6) led to a variation in the erodibility parameter E₀ of about 2.5 orders of magnitude over a mud content range of 10% and 40%, i.e., over a clay content ranging from 3% and 12% respectively. The best model results obtained from erosion setting S1_EXP2 are thus consistent with the findings of Smith et al. [17], and suggest that a rapid exponential transition may be appropriate to describe the erosion of a sand/mud mixture between non-cohesive and cohesive erosion modes in numerical hydro-sedimentary models.

5.2. Limitations of the Approach and Remaining Uncertainties

Despite successive assessments of model quality, some limitations and uncertainties concerning our modelling approach remain and have to be addressed.

5.2.1. Mud Erosion Law

A pure mud erosion law was set up based on erodimetry experiments performed on muddy sediment samples from the BoBCS. A critical shear stress for mud erosion (τ_e,mud) of 0.1 N·m⁻² was deduced. By combining this τ_e,mud value with the minimum erodibility parameter E_0,mud recommended by e.g., Winterwerp [54], i.e., 10⁻⁵ kg·m⁻²·s⁻¹, the application of the mud erosion law from the model (Section 3.3) led to good agreement between modelled erosion fluxes and those obtained in erodimetry experiments for comparable applied shear stresses (Figure 5). Such a lower critical stress for erosion when the mixture is muddier is opposite to trends most often published, characterized by an increase of the resistance to erosion when mud is added to sand (e.g., [6,7,14,15,16,17]). Other simulations were performed with higher τ_e,mud values, 0.15, 0.2, and 0.4 N·m⁻². As expected, modelled SSC was underestimated compared with observed SSC while τ_e,mud increased (even for 0.15 N·m⁻²). Another assessment of E_0,mud could have produced similar results, but we preferred to keep the shear strength provided by our experiments, the low value being justified by the fact that in our environment (erosion on a continental shelf with low bottom friction) the sediment is never remobilized at depth, and the surficial sediment remains unconsolidated.

5.2.2. Initial Condition of the Sediment and Time Variation of the Seabed

Seabed initialization was prescribed from the synthesis of sediment facies applied at the beginning of each simulation. To evaluate the influence of the sediment initialization on model results (SSC, seabed variations), the “optimum” model setting from the present study (S1_EXP2) was used again in a new simulation using the surficial sediment cover computed at the end of a one-year simulation used as spin up. We obtained similar SSC results with a RMSE of 11.3 mg·L⁻¹ over the study period (versus 10.5 mg·L⁻¹ in the original S1_EXP2 simulation). Similarities in seabed variations (thickness and composition) in the two simulations were likewise remarkable (not illustrated here). Thus, the seabed initialization prescribed at the beginning of each simulation appears to be appropriate and does not correspond to a transitional state regarding the sediment dynamics.

Model results concerning changes in the seabed highlighted pronounced gradients in the nature of the surficial sediment in most simulations, with an alternation of muddy and sandy facies depending on the intensity of forcing (e.g., shear stress, advection). Such variations in the nature of the sediment are not unrealistic, since grain size analyses of a sediment core sampled at LC station revealed a layered bed, with alternating muddy and sandy layers at different depths in the sediment. These variations are also consistent with the geographical location of LC station (Figure 1), in a zone with horizontal gradients in sediment facies.

A further validation of the model in terms of thickness or elevation would require other measurements as in situ altimetry data (e.g., [57]). Simulations make it necessary to use altimetry data, which were not available in the present study.

5.2.3. Applicability of the Sand/Mud Mixture Erosion Law

Assessment of the effect of erosion settings on the quality of model results with respect to observations would require further comparisons in other study sites where the seabed consists of both sand and mud. This would make it possible to know if an abrupt transition between non-cohesive and cohesive erosion modes systematically improves model accuracy in terms of SSC.

The use of the sand/mud mixture erosion law derived from this study requires site-specific information beforehand, in particular grain size analyses for the assessment of the mud or the clay content. Since the critical mud fraction f_mcr1, above which erosion behaviour starts to change, mainly depends on clay content, the ratio between clay and mud fractions can be used. In future works, it would be interesting to explore other mud properties than grain size and sediment fractions, such as mineralogy, to represent more accurately the key role played by cohesive sediments in erosion process, especially for transitional sand/mud mixtures between the contrasted non-cohesive/cohesive regimes.

Lastly, the formulation of erosion was based on the Patheniades-Ariathurai law, with an erosion flux proportional to the normalized excess shear stress (τ/τ_e-1). Such a formulation is very sensitive to the value of the critical shear stress for erosion, which can be difficult to estimate and highly variable in the case of sand/mud mixtures. Alternatively, a formulation of the erosion flux proportional to the excess shear stress (τ-τ_e) would reduce the sensitivity of erosion to τ_e. It would also be in agreement with the Van Kesteren-Winterwerp-Jacobs erosion law [6,58], and deserves further investigations following the pioneering work of Jacobs et al. [6].

6. Conclusions

The aim of this study was to assess the influence of the erosion law prescribed in a 3D realistic hydro-sedimentary model on sediment dynamics in the case of a seabed composed of both fine sand and mud in a slightly energetic environment, representative of continental shelves. According to the sediment model described by Le Hir et al. [1], the sediment was eroded as a mixture and was assumed to behave as pure sand below a first critical mud fraction in the surficial sediment, and as pure mud above a second one. Following hydrodynamic validation of the model and rigorous assessment of pure sand and pure mud erosion dynamics, several transition trends of erosion-related parameters (erodibility parameter, critical shear stress, and exponent in the Partheniades erosion law) were tested to describe the erosion of transitional sand/mud mixtures between the two critical fractions. Different simulations were run using linear or exponential transitions, and different critical mud fractions. The accuracy of model results regarding suspended matter dynamics was evaluated at a single point, located on the Bay of Biscay shelf, by performing comparisons with turbidity observations provided by an acoustic profiler during two typical winter months. The main conclusions of this work are:

• Using an abrupt exponential transition, e.g., an erodibility parameter decrease of 2.5 orders of magnitude over a 10% (respectively, 3%) mud (respectively, clay) content range, improves SSC model results regarding measurements, compared to results obtained with linear or less abrupt exponential transitions. This conclusion agrees with recent experimental studies in the literature on the erosion of sand/mud mixtures, which mention a drastic change in erosion mode for only a small percentage of clay added in the mixture.
• A first critical mud fraction (above which the erosion mode begins to change) of 10–20% is required to ensure a relevant model response in turbidity. By reasoning in terms of the clay fraction, the corresponding critical clay fraction ranges between 3% and 6%. Once again, this conclusion agrees with experimental studies in the literature reporting that 2% to 10% of clay minerals in a sediment mixture are sufficient to control the soil properties.
• The erosion flux of mixed sediments appears to be very sensitive to the clay fraction of the surficial sediment, and then is likely to change considerably at a given location, according to erosion and deposition events.
• The need to perform 3D simulations to account for advection, which considerably influences sediment dynamics in terms of export of resuspended sediments, sediment inflows from adjacent cells, and consequent changes in the surficial seabed (nature and thickness of deposits).

Therefore, the optimal erosion law derived from this study to describe sand/mud mixture erosion led to model results consistent with measurements and with most of the conclusions deduced from experimental studies already published. This should encourage further similar comparisons and suggests that the application of this kind of erosion setting is appropriate for hydro-sedimentary models.