Modeling the Cross‐Sectional Dynamics of Tidal Sandbanks in Sediment‐Scarce Conditions

Tidal sandbanks are large‐scale dynamic bedforms that consist of sandy sediment. They have been observed in shallow seas with varying sediment supply, including sediment‐scarce environments like the Flemish Banks, Zeeland Banks, and Norfolk Banks. However, we do not yet understand how scarcity affects sandbank evolution. Therefore, we have developed an idealized nonlinear process‐based model with the aim of studying cross‐sectional shape and migration under sediment‐scarce conditions. Scarcity is included through a non‐erodible layer from which no sediment can be entrained. Our results show that bank height and width decrease when the sediment budget decreases. The bank height is more sensitive to scarcity than bank width. Furthermore, sand scarcity decreases (and may even reverse) bank asymmetry and increases migration rate when a residual current is present. The migration rate attains a maximum for a specific sediment budget, which is controlled by the ratio of the cross‐sectional area of the sandbank and the tide‐averaged sediment flux. Our findings show that sandbank dynamics are strongly affected by scarcity, which is critical for seas with receding sediment stocks (e.g., through extraction).


Introduction
Tidal sandbanks are tide-driven marine bedforms on the bottom of shallow shelf seas with sandy beds.They are large elongated bodies of sand with crests and troughs.Their dimensions can be up to tens of meters in height, up to kilometers in width, and up to tens of kilometers in length (de Swart & Yuan, 2019;Dyer & Huntley, 1999).Their crest line is typically oriented anticlockwise with respect to the tidal flow on the Northern Hemisphere (Kenyon et al., 1981).Furthermore, sandbanks typically occur in patches of multiple parallel banks (Knaapen, 2009) and exhibit dynamic behavior on time scales of centuries and larger.These dynamics include cross-sectional growth and migration.
The presence of sandbanks shapes a unique bathymetry with ecological and economic benefits.They attract marine flora and fauna and provide a food source for birds (Atalah et al., 2013;Wyns et al., 2021).Therefore, several sandbanks in the North Sea have been designated as marine protected areas.Sandbanks can also serve as a foundation for wind farms (Velenturf et al., 2021) and those nearshore may contribute to coastal protection (Dolphin et al., 2007).Finally, sandbanks are frequently targeted for marine aggregate extraction to satisfy the rising global demand for sand (Degrendele et al., 2010;Torres et al., 2021;Van Lancker et al., 2010).
The sediment supply varies naturally across seas.Sandbanks have been observed at locations that are rich and scarce in sediments.We define an area as sediment-scarce when the sand availability is a constraint for sandbank evolution.Field observations have identified several sandbanks where the bed between banks is not covered by sand.For example, the Flemish Banks and Zeeland Banks expose a layer of hard clay in the troughs (Hademenos et al., 2019;Kint et al., 2021), and the bed between the Norfolk Banks consists of gravel (Caston, 1972).In order to understand sandbank dynamics, let us review the literature on sandbank modeling.
The formation of sandbanks has been explained as a free instability of a sandy bed subject to tidal flow (de Vriend, 1990;Hulscher et al., 1993;Huthnance, 1982b).Linear stability analysis of idealized morphodynamic models have shown how sandbank patches may grow from small topographic perturbations of the bed.They likely form with a preferred wavelength and orientation with respect to the principal tidal flow, which is known as the "fastest growing mode" (Hulscher, 1996;Hulscher et al., 1993).Bank growth is driven by horizontal circulations that result from interactions between topography, bed friction, tidal motion, and Coriolis effects, as has been observed in the field (Caston & Stride, 1970), studied theoretically (Huthnance, 1973;Pattiaratchi & Collins, 1987;Robinson, 1981;Zimmerman, 1981), and reproduced in complex numerical models (e.g., Sanay et al., 2007).The initial instability may be triggered by natural topographic variation or excavation of the seabed (Nnafie et al., 2020;Roos et al., 2008).Importantly, linear stability analysis is only applicable when bank amplitudes are low with respect to the mean water depth.
Sandbank evolution towards an equilibrium has been explained by idealised models that include nonlinear tidetopography interactions.Sandbanks may reach a static or dynamic equilibrium.Static equilibria refer to sandbanks that do not move, whereas dynamic equilibria refer to sandbanks that migrate (Roos et al., 2004) or exhibit oscillating bank ends (Yuan et al., 2017).Huthnance (1982aHuthnance ( , 1982b) ) found equilibrium cross-sections under sediment scarcity and simplified hydrodynamics.He omitted inertia and Coriolis terms, and reduced the tidal cycle to two steady flows in opposite directions.Under these simplified conditions, sediment scarcity was a prerequisite for equilibria to exist in his model.Roos et al. (2004) obtained equilibrium cross-sectional profiles for static and migrating banks without sediment scarcity under realistic tidal conditions.They assumed no along-crest variation and included depth-dependent wave stirring.Tambroni and Blondeaux (2008) also obtained equilibrium bank cross-sectional shapes under infinite sediment supply and an elliptic tidal flow, which enabled a weakly nonlinear analysis.Dynamic equilibria for bank shapes that varied in both horizontal dimensions were obtained by Yuan et al. (2017) with bank ends periodically oscillating in time.Besides equilibrium conditions, long-term  2023) for tidal dunes and on Kleinhans et al. (2002), Tuijnder et al. (2009), Dreano et al. (2010), Endo (2016), andVah et al. (2020) for fluvial dunes.
nonlinear idealised models have also enhanced our understanding of bank formation (Yuan et al., 2016), their response to sea level rise (Yuan &de Swart, 2017), andbank-breaking (van Veelen et al., 2018).Besides the equilibrium shapes by Huthnance (1982aHuthnance ( , 1982b)), the effect of sand scarcity on sandbank evolution has not been considered in any of these studies.Detailed site-specific modeling studies in combination with measurements can provide insight in the local hydrodynamics and sediment transport patterns.Examples include the studies into the Middelkerke Banks (Pan et al., 2007;Williams et al., 2000), Hinder Banks (Deleu et al., 2004), Great Yarmouth Banks (Horrillo-Caraballo & Reeve, 2008), Kwinte Bank (Brière et al., 2010;Van den Eynde et al., 2010), and the Scarweather and Nash Sands banner banks (Fairley et al., 2016;Mitchell et al., 2021).While these studies provide valuable insights in the site-specific processes, including the Hinder Banks and Kwinte Bank which experience sediment scarcity, they typically focus on dynamics on the timescale of a tidal cycle rather than the centennial timescale of sandbank evolution.Nnafie et al. (2020) did model the long-term impacts of sand extraction and land reclamation on the Belgian Continental Shelf where the Hinder and Kwinte Banks are located, but did not consider the sediment scarcity as observed in the field.
Although the reviewed studies provide a comprehensive understanding of sandbank dynamics, the effects of sediment scarcity remain unclear.Huthnance (1982aHuthnance ( , 1982b) ) implemented scarcity as a non-erodible layer and hypothesized that scarcity was required for equilibrium cross-sections to exist, but later studies with more realistic hydrodynamics (e.g., Roos et al., 2004;Yuan et al., 2017) showed that this was not the case.Regarding other bedforms, sand scarcity was shown to decrease height and width, increase wavelength and three-dimensionality in shape, and affect the migration rate of tidal (Nnafie et al., 2023;Porcile et al., 2017) and fluvial dunes (Dreano et al., 2010;Endo, 2016;Kleinhans et al., 2002;Tuijnder et al., 2009;Vah et al., 2020).However, these bedforms are driven by flow circulations in the vertical plane (Hulscher, 1996;Paarlberg et al., 2009;Vittori & Blondeaux, 2020) rather than the residual circulations in the horizontal plan that drive tidal sandbanks (Figure 1).Therefore, we do not know how scarcity affects the long-term evolution, equilibrium shape, and migration rate of tidal sandbanks under realistic tidal conditions, despite their presence in sediment-scarce areas such as the Flemish Banks, Zeeland Banks, and Norfolk Banks.
Here, we aim to understand how sediment scarcity affects the cross-sectional shape and migration rate of tidal sandbanks in morphodynamic equilibrium.We extend the idealised nonlinear process-based model by Roos et al. (2004) with a non-erodible layer from which no sediment can be entrained, inspired by Huthnance (1982aHuthnance ( , 1982b)), but drop the wind wave stirring mechanism as we will show that this is not required for equilibria to exist.Using the extended model, we simulate the evolution of a small undulation of the bed towards equilibrium under symmetric and asymmetric tidal conditions.Analogous to Roos et al. (2004), we consider sandbanks without along-crest variation.This manuscript is organized as follows.Section 2 describes our model, including governing equations and a scaling procedure.Next, Section 3 explains the solution procedure of our model.The results are presented in Section 4, and a discussion hereof is the topic of Section 5. Finally, our conclusions are summarized in Section 6.

