# Erosion Processes¶

Resuspension Erosion flux E: (Partheniades)\(E=E_0*\alpha*[\frac{\tau-\tau_{ce}}{\tau_{ce}}]^n\)

- \(E_0\) : erodability
- \(\alpha\) : correction factor
- \(\tau\) : shear stress
- \(\tau_{ce}\) : critical shear stress for erosion
- \(n\) : power coefficient applied to excess erosion

- The erosion flux is differently treated in
**version 1**(Erosion Processes : Version 1) and**version 2**(Erosion and Bedload Processes : Version 2) where erosion flux by bedload are added for sand and gravels - A lateral erosion of a wet or dry cell was introduced into the model (only for suspensions) (Lateral erosion)

## Erosion Processes : **Version 1**¶

Erosion flux is applied in proportion to the classes present, by modulating parameters depending on the nature of the superficial sediment

the

**critical erosion shear stress**of erosion \(\tau_{ce}\) depends on the content of sand / mud and on sediment consolidation.- \(\alpha\) is a correction factor
- \(\alpha=1.+corfluer1*\max(0,corfluer2-Csed_{total-surf})\) (
*corfluer1*et*corfluer2*given in*paraMUSTANGV1.txt*)

- \(\alpha=1.+corfluer1*\max(0,corfluer2-Csed_{total-surf})\) (

- \(E_0\) and \(n\) depend on the mud fraction in sediment
- \(fr_{mud} < frv1\) : sandy sediment (frv1=0.3 typically =
*frmudcr1*calculated from*coef_frmudcr1*in*paraMUSTANGV1.txt*) - \(fr_{mud} > frv2\) : muddy sediment (frv2=0.7 typically =
*frmudcr2*given in*paraMUSTANGV1.txt*)

- \(fr_{mud} < frv1\) : sandy sediment (frv1=0.3 typically =

- For muddy sediment :
- \(E_0=E0_{mud}\) (given in
*paraMUSTANGV1.txt*) - \(\tau_{ce}=x1toce_{mud}*Crel_{mud}^{x2toce_{mud}}\) (parameters given in
*paraMUSTANGV1.txt*) with \(Crel_{mud}\) which is the relative mud concentration between sand grains

- \(E_0=E0_{mud}\) (given in

- For sandy sediment :
\(E_0=E0_{sand}*E0_{sand-para}\)

- \(E0_{sand}\) given in
*paraMUSTANGV1.txt*or estimated by formulations depending on the option chosen by**E0_sand_option** *E0_sand_option*= 0 ==> E0_sand= E0_sand_cst read in this namelist*paraMUSTANGV1.txt**E0_sand_option*= 1 ==> E0_sand evaluated with Van Rijn (1984) formulation*E0_sand_option*= 2 ==> E0_sand evaluated with erodimetry formulation\(E0_{sand}=\min(0.27,1000*d50-0.01)*(\tau-\tau_{ce})^{n_{eros-sand}}\)

*E0_sand_para*is a multiplicative factor to correct erosion flux

- \(E0_{sand}\) given in
\(\tau_{ce}\) estimated from sand characteristics

For sand/mud mixing : Erosion flux depends on proportions of the mixture

- The option choice is given by
**ero_option**in*paraMUSTANGV1.txt* *ero_option*= 0 : pure mud behavior (for all particles and whatever the mixture)*ero_option*= 1 : linear interpolation between sand and mud behavior, depend on proportions of the mixture*ero_option*= 2 : formulation derived from that of J. Vareilles (2013)*ero_option*= 3 : formulations proposed by B. Mengual (2015) with exponential coefficients depend on proportions of the mixture

- The option choice is given by

Note

Several options are available but all are questionable and should be tested and used carefully

## Erosion and Bedload Processes : **Version 2**¶

The new model allows a sharper distinction of erosion regimes, non-cohesive and cohesive, by offering the possibility of adding the transport process by bedload for gravel and sand and a flux of resuspension independent of the different classes in non-cohesive regime.

Just like the version 1, the transition from one regime to another depends on the nature of the sediment through the concept of a critical mud fraction

*frmudcr1*(calculated as in Le Hir et al (2011) as a function of average diameter of the sands). From the composition of the surface layer and the*frmudcr1*associated, it is possible to know if the sediment in the superficial layer is non-cohesive or cohesive.During a time step, different iterations can occur by alternating between the two regime. Whatever the number of iterations, the flows by bedload and resuspension are gradually integrated during the time step for each class of sediment and in each mesh.

- Basic principles for calculating
**bedload**: - The horizontal fluxes (kg/m/s) are projected parallel to the shear stress and evaluated at the center of the meshes.
- The two components of the flux are taken entirely during time step, and treated as erosion. The bedload contributions of the adjacent meshes are then treated as an additional deposit.
- The flux taken is limited according to the quantity available in the active layer
- Each component is assigned to the “downstream” mesh edge: the flow is “decentered upstream”

- Basic principles for calculating
**Initialization of erosion parameters for non cohesive sediments (gravel and sand)**:initialization of the erosion parameters is done at the beginning of the simulation according to the options chosen in

*paraMUSTANGV2.txt*- \(E_0\) depending on the boolean
***E0_sand_option*** *E0_sand_option*= 0 ==> E0_sand= E0_sand_cst read in this namelist*paraMUSTANGV2.txt**E0_sand_option*= 1 ==> E0_sand evaluated with Van Rijn (1984) formulation, modified by J. Vareilles (2013)*E0_sand_option*= 2 ==> E0_sand evaluated with erodimetry formulation

\(E0_{sand}=\min(0.27,1000*d50-0.01)*(\tau-\tau_{ce})^{n_{eros-sand}}\)

*E0_sand_option*= 3 ==> formulation of Wu and Lin (2014): E0 is deduced from the equilibrium relation between the flux of erosion and the product of the settling rate by the reference concentration C_{ref} at a given height relative to the bottom (with a C_{ref} formulation specific to Wu and Lin, 2014).- E0 is then multiplicated by
*E0_sand_para*(corrective factor)

- \(E_0\) depending on the boolean
- \(\tau_{ce}\) depending on the boolean
***tau_cri_option*** *tau_cri_option*= 0 ==> \(\tau_{ce}\) calculated from the Shields formula, modified by Soulsby (1997) for small diameters.*tau_cri_option*= 1 ==> \(\tau_{ce}\) calculated according to the formulation of Wu and Lin (2014) ignoring the factor associated with masking / exposure effects. This formulation has the particularity of using a Shields mobility criterion constant of 0.03. The effects of masking / exposure can be added thanks to the activation of Booleans to inform in*paraMUSTANGV2.txt*.

- \(\tau_{ce}\) depending on the boolean
\(n=n_{eros-sand}\) given in

*paraMUSTANGV2.txt*, applicated to shear stress excess for all non-cohesives classes.

**Estimation of erosion flux**:the surface sediment is characterized by its mud fraction which is compared with

*frmudcr1*calculated as a function of average diameter of the sandsif

*frmud*> frmudcr1 : the sediment is treated in cohesif regime, and non cohesif if not.- There are two ways to deal with erosion of non-cohesive sediment, and the choice is made with the boolean
**l_eroindep_noncoh**: - if
**l_eroindep_noncoh = False**: - the erosion is treated very similar to version V1, without integrating bedload.
- All classes erode homogeneously and the erosion parameters are modulated according to the nature of the sediment.
- The procedure is the same as in a cohesive regime.

- if
- if
**l_eroindep_noncoh = True**: - the erosion is treated independently for different classes of sand and silt;
- bedload is taken into account

- if

- There are two ways to deal with erosion of non-cohesive sediment, and the choice is made with the boolean
workflow to manage erosion fluxes and reworking of sediment layers is shown in erosion routine diagram - Version V2

**non-cohesive regime**, i.e. the superficial mud fraction is less than*frmudcr1*. erosion is treated class by class independently (if*l_eroindep_noncoh*= True*)- the different classes of sand and mud are independant (the gravel being transported only by bedload).
- If the critical stress representative of the sedimentary mixture in the surface layer is greater than the shear stress near the bottom (waves + current), then an
**active layer**is created. - The thickness of the active layer is calculated from the formulation of Harris and Wiberg (1997). A theoretical active layer thickness is calculated based on the excess shear near the bottom and the representative diameter of the whole non-cohesive classes (gravel and sand). We use two coefficients
*k1HW97*and*k2HW97*, whose value is given in*paraMUSTANGV2.txt*(default 0.07 and 6, respectively, according to Harris and Wiberg, 1997). - thanks to an iterative procedure, the layers of sediment from the water / sediment interface are merged until one of the following conditions is validated:
- the thickness of fused sediment reaches the theoretical thickness of the active layer,

- a cohesive layer is encountered (
*fm*>*frmudcr1*),

- a cohesive layer is encountered (
- there is only one layer left.

- this iterative process is modulated by using a parameter
*fusion_para_activlayer*(given in*paraMUSTANGV2.txt*) : - if
*fusion_para_activlayer*=0 ==> no fusion : superficial layer = active layer - if
*fusion_para_activlayer*=1 ==> criterion (ii) used (*fm*>*frmudcr1*) - if 0 <
*fusion_para_activlayer*< frmudcr2 ==> criterion (ii) modulated (\(fm > fusion_{para_activlayer}*frmudcr1\) )

- if

- this iterative process is modulated by using a parameter

- The thickness of the active layer is calculated from the formulation of Harris and Wiberg (1997). A theoretical active layer thickness is calculated based on the excess shear near the bottom and the representative diameter of the whole non-cohesive classes (gravel and sand). We use two coefficients
- If
**key_MUSTANG_bedload**is used,**bedload fluxes**are evaluated for all non cohesive classes whode boolean*l_bedload*is equal to .TRUE. (given in variable.dat). These developments come from A. Rivier et al (2017) - The formulation chosen for calculating the fluxes carried is that of Wu and Lin (2014). This calculation depends in particular excess shear relative to the critical erosion stress.
- By activating the Boolean
*l_peph_bedload*in*paraMUSTANGV2.txt*, the critical stress of erosion relative to a given class will account of a masking / exposure effect related to the presence of other non-cohesive classes (see Wu and Lin, 2014). - These fluxes (in kg/m/s) are then projected in (x, y) according to the direction of the bottom shear stress exerted by the current.
- The
**slope effects**likely to modulate the bedload fluxes can be taken into account by informing*l_slope_effect_bedload*to .TRUE. in*paraMUSTANGV2.txt*, (from Lesser’s theory, 2004). - Finally, the fluxes carried in (x, y) are set to 0 if the neighboring mesh is exonde.

- If
**resupension fluxes**are then evaluated for each sediment variable independently (Partheniades-Ariathurai) with different parameters (\(E_0\), \(\tau_{ce}\) and \(n\) ) specific to each variable.- For
**non-cohesive classes**, the activation of a boolean*l_peph_suspension*(given in*paraMUSTANGV2.txt*) allows to take into account the effects of**masking / exposure**related to the presence of other noncohesive classes (see Wu and Lin, 2014) on the critical erosion stress estimation relative to a given class. - Still in the case of
**non-cohesive classes**if the way of calculating the parameter of erodibility \(E_0\) does not follow the formalism of Wu and Lin (2014) (i.e.,*E0_sand_option*≠ 3) and that class is also carried by bedload, so it is possible to correct the resuspension flow of the class to satisfy the total transport and the bedload / resuspension ratios described by Wu and Lin (2014) in response to a constraint forcing. This is possible thanks to the activation of the boolean*l_fsusp*in the*paraMUSTANGV2.txt*, which leads to the application of a factor to the \(E_0\) of the class concerned (*fsusp*) equal to the ratio of transport in suspension and total transport associated with the constraint regime. - For
**cohesive classes**, only transported in suspension, different options have been introduced concerning \(E_0\) and \(\tau_{ce}\). - for \(E_0\) , it is possible to apply a factor (
*E0_mud_para_noncoh*) to the parameter erosion of the*E0_mud*(defined in*paraMUSTANGV2.txt*) to make the mud more or less erodible in a non-cohesive regime. - for \(\tau_{ce}\) , there are several options with
***tau_cri_mud_option_eroindep***(given in*paraMUSTANGV2.txt*) *tau_cri_mud_option_eroindep*= 0 ==> evaluated from the relative mud concentration betwwen non cohesives grains*tau_cri_mud_option_eroindep*= 1 ==> minimum value between the average critical stress representative of sands within the mixture and the critical stress calculated from the relative mud concentration in the non-cohesive matrix.*tau_cri_mud_option_eroindep*= 2 ==> minimum value of products between critical constraints and fractions for all sands present within the mixture.*tau_cri_mud_option_eroindep*= 3 ==> minimum value between the critical stress of the finest sand and that calculated from the relative mud concentration between the non-cohesive grains.

- for \(\tau_{ce}\) , there are several options with

- for \(E_0\) , it is possible to apply a factor (

- For

- For

**sediment erosion in the active layer**. The total erosion flux including both resuspension and bedload is calculated class by class. This flux multiplied by the time allotted to erosion is compared to the sediment mass available in the active layer.- If the mass to erode is greater than or equal to the mass of the available sediment class, then the active layer is completely emptied of the concerned class.
- If not, only part of the concerned sediment class contained in the active layer is eroded
- The update of the characteristics
*post-erosion*of the surface layer (thickness, concentrations) is done via a new porosity calculation taking into account the different classes remaining. - If the new layer thickness is less than
*dzsmin*, then the layer is merged with that from below to prevent cases where one would obtain infinitely small concentrations. - The resuspension and bedload fluxes related to the iteration are added to those from previous iterations made during the same time step (including iterations in non-cohesive and cohesive regimes).

**cohesive regime**In cohesive regime, i.e., when the superficial mud fraction (

*fm*) is greater than*frmudcr1*, the physical behavior of the matrix gravel-sandy is influenced by muddy sediments. This influence grows as the mud fraction increases within the mixture. In this case the erosion of the mixture applies homogeneous.- \(E_0\), \(\tau_{ce}\) and \(n=n\) are calculated according to the respective fractions of the different classes within the mixture. The erosion parameters in a cohesive regime tend (more or less rapidly) towards those associated with a purely muddy behavior (when fm increases). Different options have been introduced to manage the type and speed of transition into a cohesive regime. Depending on the value of the
**ero_option**(to be informed in*paraMUSTANGV2.txt*), the parameters of erosion in cohesive regime are computed in different ways *ero_option*= 0
: pure mud behavior- \(E_0=E0_{mud}\) (given in
*paraMUSTANGV2.txt*) - \(\tau_{ce}=x1toce_{mud}*Crel_{mud}^{x2toce_{mud}}\) (parameters given in
*paraMUSTANGV2.txt*) with \(Crel_{mud}\) which is the relative mud concentration between non cohesives grains - \(n=n_{eros-mud}\) given in
*paraMUSTANGV2.txt*, applicated to shear stress excess for all mud variables.

- \(E_0=E0_{mud}\) (given in

*ero_option*= 1 : linear interpolation between sand and mud behavior, depend on proportions of the mixture. The transition regime applies between*frmudcr1*and a second critical fraction*frmudcr2*, to be informed in the*paraMUSTANGV2.txt*(default: 0.5).*ero_option*= 2 : same logic as case 1 except that the transition follows a power law (e.g. Carniello et al, 2012)*ero_option*= 3 : same logic as case 1 except that the transition varies exponentially and more or less fast (formulations proposed by B. Mengual (2017) and modified in 2018). The coefficient controlling the speed of the exponential transition is either equal to that prescribed by the user via the parameter*xexp_ero*(in*paraMUSTANGV2.txt*; default: 10), either, if*l_xexp_ero_cst*=.FALSE. , increases linearly depending on*frmudcr1*. So the more the average diameter of the sands increases, the more the transition between non-cohesive and cohesive regime intervenes for a large fraction of mud but the faster this transition is.

- \(E_0\), \(\tau_{ce}\) and \(n=n\) are calculated according to the respective fractions of the different classes within the mixture. The erosion parameters in a cohesive regime tend (more or less rapidly) towards those associated with a purely muddy behavior (when fm increases). Different options have been introduced to manage the type and speed of transition into a cohesive regime. Depending on the value of the
if \(\tau > \tau_{ce}\) , erosion flux is then evaluated (Partheniades)

the effective erosion of the first layer of sediment is carried out. In the case where quantity to be eroded is smaller than or equal to the sandy-muddy stock in the superficial layer, then the latter is eroded in proportion to the classes. The “post-erosion” characteristics (thickness, concentrations) of the layer are updated ensuring that its thickness remains strictly greater than

*dzsmin*. If this condition is not met, then the remaining quantities of sediment are placed in the layer underneath.- If the amount of sandy-mud sediment in the superficial layer is not enough to satisfy erosion, then all classes of sands and muds are eroded.
- If this layer is devoid of gravel then the erosion of the latter is total.
- If gravel is present, a new porosity calculation makes it possible to update the thickness of the layer and a paving effect is put in place.

If erosion time has not been fully consumed and if the sediment has more than one layer, then a new iteration will be initiated, but this time in a non-cohesive regime.

## Lateral erosion¶

- A lateral erosion process of a wet or dry cell was introduced into the model (only for suspensions).
- Parameters
*coef_erolat*,*l_erolat_wet_cell*and*coef_tenfon_lat*are used for this option., given in*paraMUSTANGV1.txt*or*paraMUSTANGV2.txt*

- Parameters

**dry cell**:\(Ero_{lat}=(\alpha_{lat} * U_{neighb}^2-\tau_{ce})*H_{neighb}\)**wet cell**(if*l_erolat_wet_cell = .TRUE.*)\(Ero_{lat}=(\beta_{lat} * \alpha_{lat} * U_{neighb}^2-\tau_{ce})*\delta_H\)with \(\alpha_{lat}\) =coef_tenfon_latwith \(\beta_{lat}\) =coef_erolatwith \(U_{neighb}\) =water current in the neigboring cellwith \(H_{neighb}\) =depth in the neighboring cellwith \(\delta_H\) =depth difference between eroded cell and the neighboring