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..

Advection Scheme

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):

(4)\[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:

(5)\[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:

(6)\[\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:

(7)\[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:

(8)\[\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 (7) 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 (7). 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.

_images/bed_interaction_parameter.pdf

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 (5) 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.

_images/bedcomposition.pdf

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:

(9)\[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:

(10)\[\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

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.

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:

(11)\[\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\).

Simulation of surface moisture

Wave runup, capillary rise from the beach groundwater, and precipitation periodically wet the intertidal beach temporally increasing the shear velocity threshold ( Fig. 3). Infiltration and evaporation subsequently dry the beach.

_images/moisture_processes.jpg

Fig. 3 Illustration of processes influencing the volumetric moisture content \(\theta\) at the beach surface.

The structure of the surface moisture module and included processes are schematized in Fig. 4. The resulting surface moisture is obtained by selecting the largest of the moisture contents computed with the water balance approach (right column) and due to capillary rise from the groundwater table (left column). The method is based on the assumption that the flow of soil water is small compared to the flow of groundwater and that the beach groundwater dynamics primarily is controlled by the water level and wave action at the seaward boundary ([RGE99], [Sch14]). Thus, there is no feedback between the processes in the right column of Fig. 4 and the groundwater dynamics described in the left column.

_images/moisture_scheme.jpg

Fig. 4 Implementation of surface moisture processes in the AeoLiS.

Runup and wave setup

The runup height and wave setup are computed using the Stockdon formulas ([SHHS06]). Their parameterization differs depending on the dynamic beach steepness expressed through the Irribaren number:

(12)\[\xi = \tan \beta /\sqrt {{H_0}/{L_0}}\]

where \({H_0}\) is the significant offshore wave height, \({L_0}\) is the deepwater wavelength, and \({\tan \beta}\) is the foreshore slope.

For dissipative conditions, \({\xi}\) < 0.3, the runup, \({R_2}\), is parameterized as,

(13)\[{R_2} = 0.043\sqrt {{H_0}{L_0}}\]

and wave setup:

(14)\[< \eta > = 0.02\sqrt {{H_0}{L_0}}\]

For \({\xi}\) > 0.3, runup is paramterized as,

(15)\[{R_2} = 1.1\left( {0.35\beta \sqrt {{H_0}{L_0}} + \frac{{\sqrt {{H_0}{L_0}\left( {0.563{\beta ^2} + 0.004} \right)} }}{2}} \right)\]

and wave setup:

(16)\[< \eta > = 0.35\xi\]

Tide- and wave-induced groundwater variations

Groundwater under sandy beaches can be considered as shallow aquifers, with only horizontal groundwater flow so that the pressure distribution is hydrostatic ([BMH98], [BSDR19], [Nie90], [RGE99]). The cross-shore flow dominates temporal variations of groundwater levels. Alongshore, groundwater table variations are typically small ([Sch14]). Although the surface moisture model can be extended over a two-dimensional grid, the groundwater simulations are performed for 1D transects cross-shore to avoid numerical instabilities at the seaward boundary and reduce computational time.

The beach aquifers is schematised as a sandy body, with saturated hydraulic conductivity, \(K\), and effective porosity, \({n_e}\). The aquifer is assumed to rest on an impermeable surface, where \(D\) is the aquifer depth. The groundwater elevation relative to the mean sea level (MSL) is denoted \(\eta\), and the shore-perpendicular x-axis is positive landwards, with an arbitrary starting point. The sand is assumed to be homogenous and isotropic. In this context, isotropy implies that hydraulic conductivity is independent of flow direction.

The horizontal groundwater discharge per unit area, \(u\), is then governed by Darcy’s law,

(17)\[u = - K\frac{{\partial \eta }}{{\partial x}}\]

and the continuity equation (see e.g., [Nie09]),

(18)\[\frac{{\partial \eta }}{{\partial t}} = - \frac{1}{{{n_e}}}\frac{\partial }{{\partial x}}((D + \eta )u)\]

where \(t\) is time.

The groundwater overheight due to runup, \({U_l}\), is computed by ([KNH94], [NDWE88]),

(19)\[\begin{split}{U_l} = \left\{ \begin{gathered}{C_l}Kf(x)\,\,\,\,{\text{if }}{x_S} \leqslant x \leqslant {x_R} \hfill \\0,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{\text{if }}x > {x_R} \hfill \\\end{gathered} \right.\end{split}\]

where \({C_l}\) is an infiltration coefficient (-), and \(f(x)\) is a function of \(x\) ranging from 0 to 1. \({x_S}\) is the horizontal location of the sum of the still water level and wave setup, and \({x_R}\) is the horizontal location of the runup limit:

(20)\[\begin{split}f(x) = \left\{ \begin{gathered} \frac{{x - {x_s}}}{{\frac{2}{3}\left( {{x_{ru}} - {x_s}} \right)}}\,\,\,\,\,\,\,\,\,\,\,\,\,if\,{x_s} < x \leqslant {x_s} + \frac{2}{3}\left( {{x_{ru}} - {x_s}} \right)\, \hfill \\ 3 - \frac{{x - {x_s}}}{{\frac{1}{3}\left( {{x_{ru}} - {x_s}} \right)}}\,\,\,\,\,if\,{x_s} + \frac{2}{3}\left( {{x_{ru}} - {x_s}} \right)\, < x < {x_{ru}} \hfill \\ \end{gathered} \right.\end{split}\]

Substitution of \(u\) (Equation (17)) in the continuity equation (Equation (18)) with the addition of \({U_l}/{n_e}\) gives the nonlinear Boussinesq equation:

(21)\[\frac{{\partial \eta }}{{\partial t}} = \frac{K}{{{n_e}}}\frac{\partial }{{\partial x}}\left( {(D + \eta )\frac{{\partial \eta }}{{\partial x}}} \right) + \frac{{{U_l}}}{{{n_e}}}\]

Capillary rise

Soil water retention (SWR) functions describe the surface moisture due to capillary transport of water from the groundwater table ([vG80]):

(22)\[\theta (h) = {\theta _r} + \frac{{{\theta _s} - {\theta _r}}}{{{{\left[ {1 + {{\left| {\alpha h} \right|}^n}} \right]}^m}}}\]

where \(h\) is the groundwater table depth, \(\alpha\) and \(n\) are fitting parameters related to the air entry suction and the pore size distribution. The parameter \(m\) is commonly parameterised as \(m = 1 - 1/n\).

The resulting surface moisture is computed for both drying and wetting conditions, i.e., including the effect of hysteresis. The moisture contents computed with drying and wetting SWR functions are denoted \({\theta ^d}(h)\) and \({\theta ^w}(h)\), respectively. When moving between wetting and drying conditions, the soil moisture content follows an intermediate retention curve called a scanning curve. The drying scanning curves are scaled from the main drying curve and wetting scanning curves from the main wetting curve. The drying scanning curve is then obtained from ([Mua74]):

(23)\[{\theta ^d}({h_\Delta },h) = {\theta ^w}(h) + \frac{{\left[ {{\theta ^w}({h_\Delta }) - {\theta ^w}(h)} \right]}}{{\left[ {{\theta _s} - {\theta ^w}(h)} \right]}}\left[ {{\theta ^d}(h) - {\theta ^w}(h)} \right]\]

where \({h_\Delta}\) is the groundwater table depth at the reversal on the wetting curve.

The wetting scanning curve is obtained from ([Mua74]):

(24)\[{\theta ^w}({h_\Delta },h) = {\theta ^w}(h) + \frac{{\left[ {{\theta _s} - {\theta ^w}(h)} \right]}}{{\left[ {{\theta _s} - {\theta ^w}({h_\Delta })} \right]}}\left[ {{\theta ^d}({h_\Delta }) - {\theta ^w}({h_\Delta })} \right]\]

where \({h_\Delta}\) is the groundwater table depth at the reversal on the drying curve.

Infiltration

Infiltration is accounted for by assuming that excess water infiltrates until the moisture content reaches field capacity, \({\theta_fc}\). The moisture content at field capacity is the maximum amount of water that the unsaturated zone of soil can hold against the pull of gravity. For sandy soils, the matric potential at this soil moisture condition is around - 1/10 bar. In equilibrium, this potential would be exerted on the soil capillaries at the soil surface when the water table is about 100 cm below the soil surface, \({\theta _{fc}} = {\theta ^d}(100)\).

Infiltration is represented by an exponential decay function that is governed by a drying time scale \(T_{\mathrm{dry}}\). Exploratory model runs of the unsaturated soil with the HYDRUS1D ([vSimrunekvSejnavG98]) hydrology model show that the increase of the volumetric water content to saturation is almost instantaneous with rising tide. The drying of the beach surface through infiltration shows an exponential decay. In order to capture this behavior the volumetric water content is implemented according to:

(25)\[\frac{{d\theta }}{{dt}} = \left( {\theta - {\theta _{fc}}} \right)\left( {{e^{ - \ln (2)\frac{{dt}}{{{T_{dry}}}}}}} \right)\]

An alternative formulation is used for simulations that does not account for ground water and SWR processes,

(26)\[\begin{split}p_{\mathrm{V}}^{n+1} = \left\{ \begin{array}{ll} p & \mathrm{if} ~ \eta > z_{\mathrm{b}} \\ p_{\mathrm{V}}^n \cdot e^{\frac{\log \left( 0.5 \right)}{T_{\mathrm{dry}}} \cdot \Delta t^n} - E_{\mathrm{v}} \cdot \frac{\Delta t^n}{\Delta z} & \mathrm{if} ~ \eta \leq z_{\mathrm{b}} \\ \end{array} \right.\end{split}\]

where \(\eta\) [m+MSL] is the instantaneous water level, \(z_{\mathrm{b}}\) [m+MSL] is the local bed elevation, \(p_{\mathrm{V}}^n\) [-] is the volumetric water content in time step \(n\), \(\Delta t^n\) [s] is the model time step and \(\Delta z\) is the bed composition layer thickness. \(T_{\mathrm{dry}}\) [s] is the beach drying time scale, defined as the time in which the beach moisture content halves.

Precipitation and evaporation

A water balance approach accounts for the effect of precipitation and evaporation,

(27)\[\frac{{d\theta }}{{dt}} = \frac{{\left( {P - E} \right)\,}}{{\Delta z}}\,\]

where \(P\) is the precipitation, \(E\) is the evaporation, and \(\Delta z\) is the thickness of the surface layer.

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.

\(E_{\mathrm{v}}\) [m/s] is the evaporation rate that is implemented through an adapted version of the Penman equation ([Shu93]):

(28)\[E_{\mathrm{v}} = \frac{m_{\mathrm{v}} \cdot R_{\mathrm{n}} + 6.43 \cdot \gamma_{\mathrm{v}} \cdot (1 + 0.536 \cdot u_2) \cdot \delta e} {\lambda_{\mathrm{v}} \cdot (m_{\mathrm{v}} + \gamma_{\mathrm{v}})} \cdot 9 \cdot 10^7\]

where \(m_{\mathrm{v}}\) [kPa/K] is the slope of the saturation vapor pressure curve, \(R_{\mathrm{n}}\) [\(\mathrm{MJ/m^2/day}\)] is the net radiance, \(\gamma_{\mathrm{v}}\) [kPa/K] is the psychrometric constant, \(u_2\) [m/s] is the wind speed at 2 m above the bed, \(\delta e\) [kPa] is the vapor pressure deficit (related to the relative humidity) and \(\lambda_{\mathrm{v}}\) [MJ/kg] is the latent heat vaporization. To obtain an evaporation rate in [m/s], the original formulation is multiplied by \(9 \cdot 10^7\).

Shear velocity threshold

The shear velocity threshold represents the influence of bed surface properties in the saturated sediment transport equation. The shear velocity threshold is computed for each grid cell and sediment fraction separately based on local bed surface properties, like moisture, roughness elements and salt content. For each bed surface property supported by the model a factor is computed to increase the initial shear velocity threshold:

(29)\[u_{\mathrm{* th}} = f_{u_{\mathrm{* th}}, \mathrm{M}} \cdot f_{u_{\mathrm{* th}}, \mathrm{R}} \cdot f_{u_{\mathrm{* th}}, \mathrm{S}} \cdot u_{\mathrm{* th, 0}}\]

The initial shear velocity threshold \(u_{\mathrm{* th, 0}}\) [m/s] is computed based on the grain size following [Bag37a]:

(30)\[u_{\mathrm{* th, 0}} = A \sqrt{ \frac{\rho_{\mathrm{p}} - \rho_{\mathrm{a}}}{\rho_{\mathrm{a}}} \cdot g \cdot d_{\mathrm{n}}}\]

where \(A\) [-] is an empirical constant, \(\rho_{\mathrm{p}}\) [\(\mathrm{kg/m^3}\)] is the grain density, \(\rho_{\mathrm{a}}\) [\(\mathrm{kg/m^3}\)] is the air density, \(g\) [\(\mathrm{m/s^2}\)] is the gravitational constant and \(d_{\mathrm{n}}\) [m] is the nominal grain size of the sediment fraction.

Moisture content

The shear velocity threshold is updated based on moisture content following [Bel64]:

(31)\[f_{u_{\mathrm{* th}}, \mathrm{M}} = \max(1 \quad ; \quad 1.8 + 0.6 \cdot \log(p_{\mathrm{g}}))\]

where \(f_{u_{\mathrm{* th},M}}\) [-] is a factor in Equation (29), \(p_{\mathrm{g}}\) [-] is the geotechnical mass content of water, which is the percentage of water compared to the dry mass. The geotechnical mass content relates to the volumetric water content \(p_{\mathrm{V}}\) [-] according to:

\[ \begin{align}\begin{aligned}:label: vol-water\\p_{\mathrm{g}} = \frac{p_{\mathrm{V}} \cdot \rho_{\mathrm{w}}}{\rho_{\mathrm{p}} \cdot (1 - p)}\end{aligned}\end{align} \]

where \(\rho_{\mathrm{w}}\) [\(\mathrm{kg/m^3}\)] and \(\rho_{\mathrm{p}}\) [\(\mathrm{kg/m^3}\)] are the water and particle density respectively and \(p\) [-] is the porosity. Values for \(p_{\mathrm{g}}\) smaller than 0.005 do not affect the shear velocity threshold ([PT90]). Values larger than 0.064 (or 10% volumetric content) cease transport ([DF10]), which is implemented as an infinite shear velocity threshold.

Roughness elements

The shear velocity threshold is updated based on the presence of roughness elements following [RGL93]:

\[ \begin{align}\begin{aligned}:label: shear-rough\\f_{u_{\mathrm{* th},R}} = \sqrt{(1 - m \cdot \sum_{k=k_0}^{n_k}{\hat{w}_k^{\mathrm{bed}}}) (1 + \frac{m \beta}{\sigma} \cdot \sum_{k=k_0}^{n_k}{\hat{w}_k^{\mathrm{bed}}})}\end{aligned}\end{align} \]

by assuming:

\[ \begin{align}\begin{aligned}:label: lambda-rough\\\lambda = \frac{\sum_{k=k_0}^{n_k}{\hat{w}_k^{\mathrm{bed}}}}{\sigma}\end{aligned}\end{align} \]

where \(f_{u_{\mathrm{* th},R}}\) [-] is a factor in Equation (29), \(k_0\) is the sediment fraction index of the smallest non-erodible fraction in current conditions and \(n_k\) is the number of sediment fractions defined. The implementation is discussed in detail in section ref{sec:roughness}.

Salt content

The shear velocity threshold is updated based on salt content following [NE81]:

(32)\[f_{u_{\mathrm{* th}},S} = 1.03 \cdot \exp(0.1027 \cdot p_{\mathrm{s}})\]

where \(f_{u_{\mathrm{* th},S}}\) [-] is a factor in Equation (29) and \(p_{\mathrm{s}}\) [-] is the salt content [mg/g]. Currently, no model is implemented that predicts the instantaneous salt content. The spatial varying salt content needs to be specified by the user, for example through the BMI interface.

Bibliography

vSimrunekvSejnavG98

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.

Bag37a

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.

BMH98

Andrew J. Baird, Travis Mason, and Diane P. Horn. Validation of a Boussinesq model of beach ground water behaviour. Marine Geology, 148(1-2):55–69, 1998. doi:10.1016/S0025-3227(98)00026-7.

Bel64

P Y Belly. Sand movement by wind. Technical Report 1, U.S. Army Corps of Engineers CERC, 1964. 38 pp.

BSDR19

Laura B. Brakenhoff, Yvonne Smit, Jasper J.A. Donker, and Gerben Ruessink. Tide-induced variability in beach surface moisture: Observations and modelling. Earth Surface Processes and Landforms, 44(1):317–330, 2019. doi:10.1002/esp.4493.

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.

KNH94

H.Y. Kang, P. Nielsen, and D.J Hanslow. Watertable overheight due to wave runup on a sandy beach. In Coastal Engineering 1994, 2115–2124. 1994.

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.

Mua74(1,2)

Y Mualem. A Conceptual Model of Hysteresis. Water Resources Research, 10(3):514–520, 1974.

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.

NDWE88

P Nielsen, GA Davis, JM Winterbourne, and G Elias. Wave setup and the watertable in sandy beaches, Technical Report Technical Memo- randum 88/1, New South Wales Public Works Department Coastal Branch. Technical Report, -, 1988.

Nie90

Peter Nielsen. Tidal dynamics of the water table in beaches. Water Resources Research, 26(9):2127–2134, 1990. doi:10.1029/WR026i009p02127.

Nie09

Peter Nielsen. Coastal and estuarine processes. World Scientific Publishing Company, 2009.

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.

RGE99(1,2)

B. Raubenheimer, R. T. Guza, and Steve Elgar. Tidal water table fluctuations in a sandy ocean beach. Water Resources Research, 35(8):2313–2320, 1999. doi:10.1029/1999WR900105.

RGL93(1,2,3,4)

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.

Sch14(1,2)

Phillip P Schmutz. Investigation of Factors Controlling the Dynamics of Beach Surface Moisture. PhD thesis, Louisiana State University, 2014.

Shu93(1,2)

W J Shuttleworth. Evaporation. In D R Maidment, editor, Handbook of Hydrology, pages 4.1–4.53. McGraw-Hill, New York, 1993.

SHHS06

Hilary F. Stockdon, Rob A. Holman, Peter A. Howd, and Asbury H. Sallenger. Empirical parameterization of setup, swash, and runup. Coastal Engineering, 53(7):573–588, 2006. URL: https://www.sciencedirect.com/science/article/pii/S0378383906000044, doi:https://doi.org/10.1016/j.coastaleng.2005.12.005.

vG80

M. Th van Genuchten. Closed-Form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. Soil Science Society of America Journal, 44(5):892–898, 1980. doi:10.2136/sssaj1980.03615995004400050002x.

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.