Model Geometry
We consider a sediment-scarce part of a shallow sea, which is unaffected by coastal boundaries (Figure 2).We define a vertical coordinate z* (the asterisk denotes an unscaled parameter), which points upwards, and two horizontal coordinates x* and y*.The x*-axis is directed perpendicular to the bank crest (cross-bank) and the y*axis is directed parallel to the bank crest (along-bank).Our study area is characterized by mean a water depth H*, a sand layer with thickness D * sand , and the presence of a non-erodible layer (e.g., rock or hard clay) at depth z * = H * ne .We define a buffer layer of thickness δ* on top of the non-erodible layer in which the effect of the nonerodible layer is already felt, which is further motivated in Section 2.2.The free water surface is located at z* = ζ* (x*, t*).Its space-averaged elevation equals z* = 0.The study area is subject to tidal forcing, which drives a tidal flow.The velocity components are u* in x*-direction and v* in y*-direction.
As sandbanks are much longer than they are wide, the dynamics over the cross-section are dominant for sandbank evolution.Therefore, we employ a topography z* = h*(x*, t*) that varies only in cross-bank direction.Note that the ambient tidal flow may approach under an angle θ with respect to the crest line in the horizontal plane.The domain length L * dom and orientation θ match the wavelength and orientation of the fastest growing mode from linear stability analysis (e.g., Hulscher et al., 1993), which is described in Section 3.1.Due to the choice of domain length, a single sandbank typically develops within the domain.We employ spatially periodic boundary conditions, such that we effectively simulate an infinite patch of identical sandbanks.
The cross-sectional profile of a sandbank is characterized by four shape parameters (Figure 2b).z * crest is the elevation of the highest point of the bank topography.z * trough is the elevation of the lowest point of the bank topography.W * bank is the width of the sandbank measured at the undisturbed bed level (z* = H*).A bank is the bank asymmetry is given by A bank = ln l * 1 / l * 2 ) , where l * 1 is the length of the stoss side and l * 2 is the length of the lee side.l * 1 and l * 2 are measured between the crest and the bank toes with positions x * 1 and x * 2 .The bank toes are located at the points where h * = H * ne δ * , or if these do not exist, at the lowest point in the topography.
The difference between z * crest and z * trough is referred to as the bank height.Furthermore, the migration speed, positive in x*-direction, is denoted by c * mig .Its evaluation will be addressed in Section 3.6.

