# Model description¶

The model approach of [dVvTdVvR+14] is extended to compute the spatiotemporal varying sediment availability through simulation of the process of beach armoring. For this purpose the bed is discretized in horizontal grid cells and in vertical bed layers (2DV). Moreover, the grain size distribution is discretized into fractions. This allows the grain size distribition to vary both horizontally and vertically. A bed composition module is used to compute the sediment availability for each sediment fraction individually. This model approach is a generalization of existing model concepts, like the shear velocity threshold and critical fetch, and therefore compatible with these existing concepts.

A 1D advection scheme is adopted in correspondence with [dVvTdVvR+14] in which $$c$$ [$$\mathrm{kg/m^2}$$] is the instantaneous sediment mass per unit area in transport:

(1)$\frac{\partial c}{\partial t} + u_z \frac{\partial c}{\partial x} = E - D$

$$t$$ [s] denotes time and $$x$$ [m] denotes the cross-shore distance from a zero-transport boundary. $$E$$ and $$D$$ [$$\mathrm{kg/m^2/s}$$] represent the erosion and deposition terms and hence combined represent the net entrainment of sediment. Note that Equation (1) differs from Equation 9 in [dVvTdVvR+14] as they use the saltation height $$h$$ [m] and the sediment concentration $$C_{\mathrm{c}}$$ [$$\mathrm{kg/m^3}$$]. As $$h$$ is not solved for, the presented model computes the sediment mass per unit area $$c = h C_{\mathrm{c}}$$ rather than the sediment concentration $$C_{\mathrm{c}}$$. For conciseness we still refer to $$c$$ as the sediment concentration.

The net entrainment is determined based on a balance between the equilibrium or saturated sediment concentration $$c_{\mathrm{sat}}$$ [$$\mathrm{kg/m^2}$$] and the instantaneous sediment transport concentration $$c$$ and is maximized by the available sediment in the bed $$m_{\mathrm{a}}$$ [$$\mathrm{kg/m^2}$$] according to:

(2)$E - D = \min \left ( \frac{\partial m_{\mathrm{a}}}{\partial t} \quad ; \quad \frac{c_{\mathrm{sat}} - c}{T} \right )$

$$T$$ [s] represents an adaptation time scale that is assumed to be equal for both erosion and deposition. A time scale of 1 second is commonly used ([dVvTdVvR+14]).

The saturated sediment concentration $$c_{\mathrm{sat}}$$ is computed using an empirical sediment transport formulation (e.g. [Bag37b]):

(3)$q_{\mathrm{sat}} = \alpha C \frac{\rho_{\mathrm{a}}}{g} \sqrt{\frac{d_{\mathrm{n}}}{D_{\mathrm{n}}}} \left ( u_z - u_{\mathrm{th}} \right )^3$