Governing Equations
The hydrodynamics are governed by the depth-averaged shallow water equations, according to Equations 1 and 2 are the momentum equations in x* and y*-direction respectively, and Equation 3 is the continuity equation.f* = 2Ω* sin φ denotes the Coriolis parameter (with the angular frequency of the Earth's rotation Ω* = 7.292 × 10 5 rad s 1 and latitude ϑ).r * = C D linearization (e.g., Huthnance, 1982aHuthnance, , 1982b;;Roos et al., 2004;Zimmerman, 1982) with drag coefficient C D and tidal velocity scale U*.Furthermore, g* = 9.81 m s 2 is the gravitational acceleration and P * x ,P * y ) are spatially uniform yet time-dependent forcing terms, to be specified in Section 2.4.
We model suspended sediment transport via an advection equation (Schuttelaars & de Swart, 1996), expressed as where c*(x*, t*) is the depth-averaged suspended sediment concentration in m (strictly m 3 m 2 ), c * e (x * ,t * ) is the entrainment concentration in m, and γ* is a deposition coefficient in s 1 .The left-hand side of Equation 4 describes the spatiotemporal relaxation of the sediment concentration, similar to the spatial relaxation under steady flow as in Kroy et al. (2002) and Andreotti et al. (2012) but with an extra temporal relaxation term due to the oscillatory flow in our model.The right-hand side models the exchange between the bed and the water column.Sediment pick-up is controlled by the entrainment concentration c * e , according to Here, α * s is a sediment transport coefficient and, as a novel element in this study, μ s is a sediment transport limiter based on sediment scarcity, according to (6) The limiter is set at 1 when sediment is available and at 0 when the non-erodible layer is exposed.In the buffer layer, μ s varies smoothly between 0 and 1 with ∂μ s /∂h* = ∂ 2 μ s /∂h* 2 = 0 at the transition points (Figure 3), which improves the numerical stability of our solution procedure (Section 3).
Finally, the bed evolution is driven by the difference between sediment entrainment and deposition plus the gravity-driven down-slope transport, according to ϵ p = 0.4 is the bed porosity and λ* is a bed slope coefficient.Importantly, in all the above equations we have assumed uniformity in the y*-direction, that is, ∂/∂y* = 0.

Scaling Procedure
The model will be built using dimensionless quantities and equations.To this end, we introduce dimensionless coordinates, without an asterisk, according to Herein, σ* is the angular frequency of the M 2 -tide in rad/s, L* = U*/σ* is the tidal excursion length in m, and ) is the morphological time scale, approximately 207 years under the conditions studied here.The morphological time scale is much larger than the time scale of tidal motion (2π/σ* ≈ 12.42 hr).
Therefore, we can ignore bed evolution when resolving the tidal dynamics, known as the quasi-stationary approach, so that h is a function of x and τ only.Furthermore, the dimensionless variables are and the dimensionless parameters are Using the dimensionless coordinates, variables and parameters, our set of model equations becomes The angular brackets 〈•〉 in Equation 16denote averaging over a tidal cycle.Finally, we have adopted the rigid lid approach in which the contribution of the free surface elevation to the local water depth is assumed to be negligible.This assumption is supported by the small Froude Number Fr 2 = U* 2 /(g*H*) < 10 2 in this study.

Basic State: Tidal Flow Over a Flat Bed
As a reference state, needed to specify our forcing P x and P y , we consider the tidal flow over a flat bed, that is, unaffected by bank topography.Our desired reference tidal flow is a spatially uniform yet time-dependent M2tide, which may be supplemented by a unidirectional current (M0) and an M4-overtide.This tidal flow over a horizontal bed is prescribed as where the complex Fourier coefficients ( Ũ0p , Ṽ0p ), with accent ̃denoting a temporal expansion, are chosen such that the tidal flow satisfies (u 0 ,v 0 ) = U( sin θ, cos θ) j 0 + j 2 cos t + j 4 cos 2t φ 4 )).
(18) j 0 , j 2 , and j 4 are the dimensionless amplitudes of the M0, M2, and M4 tidal constituents respectively, and φ 4 is the phase difference between the M2 and M4 tidal constituents.The prescribed tidal flow is included in the momentum equations Equations 11 and 12 through forcing terms P x and P y , which read The coefficients ( Pxp , Pyp ) are obtained by substituting Equations 17 and 19 in Equations 11 and 12, with ∂/∂x = 0.This leads to 3. Solution Procedure