in which $$q_{\mathrm{sat}}$$ [kg/m/s] is the equilibrium or saturated sediment transport rate and represents the sediment transport capacity. $$u_z$$ [m/s] is the wind velocity at height $$z$$ [m] and $$u_{\mathrm{th}}$$ the velocity threshold [m/s]. The properties of the sediment in transport are represented by a series of parameters: $$C$$ [–] is a parameter to account for the grain size distribution width, $$\rho_{\mathrm{a}}$$ [$$\mathrm{kg/m^3}$$] is the density of the air, $$g$$ [$$\mathrm{m/s^2}$$] is the gravitational constant, $$d_{\mathrm{n}}$$ [m] is the nominal grain size and $$D_{\mathrm{n}}$$ [m] is a reference grain size. $$\alpha$$ is a constant to account for the conversion of the measured wind velocity to the near-bed shear velocity following Prandtl-Von Kármán’s Law of the Wall: $$\left(\frac{\kappa}{\ln z / z'} \right)^3$$ in which $$z'$$ [m] is the height at which the idealized velocity profile reaches zero and $$\kappa$$ [-] is the Von Kármán constant.

The equilibrium sediment transport rate $$q_{\mathrm{sat}}$$ is divided by the wind velocity $$u_z$$ to obtain a mass per unit area (per unit width):

$c_{\mathrm{sat}} = \max \left ( 0 \quad ; \quad \alpha C \frac{\rho_{\mathrm{a}}}{g} \sqrt{\frac{d_{n}}{D_{n}}} \frac{\left ( u_z - u_{\mathrm{th}} \right )^3}{u_z} \right )$

in which $$C$$ [–] is an empirical constant to account for the grain size distribution width, $$\rho_{\mathrm{a}}$$ [$$\mathrm{kg/m^3}$$] is the air density, $$g$$ [$$\mathrm{m/s^2}$$] is the gravitational constant, $$d_{\mathrm{n}}$$ [m] is the nominal grain size, $$D_{\mathrm{n}}$$ [m] is a reference grain size, $$u_z$$ [m/s] is the wind velocity at height $$z$$ [m] and $$\alpha$$ [–] is a constant to convert from measured wind velocity to shear velocity.

Note that at this stage the spatial variations in wind velocity are not solved for and hence no morphological feedback is included in the simulation. The model is initially intended to provide accurate sediment fluxes from the beach to the dunes rather than to simulate subsequent dune formation.

## Multi-fraction Erosion and Deposition¶

The formulation for the equilibrium or saturated sediment concentration $$c_{\mathrm{sat}}$$ (Equation equilibrium-transport) is capable of dealing with variations in grain size through the variables $$u_{\mathrm{th}}$$, $$d_{\mathrm{n}}$$ and $$C$$ ([Bag37b]). However, the transport formulation only describes the saturated sediment concentration assuming a fixed grain size distribution, but does not define how multiple fractions coexist in transport. If the saturated sediment concentration formulation would be applied to each fraction separately and summed up to a total transport, the total sediment transport would increase with the number of sediment fractions. Since this is unrealistic behavior the saturated sediment concentration $$c_{\mathrm{sat}}$$ for the different fractions should be weighted in order to obtain a realistic total sediment transport. Equation (2) therefore is modified to include a weighting factor $$\hat{w}_k$$ in which $$k$$ represents the sediment fraction index:

(4)$E_k - D_k = \min \left ( \frac{\partial m_{\mathrm{a},k}}{\partial t} \quad ; \quad \frac{\hat{w}_k \cdot c_{\mathrm{sat},k} - c_k}{T} \right )$

It is common to use the grain size distribution in the bed as weighting factor for the saturated sediment concentration (e.g. [Delft3DFManual14], section 11.6.4). Using the grain size distribution at the bed surface as a weighting factor assumes, in case of erosion, that all sediment at the bed surface is equally exposed to the wind.

Using the grain size distribution at the bed surface as weighting factor in case of deposition would lead to the behavior where deposition becomes dependent on the bed composition. Alternatively, in case of deposition, the saturated sediment concentration can be weighted based on the grain size distribution in the air. Due to the nature of saltation, in which continuous interaction with the bed forms the saltation cascade, both the grain size distribution in the bed and in the air are likely to contribute to the interaction between sediment fractions. The ratio between both contributions in the model is determined by a bed interaction parameter $$\zeta$$.

The weighting of erosion and deposition of individual fractions is computed according to:

\begin{split}\begin{align} \hat{w}_k &= \frac{w_k}{ \sum_{k=1}^{n_{\mathrm{k}}}{w_k} } \\ \mathrm{where} \quad w_k &= (1 - \zeta) \cdot w^{\mathrm{air}}_k + (1 - \hat{S}_k) \cdot w^{\mathrm{bed}}_k \end{align}\end{split}

in which $$k$$ represents the sediment fraction index, $$n_{\mathrm{k}}$$ the total number of sediment fractions, $$w_k$$ is the unnormalized weighting factor for fraction $$k$$, $$\hat{w}_k$$ is its normalized counterpart, $$w^{\mathrm{air}}_k$$ and $$w^{\mathrm{bed}}_k$$ are the weighting factors based on the grain size distribution in the air and bed respectively and $$\hat{S}_k$$ is the effective sediment saturation of the air. The weighting factors based on the grain size distribution in the air and the bed are computed using mass ratios:

(5)$w^{\mathrm{air}}_k = \frac{c_k}{c_{\mathrm{sat},k}} \quad ; \quad w^{\mathrm{bed}}_k = \frac{m_{\mathrm{a},k}}{\sum_{k=1}^{n_{\mathrm{k}}}{m_{\mathrm{a},k}}}$

The sum of the ratio $$w^{\mathrm{air}}_k$$ over the fractions denotes the degree of saturation of the air column for fraction $$k$$. The degree of saturation determines if erosion of a fraction may occur. Also in saturated situations erosion of a sediment fraction can occur due to an exchange of momentum between sediment fractions, which is represented by the bed interaction parameter $$\zeta$$. The effective degree of saturation is therefore also influenced by the bed interaction parameter and defined as:

$\hat{S}_k = \min \left ( 1 \quad ; \quad (1 - \zeta) \cdot \sum_{k=1}^{n_{\mathrm{k}}} w_k^{\mathrm{air}} \right )$

When the effective saturation is greater than or equal to unity the air is (over)saturated and no erosion will occur. The grain size distribution in the bed is consequently less relevant and the second term in Equation (5) is thus minimized and zero in case $$\zeta = 0$$. In case the effective saturation is less than unity erosion may occur and the grain size distribution of the bed also contributes to the weighting over the sediment fractions. The weighting factors for erosion are then composed from both the grain size distribution in the air and the grain size distribution at the bed surface. Finally, the resulting weighting factors are normalized to sum to unity over all fractions ($$\hat{w}_k$$).

The composition of weighting factors for erosion is based on the saturation of the air column. The non-saturated fraction determines the potential erosion of the bed. Therefore the non-saturated fraction can be used to scale the grain size distribution in the bed in order to combine it with the grain size distribution in the air according to Equation (5). The non-saturated fraction of the air column that can be used for scaling is therefore $$1 - \hat{S}_k$$.

For example, if bed interaction is disabled ($$\zeta = 0$$) and the air is 70% saturated, then the grain size distribution in the air contributes 70% to the weighting factors for erosion, while the grain size distribution in the bed contributes the other 30% (Figure Fig. 1, upper left panel). In case of (over)saturation the grain size distribution in transport contributes 100% to the weighting factors and the grain size distribution in the bed is of no influence. Transport progresses in downwind direction without interaction with the bed.

Fig. 1 Contributions of the grain size distribution in the bed and in the air to the weighting factors $$\hat{w}_k$$ for the equilibrium sediment concentration in Equation (4) for different values of the bed interaction parameter.

To allow for bed interaction in saturated situations in which no net erosion can occur, the bed interaction parameter $$\zeta$$ is used (Figure Fig. 1). The bed interaction parameter can take values between 0.0 and 1.0 in which the weighting factors for the equilibrium or saturated sediment concentration in an (over)saturated situation are fully determined by the grain size distribution in the bed or in the air respectively. A bed interaction value of 0.2 represents the situation in which the grain size distribution at the bed surface contributes 20% to the weighting of the saturated sediment concentration over the fractions. In the example situation where the air is 70% saturated such value for the bed interaction parameter would lead to weighting factors that are constituted for $$70\% \cdot (100\% - 20\%) = 56\%$$ based on the grain size distribution in the air and for the other 44% based on the grain size distribution at the bed surface (Figure Fig. 1, upper right panel).

The parameterization of the exchange of momentum between sediment fractions is an aspect of saltation that is still poorly understood. Therefore calibration of the bed interaction parameter $$\zeta$$ is necessary. The model parameters in Equation equilibrium-transport can be chosen in accordance with the assumptions underlying multi-fraction sediment transport. $$C$$ should be set to 1.5 as each individual sediment fraction is well-sorted, $$d_{\mathrm{n}}$$ should be chosen equal to $$D_{\mathrm{n}}$$ as the grain size dependency is implemented through $$u_{\mathrm{th}}$$. $$u_{\mathrm{th}}$$ typically varies between 1 and 6 m/s for sand.

## Simulation of Sediment Sorting and Beach Armoring¶

Since the equilibrium or saturated sediment concentration $$c_{\mathrm{sat},k}$$ is weighted over multiple sediment fractions in the extended advection model, also the instantaneous sediment concentration $$c_k$$ is computed for each sediment fraction individually. Consequently, grain size distributions may vary over the model domain and in time. These variations are thereby not limited to the horizontal, but may also vary over the vertical since fine sediment may be deposited on top of coarse sediment or, reversely, fines may be eroded from the bed surface leaving coarse sediment to reside on top of the original mixed sediment. In order to allow the model to simulate the processes of sediment sorting and beach armoring the bed is discretized in horizontal grid cells and vertical bed layers (2DV; Figure Fig. 2).

The discretization of the bed consists of a minimum of three vertical bed layers with a constant thickness and an unlimited number of horizontal grid cells. The top layer is the bed surface layer and is the only layer that interacts with the wind and hence determines the spatiotemporal varying sediment availability and the contribution of the grain size distribution in the bed to the weighting of the saturated sediment concentration. One or more bed composition layers are located underneath the bed surface layer and form the upper part of the erodible bed. The bottom layer is the base layer and contains an infinite amount of erodible sediment according to the initial grain size distribution. The base layer cannot be eroded, but can supply sediment to the other layers.

Fig. 2 Schematic of bed composition discretisation and advection scheme. Horizontal exchange of sediment may occur solely through the air that interacts with the bed surface layer. The detail presents the simulation of sorting and beach armoring where the bed surface layer in the upwind grid cell becomes coarser due to non-uniform erosion over the sediment fractions, while the bed surface layer in the downwind grid cell becomes finer due to non-uniform deposition over the sediment fractions. Symbols refer to Equations (1) and (2).

Each layer in each grid cell describes a grain size distribution over a predefined number of sediment fractions (Figure Fig. 2, detail). Sediment may enter or leave a grid cell only through the bed surface layer. Since the velocity threshold depends among others on the grain size, erosion from the bed surface layer will not be uniform over all sediment fractions, but will tend to erode fines more easily than coarse sediment (Figure Fig. 2, detail, upper left panel). If sediment is eroded from the bed surface layer, the layer is repleted by sediment from the lower bed composition layers. The repleted sediment has a different grain size distribution than the sediment eroded from the bed surface layer. If more fines are removed from the bed surface layer in a grid cell than repleted, the median grain size increases. If erosion of fines continues the bed surface layer becomes increasingly coarse. Deposition of fines or erosion of coarse material may resume the erosion of fines from the bed.

In case of deposition the process is similar. Sediment is deposited in the bed surface layer that then passes its excess sediment to the lower bed layers (Figure Fig. 2, detail, upper right panel). If more fines are deposited than passed to the lower bed layers the bed surface layer becomes increasingly fine.

## Simulation of the Emergence of Non-erodible Roughness Elements¶

Sediment sorting may lead to the emergence of non-erodible elements from the bed. Non-erodible roughness elements may shelter the erodible bed from wind erosion due to shear partitioning, resulting in a reduced sediment availability ([RGL93]). Therefore the equation of [RGL93] is implemented according to:

(6)$u_{\mathrm{* th, R}} = u_{\mathrm{* th}} \cdot \sqrt{ \left( 1 - m \cdot \sum_{k=k_0}^{n_{\mathrm{k}}}{w_k^{\mathrm{bed}}} \right) \left( 1 + \frac{m \beta}{\sigma} \cdot \sum_{k=k_0}^{n_{\mathrm{k}}}{w_k^{\mathrm{bed}}} \right) }$

in which $$\sigma$$ is the ratio between the frontal area and the basal area of the roughness elements and $$\beta$$ is the ratio between the drag coefficients of the roughness elements and the bed without roughness elements. $$m$$ is a factor to account for the difference between the mean and maximum shear stress and is usually chosen 1.0 in wind tunnel experiments and may be lowered to 0.5 for field applications. The roughness density $$\lambda$$ in the original equation of [RGL93] is obtained from the mass fraction in the bed surface layer $$w_k^{\mathrm{bed}}$$ according to:

$\lambda = \frac{\sum_{k=k_0}^{n_{\mathrm{k}}}{w_k^{\mathrm{bed}}}}{\sigma}$

in which $$k_0$$ is the index of the smallest non-erodible sediment fraction in current conditions and $$n_{\mathrm{k}}$$ is the total number of sediment fractions. It is assumed that the sediment fractions are ordered by increasing size. Whether a fraction is erodible depends on the sediment transport capacity.

## Simulation of the Hydraulic Mixing, Infiltration and Evaporation¶

As sediment sorting due to aeolian processes can lead to armoring of a beach surface, mixing of the beach surface or erosion of course material may undo the effects of armoring. To ensure a proper balance between processes that limit and enhance sediment availability in the model both types of processes need to be sufficiently represented when simulating spatiotemporal varying bed surface properties and sediment availability.

A typical upwind boundary in coastal environments during onshore winds is the water line. For aeolian sediment transport the water line is a zero-transport boundary. In the presence of tides, the intertidal beach is flooded periodically. Hydraulic processes like wave breaking mix the bed surface layer of the intertidal beach, break the beach armoring and thereby influence the availability of sediment. Moreover, the hydraulic processes periodically wet the intertidal beach temporally increasing the shear velocity threshold. Infiltration and evaporation subsequently dry the beach.

In the model the mixing of sediment is simulated by averaging the sediment distribution over the depth of disturbance ($$\Delta z_{\mathrm{d}}$$). The depth of disturbance is linearly related to the breaker height (e.g. [Kin51], [Wil71], [MAROHare07]). [MAROHare07] proposes an empirical factor $$f_{\Delta z_{\mathrm{d}}}$$ [-] that relates the depth of disturbance directly to the local breaker height according to:

$\Delta z_{\mathrm{d}} = f_{\Delta z_{\mathrm{d}}} \cdot \min \left ( H \quad ; \quad \gamma \cdot d \right )$

in which the offshore wave height $$H$$ [m] is taken as the local wave height maximized by a maximum wave height over depth ratio $$\gamma$$ [-]. $$d$$ [m] is the water depth that is provided to the model through an input time series of water levels. Typical values for $$f_{\Delta z_{\mathrm{d}}}$$ are 0.05 to 0.4 and 0.5 for $$\gamma$$.

The drying of the beach is simulated by simplified functions for infiltration and evaporation. Infiltration is represented by an exponential decay function that is governed by a drying time scale $$T_{\mathrm{dry}}$$. Evaporation is simulated using an adapted version of the Penman-Monteith equation ([Shu93]) that is governed by meteorological time series of solar radiation, temperature and humidity.

Bibliography

 [Bag37] RA Bagnold. The size-grading of sand by wind. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, pages 250–264, 1937.
 [Bag37b] (1, 2) RA Bagnold. The transport of sand by wind. Geographical journal, pages 409–438, 1937.
 [Bel64] P Y Belly. Sand movement by wind. Technical Report 1, U.S. Army Corps of Engineers CERC, 1964. 38 pp.
 [dVvTdVvR+14] (1, 2, 3, 4) S de Vries, J S M van Thiel de Vries, L C van Rijn, S M Arens, and R Ranasinghe. Aeolian sediment transport in supply limited situations. Aeolian Research, 12:75–85, 2014. doi:10.1016/j.aeolia.2013.11.005.
 [DF10] I Delgado-Fernandez. A review of the application of the fetch effect to modelling sand supply to coastal foredunes. Aeolian Research, 2:61–70, 2010. doi:10.1016/j.aeolia.2010.04.001.
 [Kin51] C. A. M. King. Depth of disturbance of sand on sea beaches by waves. Journal of Sedimentary Petrology, 21(3):131–140, 1951. URL: http://archives.datapages.com/data/sepm/journals/v01-32/data/021/021003/pdfs/0131.pdf.
 [MAROHare07] (1, 2) G Masselink, N Auger, P Russell, and T O’Hare. Short-term morphological change and sediment dynamics in the intertidal zone of a macrotidal beach. Sedimentology, 54:39–53, 2007. doi:10.1111/j.1365-3091.2006.00825.x.
 [NE81] W G Nickling and M Ecclestone. The effects of soluble salts on the threshold shear velocity of fine sand. Sedimentology, 28:505–510, 1981.
 [PHN13] S. D. Peckham, E. W. H. Hutton, and B. Norris. A component-based approach to integrated modeling in the geosciences: the design of CSDMS. Computers and Geosciences, 53:3–12, 2013. doi:10.1016/j.cageo.2012.04.002.
 [PT90] K. Pye and H. Tsoar. Aeolian Sand and Sand Dunes. Unwin Hyman, London, 1990.
 [RGL93] (1, 2, 3) MR Raupach, DA Gillette, and JF Leys. The effect of roughness elements on wind erosion threshold. Journal of Geophysical Research: Atmospheres, 98(D2):3023–3029, 1993. doi:10.1029/92JD01922.
 [Shu93] W J Shuttleworth. Evaporation. In D R Maidment, editor, Handbook of Hydrology, pages 4.1–4.53. McGraw-Hill, New York, 1993.
 [Wil71] A. T. Williams. An analysis of some factors involved in the depth of disturbance of beach sand by waves. Marine Geology, 11(3):145–158, 1971. doi:10.1016/0025-3227(71)90003-X.
 [Delft3DFManual14] Delft3D-FLOW Manual. Delft3D - 3D/2D modelling suite for integral water solutions - Hydro-Morphodynamics. Deltares, Delft, May 2014. Version 3.15.34158.
 [SimunekSejnavG98] J. Šimůnek, M. Šejna, and M. Th. van Genuchten. The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably- saturated media. International Ground Water Modeling Center, Colorado School of Mines, Golden, Colorado, version 1.0. igwmc - tps - 70 edition, 1998. 186pp.