Initial Topography
The topography is expanded as a truncated spatial Fourier series, according to Herein, the accent ˇdenotes that ȟm are coefficients for spatial modes with corresponding subscript m.The total number of modes included is M and k m = mk min is the topographic wave number corresponding to each mode with k min = 2π/L dom being the minimum wave number.The initial bank topography is a sinusoidal shape with amplitude H init .Its corresponding topographic coefficients are ȟ0 = 1, ȟ1 (0) = ȟ 1 (0) = H init / 2, and ȟm (0) = 0 for |m| ≥ 2. H init is always smaller than D sand such that the non-erodible layer is never exposed at the start of the simulation.
Following Roos et al. (2004), we conduct a linear stability analysis to obtain the fastest growing mode (e.g., Hulscher et al., 1993).The length and orientation of the fastest growing mode serve as input for the domain length L dom and bank orientation θ, as pointed out in Section 2.1.Furthermore, we compute a linear growth rate ω m for each spatial mode m in Equation 22.These will be used later in the semi-implicit time stepping, which is described in Section 3.5.

Fourier Expansion
We define the vector of the hydrodynamic and sediment transport unknowns ϕ = (u,v,ζ,c e ,c), and expand ϕ into P temporal and M spatial Fourier modes, according to with complex spatiotemporal modes φpm = ûpm , ûpm , ζpm , ĉe,pm , ĉpm ), denoted by accent ̂and subscript pm.The expansion enables a spectral solution procedure of the hydrodynamics and a pseudospectral solution procedure of the sediment transport.

Hydrodynamic Solution
The hydrodynamic solution procedure follows Roos et al. (2004).We derive the cross-bank water flux ξ(t) by spatially integrating the continuity Equation 13 such that By substituting Equations 22-24 in the momentum Equations 11 and 12, followed by taking the spatial average of Equation 11 and multiplying Equation 12 by h, the hydrodynamic solution in Fourier space satisfies.
The brackets with subscript {•} 0 denote a spatial average, {hv} pm denotes a spatial convolution sum, and {ξv} pm denotes a temporal convolution sum.Equations 25 and 26 can be solved sequentially if f = 0 or iteratively if f ≠ 0.

Suspended Sediment Concentration
The suspended sediment concentration is obtained using a pseudo-spectral method.First, the entrainment limiter μ s is computed in physical space using Equation 6.Then, the entrainment concentration c e (Equation 14) is calculated in physical space, based on μ s and the hydrodynamic solution.The entrainment concentration c e is expanded following Equation 23.The concentration c is then obtained by solving the expanded advection equation in Fourier space, which reads in which {cu} pm represents a spatiotemporal convolution sum.

Bed Evolution
The bed evolution (Equation 16) over a morphological time step Δτ is computed using discrete semi-implicit time stepping (Roos et al., 2004), expressed as in which Bm = γ( ĉe,0m ĉ0m ) k m 2 λ{ĉ e,0 h} m is the growth of each spatial mode.Note that the tide-averaged complex Fourier coefficients are used for entrainment and deposition, which are denoted by subscript p = 0. ω m is the linear growth rate of spatial mode m (see Section 3.1).This semi-implicit time stepping initially favors the growth of lower order modes, which enhances the model stability.It does not affect the equilibrium shapes of the sandbanks which are the focus of this study.The morphological time step Δτ is adaptively limited to improve numerical stability under three conditions: (a) when the growth would erode the non-erodible layer, (b) when growth exceeds a critical value max ⃒ ⃒B m ⃒ ⃒ )> Bcrit , or (c) in the period following the first erosion of the buffer layer.
Specifically, the time step is multiplied by 0.1 when max(h δ ) = 0.The multiplier then increases linearly until it equals 1 at max(h δ ) = 0.5.

Equilibrium Conditions
The root-mean-square height of the bedforms h rms is a measure of the potential energy of the bedform (Garnier et al., 2006;Yuan et al., 2017).It is defined as with ȟ0 = 1 as the mean water depth.The relative change in potential energy Γ, scaled as Γ = Γ * T * mor , indicates how close the topography is to a morphodynamic equilibrium (Garnier et al., 2006;Yuan et al., 2017), and is computed via We consider the topography in equilibrium when |Γ| is smaller than a critical value Γ crit = 10 2 over a time span of Δτ = 5.The critical value has been selected for this study based on an equivalent dimensional value of 5 • 10 5 yr 1 and is applied to static and dynamic equilibria.Following Roos et al. (2004), we estimate the migration rate from the growth rate Bm (Section 3.5), according to This quantity is real and identical for all modes in a shape-preserving dynamic equilibrium.In our numerical implementation, we take the average over the first 8 modes (1 ≤ |m| ≤ 8) when the equilibrium criterion is satisfied.
In order to validate the solution procedure and numerical implementation, we compared our equilibria with unlimited sediment supply against equilibria in Roos et al. (2004) using their parameter settings, including their wind wave stirring mechanism which we added only for this comparison.Their equilibrium shapes and migration rates were successfully reproduced in our simulations.

Overview of Simulations
We investigate the impact of sand scarcity on sand bank equilibria under conditions that are typical for the North Sea, as listed in Table 1.Our model runs consist of two categories.First, we conduct a detailed study of two representative cases A and B. Case A represents symmetric tidal forcing by an M2-tide only.Case B represents asymmetric tidal forcing with an M2-tide and a residual current (Table 2).We vary the thickness of the sand layer D sand over eight values, being 0.05, 0.075, 0.10, 0.125, 0.15, 0.20, 0.25, and ∞.
Within the second group of runs, we vary the bottom friction parameter r, the deposition coefficient γ, and the bed slope coefficient λ within case B (Table 3).For each parameter change, the equilibrium dynamics are evaluated at D sand = 0.05, 0.10, 0.15, and ∞.Only one parameter is changed each time, while the others are kept constant.These runs are conducted to study the impact of D sand on bank equilibria relative to other morphodynamic parameters.Two runs (A: D sand = 0.20 and B-γ high : D sand = 0.15) out of 44 did not achieve an equilibrium.They evolved toward an equilibrium, but a numerical instability occurred before an equilibrium was reached.These simulations are excluded from our analysis.Furthermore, the color scheme in all figures is linked to D sand , ranging from teal (D sand = 0.05) to red (D sand = ∞).

Evolution Toward Equilibrium
The presence of a non-erodible layer changes the morphodynamic evolution toward an equilibrium from the moment that the trough exposes the non-erodible layer.The evolution of two banks with D sand = 0.10, representative for scarcity, and ∞ is compared for cases A and B (Figure 4a).The crest and trough elevation develop identically until τ ≈ 3, when the trough exposes the non-erodible layer in the runs with D sand = 0.10.The trough elevation remains constant z = 1.1, whereas the trough continues to erode in the runs without non-erodible layer.The crest elevation is not immediately affected.The crest elevation increases at the same pace for both runs until τ ≈ 6.The non-erodible layer limits the crest height from this moment onward.The runs without non-erodible layer attain a deeper trough and a higher crest than the runs with non-erodible layer (Figure 4b).The moment that the bank trough hits the non-erodible layer also triggers a jump in the change in potential energy Γ as higher order topographic modes are triggered (Figure 4c).The change in potential energy gradually decreases when the sandbank evolves toward an equilibrium.
Our results also show that neither wind wave stirring nor limited sand supply is required to obtain equilibria.Huthnance (1982aHuthnance ( , 1982b) )  conditions studied here, neither is required but both processes may affect the equilibria and the evolution toward equilibria.

Cross-Sectional Shape of Sandbanks in Equilibrium
Our model results show that under symmetric tidal conditions (case A), sandbanks in equilibrium decrease in height with increasing sediment scarcity.We present images of the cross-sectional shape in Figure 5 and show how bank properties depend on D sand in Figure 6.With increasing scarcity (decreasing D sand ), the crest elevation z crest decreases and the troughs become shallower.The depth of the bank troughs is identical to the depth of the non-erodible layer, except for D sand = ∞, which does not feature a non-erodible layer.
The bank width is found to be rather insensitive to sand scarcity in case A. The width varies between 0.32 (D sand = 0.05) and 0.38 (D sand = ∞), which is much less than the modeled variation in crest elevation.Figure 5 confirms visually that the variation in height is more profound than the variation in width with the bank slopes decreasing in steepness as D sand decreases.Sediment scarcity leads to sandbanks with a smaller cross-sectional area, which is achieved by a reduction in height rather than width under symmetrical tidal conditions.Under asymmetric tidal conditions, bank height and width are both affected by sediment scarcity.The cross-sectional shape under D sand = ∞ is both higher and wider than for any other thickness of the sand layer (Figure 5).The crosssectional profile decreases in height and width when D sand decreases from ∞ to 0.15.However, it only decreases in height when D sand decreases from 0.15 to 0.05.These observations are confirmed by the bank properties as function of D sand (Figures 6a and 6b).The crest height varies more strongly between D sand = 0.05 (where z crest = 0.85) and 0.15 ( 0.51) than between 0.15 and ∞ ( 0.39).On the other hand, the bank width decreases between D sand = ∞ (where W bank = 0.55) and 0.15 (0.37), but does not change significantly between 0.15 and 0.05 (0.34).Thus, banks initially decrease in height and width when sand is moderately scarce (D sand > 0.15) under asymmetric tidal conditions, but only decrease in height when sand is very scarce (D sand < 0.15).
Furthermore, our results show that bank asymmetry is sensitive to sand scarcity and may even be reversed under the most sediment-scarce conditions.Bank asymmetry arises under asymmetric tidal conditions.Therefore, it is only observed under case B. Sandbanks with unlimited sediment supply have a steep lee side and a gentle stoss side (Figure 5).However, the bank becomes increasingly symmetric (A bank closer to 0) as D sand decreases (Figure 6c) and is reversed for D sand ≤ 0.10, that is, the stoss side is steeper than the lee side.A visual comparison of the cross-sectional shapes shows that the slope of lee side flattens as the thickness of the sand layer decreases, whereas the slope of the stoss side varies less.

Migration
Sand scarcity increases the migration rates of sandbanks.Like bank asymmetry, migration is driven by an asymmetry in the tidal forcing and, therefore, only applies to case B. Figure 6d shows that all banks with a non-erodible layer migrate faster than the bank with an infinite sand layer.The migration rate for D sand = ∞ is 0.07, which is equivalent to 1.7 m/year under the conditions tested here.The migration rates under a φ 4 = 0 and j 2 is decreased when j 4 is increased, such that j 0 + j 2 + j 4 = 1 under all conditions.sediment-scarce conditions start at 0.13 (3.1 m/year) for D sand = 0.25.The maximum migration rate of 0.17 (4.2 m/year) is attained when D sand = 0.15.The migration rate decreases when D sand becomes thicker or thinner.These results suggest that there is an optimal thickness of the sand layer for which the migration rate attains a maximum.

Sensitivity to Morphodynamic Parameters
Our results on the cross-sectional shape (Section 4.3) and migration rate (Section 4.4) have so far been based on two cases under specific conditions.In this section, we explore whether the findings from Case B cases hold when we vary model parameters.We have selected case B for this purpose, as it covers all shape parameters and migration, whereas case A does not include bank asymmetry nor migration.Specifically, we add an M4 tidal component j 4 , as defined in Equation 18, and vary r, γ, and λ according to Table 3.For each simulation in this sensitivity analysis, we also recompute the preferred wavelength and orientation of the fastest growing mode according to Section 3.1.The variations in cross-sectional shape are presented in Figure 7.The relationships between D sand and A bank as well as c mig under varied morphodynamic conditions are shown in Figure 8.
Although variations in morphodynamic conditions change the bank shape, sand scarcity remains a key driver for bank height (Figure 7).An M4-tide decreases the crest height and flattens the slope on the stoss side.The bottom friction r, deposition γ, and bed slope λ coefficients all affect the steepness of the bank slope.An increase in r and γ, or a decrease in λ leads to steeper slopes, which results in higher and narrower banks.Conversely, a decrease in r and γ, or an increase in λ flattens bank slopes, leading to lower and wider banks.The most profound change is observed for a decrease in the deposition coefficient, which leads to significantly lower and wider banks.A decrease in D sand leads to a reduction in crest height for all tested conditions, which is consistent with results from cases A and B. Furthermore, the banks decrease in width when the sand thickness decreases from ∞ to 0.15 under all conditions, but do not further decrease when D sand decreases to 0.05.This is also consistent with results from case B.
The sensitivity analysis further demonstrates that the thickness of the sand layer affects bank asymmetry (Figure 8).The bank asymmetry parameter A bank is inversely correlated to sand scarcity under all tested morphodynamic conditions.The asymmetry always changes sign when sand is sufficiently scarce.While these observations are consistent across all runs, the value of A bank is sensitive to the morphodynamic conditions.The maximum asymmetry is attained for conditions γ high and λ low , whereas the most symmetric banks are obtained under γ low .These banks are characterized by steeper and flatter slopes respectively.
The impact of the sand layer on the migration rate is also consistent across all tested morphodynamic conditions (Figure 8).The migration rate of banks under sediment-scarce conditions (D sand ≠ ∞) is typically at least twice as fast as the migration rate under sediment-rich conditions (D sand = ∞).The migration rate attains a maximum when D sand = 0.10 or D sand = 0.15 depending on the conditions, but never at D sand = 0.05 or D sand = ∞.The migration rate itself also varies significantly across conditions.The fastest migration is found under conditions j 4 and the slowest migration under γ low .Our results show that the effect of sand scarcity on sandbank equilibria is qualitatively consistent, but that the cross-sectional shape and migration rate can vary based on exact morphodynamic conditions.

Drivers of Sandbank Migration
Our results show that the migration rate increases when the non-erodible layer is exposed and that a maximum migration rate is achieved for a specific thickness of the sand layer (Figure 9a).We explain these results by identifying the contributions to the migration rate.In a shape-preserving dynamic equilibrium, the topography satisfies We substitute the concentration and bed evolution Equations 15 and 16 into Equation 32 and integrate the result with respect to x.By evaluating the resulting function in the trough, where ∂h/∂x = 0, we obtain the integration constants z trough and 〈cu〉 trough , such that  33between the bank toes located at x 1 and x 2 (Figure 2), we find that from which c mig = (q tide + q slope )/A cs can be derived.The left-hand side of Equation 34 represents the crosssectional area of the migrating bank.The terms on the right-hand side represent the bank-aggregated contributions of the tide-driven and the gravity-driven down-slope sediment fluxes respectively.
Next, let us analyze the dependency of A cs , q tide , and q slope on D sand .We find that A cs and q tide both increase with increasing D sand but at different rates (Figure 9b).A cs depends linearly on D sand , whereas q tide varies more strongly when D sand < 0.15 than for conditions with more sediment.The bank-aggregated gravity-driven sediment flux hardly contributes to bank migration, because the flux is directed in opposite directions on both slopes.The contrasting responses of A cs and q tide to scarcity explain the maximum migration rate observed at D sand = 0.15 in case B. The low migration rate under conditions with D sand = ∞ is caused by the large crosssectional area that follows from its deep trough.At the same time, q tide is only slightly higher than under sediment-scarce conditions.

Drivers of Bank Asymmetry
Bank asymmetry is driven by an asymmetry in the tidal forcing, which leads to a tide-averaged sediment transport flux through the interaction with the topography.When bank evolution is not restricted by scarcity, our results show that asymmetry develops over time (Figure 10a, red line), and that asymmetry correlates positively with bank height (Figure 10b, red line).When sediment is scarce, the asymmetry in the topography decreases and may reverse under very scarce conditions.
We hypothesize that the reduction in asymmetry is caused by lower sandbanks and restricted entrainment when the non-erodible layer is exposed.Sediment scarcity limits bank height, which correlates with a decrease in asymmetry.Although the correlation does not physically explain the complex relation, it is consistent with stronger tide-topography interactions (e.g., residual circulations) associated with higher sandbanks.
Furthermore, the sediment transport flux is disturbed when the non-erodible layer is exposed.Our results show that the asymmetry immediately drops or stalls when the non-erodible layer is exposed, both as function of time and height (Figure 10).The bank may partially recover from the initial drop in asymmetry with greater recoveries observed for thicker sand layers.Although the non-erodible layer affects the full topography, it disturbs the stoss side more than the lee side, because the tide-averaged sediment transport flux is directed from the trough to the crest on the stoss side and from the crest to the trough on the lee side.The non-erodible layer limits entrainment in the troughs, but not on the crests.Therefore, the trough-to-crest transport on the stoss side is predominantly affected.Our model results show that this leads to steeper slopes on the stoss side and flatter slopes on the lee side.
In order to elevate our understanding of bank asymmetry, we recommend further research into the complex relation between scarcity, bank height, and asymmetry.

Limitations in Bank Shape
Our equilibrium shapes are constrained by our modeling approach in two ways.First, the domain length is fixed in our simulations at the wavelength of the fastest growing mode (Section 3.1).Therefore, the crest-to-crest distance cannot increase in our simulations, as has been observed for tidal and fluvial dunes (Nnafie et al., 2020;Porcile et al., 2020).Second, our model cannot simulate along-bank variations in topography, whereas tidal and fluvial dunes develop more irregular and isolated crestlines, including barchans, when sediment is scarce (Endo, 2016;Kleinhans et al., 2002;Nnafie et al., 2020).Findings from tidal and fluvial dunes cannot be applied to tidal sandbanks due to their distinct flow-topography interactions, but may inspire further research into bank shapes under scarcity.

Limitations of the Sediment Transport Formulation
This study considers sediment of a uniform grain size that is transported in suspension.However, field conditions can differ in three ways: (a) they typically feature a range of sediment grain sizes (Walgreen et al., 2004); (b) sediment transport occurs through a mix of suspended and bed load (Vis-Star et al., 2007); (c) sediment is only entrained above a critical velocity (Idier et al., 2009;Vis-Star et al., 2007).The sediment transport coefficient α * s can be adjusted in our model to reflect different grain sizes, but it cannot model multiple grain sizes simultaneously.Furthermore, bed load transport can be approximated using a high deposition coefficient γ*.The sensitivity analysis (Section 4.5) showed that our findings hold when γ* is raised.Nonetheless, a bed load transport formula would provide greater insight in the dynamics under bed load transport.Finally, a critical velocity for entrainment is not implemented.We recommend that these extensions are studied in future research.

Challenges in Comparing With Field Observations
Our idealized model provides generic insight into the impact of sand scarcity on sandbank dynamics under conditions that are typical for the North Sea.Our results are relevant to the Flemish Banks, Zeeland Banks, and Norfolk Banks as well as other sandbanks that experience scarcity, especially when they have receding sediment stocks (e.g., through extraction).However, a reproduction of observations is not the goal of this study for two reasons.First, a comprehensive repository of sediment budgets and characteristics of sandbank fields is currently unavailable, particularly for sediment-scarce areas.Voxel modeling (Hademenos et al., 2019) and regional studies of bank characteristics (e.g., Knaapen, 2009;Liu et al., 2007) are a first step toward this repository, but further research is required.Second, complex site-specific models are better suited toward reproducing observations, as they can include more complex topographies and site-specific conditions.In contrast, the strength of our idealized approach lies in the process-based understanding that it provides, particularly over long time scales.

Conclusions
We have developed an idealized nonlinear process-based model for the finite amplitude evolution of tidal sandbanks under sediment-scarce conditions.Our model includes tide-topography interactions (through bottom friction and the Coriolis effect), suspended sediment transport, and captures long-term nonlinear dynamics of topographies that vary in one horizontal direction only.As a novelty, a non-erodible layer limits the thickness of the sand layer.Starting from small-amplitude bed undulations, we run our model until an equilibrium bank shape is reached.
Our model results show that the thickness of the sand layer affects the cross-sectional shape of sandbanks in equilibrium.Under symmetric tidal conditions, we find that sand scarcity decreases the crest elevation and, to a lesser extent, the bank width.When the tidal conditions are asymmetric, the bank width is predominantly affected under moderately scarce conditions and the bank height is predominantly affected under very scarce conditions.Furthermore, the asymmetry gradually decreases with increasing scarcity, and may reverse under very scarce conditions.Banks with unlimited sediment supply have steep lee sides and gentle stoss sides, whereas banks under very scarce conditions have stoss sides that are steeper than the lee side.
The presence of a non-erodible layer increases migration rates of sandbanks under asymmetric tidal conditions.All runs with sediment-scarce conditions led to higher migration rates than our runs with unlimited sand availability.The migration rate attains a maximum at a specific thickness of the sand layer, which depends on the morphodynamic parameters.Thicker or thinner sand layers both decrease migration rates due to distinct responses of the bank volume and the tide-averaged sediment flux to scarcity.The bank volume decreases linearly and the tide-averaged sediment flux decreases nonlinearly in response to a thinning sand layer.
Our conclusions about the cross-sectional shape and migration rate continue to hold when morphodynamic parameters are varied.We systematically varied the tidal constituents, bottom friction, settling velocity, and the slope-driven transport parameter.The shape and migration rates responded qualitatively similar to variations in the thickness of the sand layer.A direct comparison of our model results with observations would require a comprehensive dataset of sediment budget and bank characteristics, which does not yet exist.Nevertheless, our findings provide new process-based insights into the dynamics of sandbanks in sediment-scarce areas, such as the Flemish Banks, Zeeland Banks, and Norfolk Banks.

Figure 1 .
Figure 1.The classification and response to scarcity of tidal sandbanks, tidal dunes/sandwaves, and fluvial/river dunes.Top row: the flow pattern (forcing and residual circulations) and the wavelength of each bedform.The tidal sandbank is shown in top view, because of the horizontal residual circulations and the oblique crest orientation (here Northern Hemisphere).The tidal and fluvial dunes are shown in side view motivated by their vertical flow circulations.Bottom row: the effect of scarcity on each bedform, based on Porcile et al. (2017) and Nnafie et al. (2023) for tidal dunes and on Kleinhans et al. (2002), Tuijnder et al. (2009), Dreano et al. (2010), Endo (2016), and Vah et al. (2020) for fluvial dunes.

Figure 2 .
Figure 2. Model schematization.(a) The topography without sandbank consists of a sand layer (yellow) at depth z = H* and with thickness D * sand on top of a nonerodible layer (gray) at depth z = H * ne .δ* is the thickness of the numerical buffer layer.(b) A sandbank topography h*(x, t) with shape parameters z * crest , z * trough , W * bank , and A bank = ln l * 1 / l * 2 ) and with migration rate c * mig .The ambient tidal flow approaches under an angle θ with respect to the y*-axis.The y*-axis points into the figure and is not visible on the side view.

Figure 3 .
Figure 3. Sediment limiter μ s as function of water depth, according to Equation 6.It ensures a smooth transition from unrestricted entrainment (μ s = 1) at h * = H *

Figure 4 .
Figure 4. (a) Evolution of the crest and trough elevation as function of the morphological time coordinate τ.The dashed red line represents the evolution of a bank with D sand = ∞ and the solid blue line represents a bank with D sand = 0.10.The black line denotes the undisturbed bed level.(b) Equilibrium cross-sections z eq of both banks.(c) Change in potential energy Γ over time.

Figure 5 .
Figure 5. Equilibrium cross-sections of cases A and B for D sand = ∞, 0.15, 0.10, 0.05, as denoted by the line color.The solid black line represents the undisturbed bed level.

Figure 6 .
Figure 6.Influence of layer thickness D sand on properties of the equilibrium bank shape and migration rate: (a) crest and trough elevation z crest and z trough ; (b) bank width W bank ; (c) asymmetry A bank ; and (d) migration rate c mig .These properties are defined in Section 2.1.The dotted line represents the undisturbed bed level in (a) and symmetry in (c).We did not obtain an equilibrium profile for A: D sand = 0.20.

Figure 7 .
Figure7.Sensitivity of equilibrium bank cross-sections under case B conditions to selected model parameters: j 4 added M4 tidal component; r increase/decrease bottom friction; γ increase/decrease deposition coefficient; λ increase/decrease bed slope coefficient.We did not obtain an equilibrium profile for B-γ high : D sand = 0.15.

Figure 8 .
Figure 8. Sensitivity of (a) bank asymmetry A bank and (b) migration rate c mig under case B. The colors refer to the thickness of the sand layer.We did not obtain an equilibrium profile for B-γ high : D sand = 0.15.

Figure 9 .
Figure 9. (a) Change of the migration rate c mig of the first topographic mode (|m| = 1) over morphological time coordinate τ.The crosses (×) mark the equilibria.They differ slightly from Figure 6d, because those were obtained by averaging over 8 modes (1 ≤ |m| ≤ 8).(b) Tide-driven component q tide and gravity-driven down-slope component q slope of the tide-averaged sediment flux and the cross-sectional area A cs of the sandbank as function of D sand .These three components affect the migration rate according to Equation 34.

Figure 10 .
Figure 10.(a) Change of bank asymmetry A bank over morphological time coordinate τ.(b) Relation between bank height, defined as z crest z trough , and bank asymmetry over time.The square (□), triangle (△), and cross (×) mark τ = 4, 8, and the equilibrium respectively.Both graphs (a) and (b) are smoothed over time using a 15-point moving average (Δτ ≤ 0.15) in order to reduce oscillations that follow from detecting migrating bank toes over a discrete grid.

Table 1
Overview of Physical and Numerical Model Parameters That Reflect North Sea Conditions

Table 2
Tidal Flow Components As Well As Wavelength and Orientation of the Fastest Growing Mode for Case A and B VAN VEELEN ET AL.

Table 3
List of Morphodynamic Parameters That Are Changed in the Sensitivity Analysis