Modeling Form Roughness Induced by Tidal Sand Waves

Tide‐dominated sandy shelf seas, such as the Dutch North Sea, are covered by sand waves. Yet, basin‐scale hydrodynamic models do not include any sand wave information because their grid sizes are too coarse to resolve sand waves individually. We explore the possibility of parametrizing the effects of sand waves on the larger‐scale tidal flow by means of a form roughness. Specifically, our aim is to see to what extent the flow over a sand wave field can be reproduced by that over a flat seabed with an increased effective roughness (accounting for both grain and form roughness). To do so, we use two process‐based hydrodynamic models: a second order perturbation approach, and Delft3D. Both models demonstrate that the presence of sand waves causes amplitude decrease and phase shift of the tidal flow. We explore the dependencies of form roughness on different sand wave characteristics (wavelength, height and asymmetry). Shorter and higher sand waves cause a higher form roughness, while our analysis does not reveal any dependency on sand wave asymmetry. Notably, the consideration of a tidal flow, characterized by several tidal constituents, each represented by an amplitude and a phase, results in a more complex form roughness analysis than in a fluvial setting, where the flow is unidirectional and steady. We thus obtain an amplitude‐based form roughness and a phase‐based form roughness, each yielding a different value, yet displaying the same qualitative dependencies.


Introduction
Tidal sand waves are rhythmic bedforms usually found in tide-dominated sandy shelf seas.They are characterized by wavelengths of 100-1,000 m, heights of 1-10 m and migration rates of 1-10 m/yr (van der Meijden et al., 2023;van Dijk & Kleinhans, 2005).Particularly their height and migration rate make the understanding of sand wave dynamics of practical importance.Migrating sand waves may, for instance, affect water depth in shipping routes and expose pipelines, cables or monopile foundations (Németh et al., 2003).Seas covered by sand waves vary in size, from hundreds of km 2 of the Gulf of Cádiz, Spain (Hanquiez et al., 2007), to thousands of km 2 of the Dutch Continental Shelf in the North Sea (van der Spek et al., 2022), or tens of thousands of km 2 of the northern South China Sea (Zhou et al., 2018).
Given the temporal and spatial scales involved in sand wave evolution, process-based numerical models are a suitable method to study sand wave dynamics.Hulscher (1996) explained sand wave formation as a free instability of the flat seabed subject to a symmetric tidal flow using a linear perturbation approach.Such a model has been expanded to explain the influence of, for example, tidal asymmetry (Besio et al., 2004), grain sorting (van Oyen & Blondeaux, 2009), non-erodible bed layers (Blondeaux et al., 2016), storms (Campmans et al., 2017), and benthic organisms (Damveld et al., 2019).Furthermore, Blondeaux and Vittori (2009) extended a perturbation approach up to second order to study the formation of 3D sand wave patterns.Yet, perturbation approaches rely on bedforms being low in amplitude compared to the mean water depth, thus they can only represent the sand wave formation stage.
Alternatively, Borsje et al. (2013) used the process-based nonlinear model Delft3D to study the effects that different turbulence formulations may have on sand wave formation.An advantage of Delft3D is that it allows for a more complex representation of physical processes such as turbulence and bottom friction.Delft3D also allows for the study of sand waves in the finite-amplitude regime.For example, van Gerwen et al. (2018) studied the effects of tidal asymmetry on sand wave height, but overestimated sand wave heights.Krabbendam et al. (2021) studied the evolution of a realistic sand wave bathymetry.Other nonlinear models (e.g., Campmans et al., 2018;Németh et al., 2007;van den Berg et al., 2012) can also reproduce sand wave evolution, but they overestimate sand wave heights.
The effects of sand waves on larger-scale hydrodynamics have received less attention than their morphodynamic evolution.A better understanding of form roughness induced by sand waves is needed, for example, to better model the hydrodynamics of basin-scale domains, such as in the Dutch Continental Shelf Model (DCSM) of the North Sea (Zijl et al., 2023), or the MIKE model configuration of the northern South China Sea (Zhou et al., 2018).The grid size of such models is too coarse to resolve the individual sand waves, with grid sizes of the order of kilometers for the DCSM (Grasmeijer et al., 2022).The DCSM is calibrated against observed water levels by adjusting the bed roughness (Zijl et al., 2023).Alternatively, the bed roughness in the MIKE model is adjusted based on mean water depth (Zhou et al., 2018).However, these bed roughness adjustments are not based on the local topography (which is not uniform, see Figure 1 for an example of the North Sea bathymetry), and thus cannot represent the effects of bedforms.Furthermore, the sea surface calibration of the DCSM yields significant errors in the modeled current velocities (Zijl et al., 2023), yet velocities are important for many applications (e.g., sediment transport).Hence, a bedform-based roughness parametrization would introduce a more realistic and physics-based element into basin-scale hydrodynamic models.A physics-based approach would give more reliable results, specially in scenarios where calibration is not possible because they have never taken place before (e.g., large-scale sand mining in sand wave fields).
Existing form roughness parametrizations rely on the assumption that the total roughness effectively experienced by the larger-scale flow can be written as the superposition of two contributions: grain (or background) roughness, and form roughness due to the presence of bedforms (e.g., Lefebvre & Winter, 2016;van Rijn, 1993).With this approach, one can model the flow generated over a sand wave field by replacing the bedform topography with an increased roughness over a flat bed (Figure 2).Importantly, roughness can be quantified in different ways depending on the model formulation, for example, by means of slip parameter (as in e.g., Hulscher, 1996) or a friction coefficient (as in e.g., Borsje et al., 2013).The latter is usually computed from the Chézy coefficient, the Manning coefficient or the Nikuradse roughness length.
Until now, studies on bedform-induced roughness have focused on ripples or river dunes (e.g., Lefebvre & Winter, 2016;Soulsby, 1997;van Rijn, 1993, and references therein).Based on laboratory experiments, van Rijn (1993) developed a formulation for the ripple-induced form roughness that showed good agreement with measurements from the Mississippi River.Lefebvre et al. (2014) studied form roughness induced by estuarine dunes and its temporal evolution within a tidal cycle.They found a higher form roughness when the tidal flow faced the stoss side of the sand wave, which corresponded to a larger flow separation zone over the steeper lee side of the sand wave.Their work was later expanded by Lefebvre and Winter (2016), who developed a form roughness formula for river dunes that accounted for the effects of the lee-side angle.For typical ripple dimensions, the parametrization of van Rijn (1993) yields a form roughness around 10 times larger than grain roughness.The parametrization proposed by Lefebvre and Winter (2016) yields a form roughness that depends on dune conditions, ranging from 10% of the grain roughness, up to four times the grain roughness.These values emphasize the importance of considering form roughness.
The existing form roughness parametrizations have been obtained under the assumption of steady and unidirectional flow, which is in the downstream direction for the case of river dunes, and in the direction of the instantaneous flow for the case of ripples.Sand waves, on the other hand, are located in tide-dominated regions, characterized by a tidal flow and elevation.The tide, unlike a unidirectional flow, is characterized by several tidal constituents, each expressed in terms of an amplitude and a phase.Thus, the existence of sand waves may affect the amplitude and phase of each tidal constituent differently.A comprehensive study of tidal sand wave-induced roughness demands a more refined analysis, including the representation of the tidal flow.Motivated by the above knowledge gaps, we perform a process-based modeling study on the sand wave-induced form roughness in a setting characterized by a tidal flow.The aims of this study are twofold: (a) to extend the concept of form roughness to a tidal setting and (b) to assess how sand wave-induced form roughness depends on bedform characteristics (i.e., height, wavelength, asymmetry).
To this end, we use two numerical models in this study: a second order perturbation approach analysis (PA), and Delft3D (D3D).The use of two models, each with a different physical representation of the processes involved, allows us to benefit from the strong points of each model, and to verify the obtained results.The second order perturbation approach analysis allows for the study of form roughness induced by low-amplitude sand waves on sand wave-averaged scale when forced with an oscillatory flow.Second order terms in the ratio of bedform height and mean water depth are needed since this is the lowest order at which the influence of sand waves on the spatially uniform component of the flow is felt.Delft3D simulations allow for more complex dynamics, such as tidal wave propagation, and finite-amplitude sand waves.Delft3D is used under the hydrostatic assumption.Both models are used with a fixed morphology (morphostatic) because we study the effects of bedforms on the hydrodynamics at tidal timescales (hours, days).The feedback on the morphology is beyond the scope of this study.
This study is the first to consider the effects of form roughness induced by sand waves in a tidal setting, and also the first to use a second order expansion to investigate the flow at sand wave-scale.The inclusion and comparison of two different modeling methodologies is also an innovation compared to previous sand wave studies.This paper is organized as follows.Section 2 includes a description on how form roughness is quantified, together with the descriptions of the two models used.Section 3 contains the results, which are discussed in Section 4. The final section contains the conclusions.

Modeling Approaches
We use two hydrodynamic models: a second order perturbation approach analysis and Delft3D.The use of two models, each with a different physical representation of the processes involved, allows us to verify the obtained results.This study is the first to consider a second order perturbation approach in a hydrodynamic sand wave framework.The Delft3D setting used is an extension of that of van Gerwen et al. (2018).
Both models are used in a 2DV setting, and forced with a tidal flow.The sand wave profiles used are constant in time and regular in space, the latter being a schematization of the generally non-uniform sand wave patterns observed (see e.g., van der Mark et al., 2008, and Figure 1).
The seabed is described in both models as (1) Here, H corresponds to the mean water depth, and h(x) to the sand wave topography.
Bottom friction is implemented differently in both models.The second order perturbation approach uses a linear slip condition with slip parameter s, while Delft3D uses a nonlinear friction law with friction coefficient f.A linear slip condition in the second order perturbation approach implies a simplified representation of roughness.We choose to keep that model as simplified as possible, yet being able to include the effects of sand waves.
In the following, we present the two models, the model settings and sand wave topography, and the method used to quantify form roughness in both models.The main differences between the two models are highlighted in Table 1.

Second Order Perturbation Approach (PA)
Consider a tidal flow in an offshore part of a shallow shelf sea of mean water depth H, far away from coastal boundaries.We make use of a Cartesian 2DV coordinate system (x, z), with horizontal coordinate x, and the vertical coordinate z pointing upward (Figure 3a).The sea surface elevation is denoted by z = ζ(x, t), with z = 0 the undisturbed water level.As already explained in Equation 1, the seabed is denoted by z = H + h(x), where h(x) represents the sand wave topography (with zero spatial average).The horizontal and vertical flow velocities are denoted by u(x, z, t) and w(x, z, t), respectively.
Note.Symbol definitions are given in the main text.
Unsteady flow is described by the shallow water equations, that is Here, g is the gravitational acceleration, and the eddy viscosity A v is assumed constant and uniform.The tidal forcing is decomposed into a sand wave-scale component (first term in the right-hand side of Equation 2a), and a tidal wave-scale component (F in Equation 2a).The latter describes a barotropic pressure gradient forcing representing the spatially invariant part of the surface elevation gradient.Thus, F is a uniform and oscillatory signal, , with pω and Fp the frequency and (complex) amplitude of each tidal component, chosen such as to achieve a prescribed depth-averaged flow signal over a flat bed.We truncate the forcing at P tidal harmonics.
At the free surface, we impose a dynamic and a kinematic boundary condition.The former imposes zero shear stress, while the latter imposes the particles to follow the water surface.A scaling procedure as performed by Campmans et al. (2017) (not shown here) reveals that the rigid lid approximation may be applied, that is, the sea surface boundary conditions may be imposed at z = 0 instead of z = ζ(x, t).Both conditions thus read Similarly, the partial slip boundary condition and a kinematic condition imposing the particles to follow the bed profile are imposed at the bed, that is, Here, ρ is the density of water, and τ bed is the bed shear stress, which is parametrized with the slip parameter s.We impose horizontal spatially periodic boundary conditions.
Extending the linear stability analysis of Hulscher (1996), the unknowns are expanded around a basic flow, that is, that corresponding to a flat bed, up to second order, and higher order terms (h.o.t.) are neglected.Here, ϕ represents any property of the system (e.g., flow velocity, sea surface elevation).The basic and first order flow signals are denoted by ϕ 0 , ϕ 1 , respectively.Second order terms are needed since this is the lowest order at which the influence of sand waves on the spatially uniform flow is felt.Note that higher order terms are assumed to be much smaller than ϕ 0 , that is, The sand wave topography is described by Here, h sw is the sand wave height, and k sw is the topographic wavenumber.We impose horizontal periodicity by taking the length of the domain as a multiple of the sand wave wavelength λ sw = 2π/k sw .
The expansion in Equation 5 is applied to Equations 2-4.Equations for the basic and first order flow are given in Appendix A. The second order flow is only solved for its spatially uniform component, that is, the sand waveaveraged flow.We use the notation 〈⋅〉 = 1 λ sw ∫ λ sw 0 ⋅ dx to denote sand wave-averaged quantities.The resulting second order, sand wave-averaged momentum and mass conservation equations read.
From the right-hand side of Equation 7a, we infer that the first order advective term acts as a forcing in the second order flow equations.
The second order boundary conditions at the free surface read Note that due to mass conservation (Equation 7b) and the kinematic boundary condition at z = 0, it follows that 〈w 2 〉 ≡ 0 for any z.Thus, the second order boundary conditions at the seabed read which have been expanded around the flat bed, such that they are applied at z = H instead of z = H + h 1 (x).
Note that Equation 9b is automatically satisfied due to mass conservation, since ⟨u 1 , where in the first equality we have applied the first order continuity equation, that is, ∂u 1 ∂x + ∂w 1 ∂z = 0.
Given the temporal dependency of the forcing, we expand the sand wave-averaged part of the second order flow solution as Analogous to the basic and first order flows (Appendix A), Equation 7is solved for the truncated Fourier components (P = 2) using the Fourth Order Runge-Kutta method, and a central differences scheme for the zderivatives of the forcing terms.We performed a sensitivity analysis on the order at which the solution is truncated (not shown here), for an M2 forcing, and the differences in form roughness values obtained with P = 2 and P = 4 are less than 0.5%.
We use the depth-averaged, sand wave-averaged flow 〈U PA 〉 to compute a form roughness (see Section 2.4 for more details), with 〈U 0 〉 and 〈U 2 〉 the depth-averaged, sand wave-averaged contributions of the basic and second order flow, respectively, defined as

Delft3D
We use Delft3D in a 2DV configuration based on that of van Gerwen et al. (2018).We solve the 2DV shallow water equations in the (x, σ)-coordinate system, where σ = 0, 1 correspond to the seabed and sea surface, respectively.The seabed, which is imposed to be static, consists of a central region of length L sw = 2.1 km with sand waves, surrounded by a flat seabed of length L flat = 18.54 km at each side (Figure 3b).The total domain length is thus L = 2L flat + L sw = 39.18 km.The transition from flat bed to full amplitude sand waves follows a sinusoidal envelope shape of horizontal extent L env = 100 m (within L sw ).The horizontal grid spacing is nonuniform.It takes a maximum value of Δx max = 1,542 m at x = 0, L, and gradually decreases to Δx min = 2 m, which holds everywhere in the sand wave region (Figure 3d).The thickness of each σ-layer with respect to the local water depth gradually increases from 5 • 10 4 at the bed, to 3.28 • 10 2 at the sea surface.
The model has two open boundaries, at x = 0, L, where the tidal wave is forced using Riemann invariants, R, which prescribe a linear combination of the horizontal depth-averaged velocity U and the sea surface elevation ζ, Imposing a Riemann invariant minimizes the reflection of the tidal wave at the boundaries.Here, R + = 1.2 m/s is imposed at x = 0, and R = 0 m/s at x = L, such that we obtain a tidal wave propagating in the positive x-direction.
The model yields a local approximation of a Kelvin wave, thus being representative of a region far from the coast, and more realistic than the standing wave imposed in van Gerwen et al. (2018).Given that the model is 2DV, Coriolis effects and the amplitude decay away from the coast are neglected, which is justified by the length of the domain being much smaller than the Rossby radius of deformation.In the remainder of the text we will use the terms up-wave and down-wave to describe horizontal positions with a lower or higher x-value than a reference, respectively.
The bed shear stress is computed at the seabed using a nonlinear friction law, Journal of Geophysical Research: Earth Surface 10.1029/2023JF007610 PORTOS-AMILL ET AL.
Here, f = g/C 2 is the friction coefficient, commonly expressed in terms of the gravitational acceleration g and the Chézy coefficient C. Note that roughness may be expressed in terms of many other different parameters (e.g., roughness length, Chézy coefficient, Manning coefficient).We choose f as it provides the most accessible formulation of the seabed boundary condition (Equation 14).

Model Settings and Sand Wave Topography
Both models are forced with an M2 tidal signal of frequency ω = 1.4 • 10 4 rad/s.The amplitudes of the forcings in both models are set such that the resulting flows are comparable for a flat bed situation.Delft3D is forced with the Riemann boundary conditions, which generate a tidal wave with depth-averaged flow amplitude of 0.59 m/s at the center point of the domain.We choose the amplitude F±2 of the forcing in the second order perturbation approach such that the depth-averaged basic flow also has an amplitude of 0.59 m/s, yielding F±2 = 5 ⋅ 10 5 m/s 2 .
Both models have a mean water depth of H = 30 m.The grain roughness values are taken equal to previous studies.The grain-induced slip parameter s g for the second order perturbation approach is 1.0 • 10 2 m/s (as in Campmans et al., 2017), and the grain-induced friction coefficient f g in Delft3D is 2.3 • 10 3 (as in Borsje et al., 2013).The eddy viscosity A v in the second order perturbation approach takes the value 0.04 m 2 /s (as in Campmans et al., 2017).
The default sand wave shape used is symmetric, with h sw = 4 m and λ sw = 350 m (black line in Figure 4).With both models, we perform simulations with different sand wave height (h sw ), and wavelength (λ sw ).The wavelengths chosen for the Delft3D simulations are such that there is always an integer number of sand waves, N sw ≥ 2, within the region of interest.Additionally, the analysis performed with Delft3D simulations allows to study the influence of sand wave asymmetry.The sand wave topographic profile is expressed in terms of the sand wave height h sw , wavenumber k sw , and asymmetry A sw as Here, the parameter A sw represents the difference between stoss-side and lee-side lengths (L s and L l , respectively) relative to the sand wave wavelength, The stoss side is the sand wave side with a milder slope, and sand waves are always flood oriented.Note that flow separation is not expected given the lee-side angles considered.The most asymmetric sand wave profile considered (A sw = 0.7) corresponds to a lee-side angle of 3°, and flow separation is expected only for angles higher than 10° (Lefebvre et al., 2014).

Quantification of Form Roughness
Since both models have a different representation of the tidal flow (Table 1), each requires a different method to quantify form roughness.However, the overall principle is identical.
Bed roughness is incorporated in the second order perturbation approach (PA) by increasing the slip parameter s used for the dynamic boundary condition at the seabed (Equation 4).The sand wave-averaged, depthaveraged flow over a sand wave field (with the grain-induced slip parameter s g ) is compared with that over a flat bed (with an increased "effective" slip parameter, s eff > s g ).The additional value needed for the slip parameter such that both flow signals have the same effective amplitude or M2 phase expresses the form roughness, that is, Importantly, this leads to two form-induced slip parameter values, amplitude-based and phase-based, which are not necessarily identical.The effective amplitude is defined as half of the difference between maximum and minimum values of the depth-averaged flow velocity.The M2 phase is obtained from the Fourier transform of the depth-averaged flow velocity.This will be further addressed in Section 3.
For the Delft3D simulations, we follow a similar approach, but we do not consider the sand wave-averaged signal because of the propagating nature of the tidal wave.Instead, we consider the depth-averaged flow at a position down-wave of the sand wave field (red dot in Figure 3b, x = 20.75 km).We compare it with the depth-averaged flow at the same x-position obtained with a simulation with a flat bed and an increased "effective" friction within the region of interest (L flat < x < L flat + L sw , Figure 3c).A similar analysis is also done with the sea surface elevation instead of the depth-averaged flow.As before, the added friction needed for both signals (depth-averaged flow or sea surface elevation) to have the same effective amplitude or M2 phase is the form roughness: This produces four roughness values: amplitude-based and phase-based, for both depth-averaged flow and sea surface elevation.The spatial ramp-up toward an increased friction is done using the same sinusoidal shape as for the envelope of the sand wave field (Figure 3c).
Because of the different roughness parametrizations, form roughness expressed in terms of a slip parameter (PA) and friction coefficient (Delft3D) should only be compared qualitatively, not quantitatively.

Effects of Sand Waves on Tidal Flow
We first examine the effects that a sand wave field exerts on tidal flow by studying the temporal evolution of the depth-averaged flow obtained with both models.
The presence of sand waves results in a decrease in amplitude and a phase shift of the depth-averaged flow, as shown in Figure 5.The top two panels of Figure 5 refer to the results obtained with the second order perturbation approach, and the bottom five to the Delft3D results.To mimic the effective roughness experienced by the presence of sand waves (black lines in Figures 5b and 5g), a flat bed with an increased roughness is imposed (gray dashed lines in Figures 5b and 5g).The example shown is obtained following the procedure outlined in Section 2.4: by imposing an effective roughness such that both signals have the same velocity amplitude.Note, however, that even if both signals have the same amplitude, they differ in phase (visible in Figure 5b).
Moreover, the flow obtained from the second order perturbation approach is purely oscillatory, while Deflt3D is able to represent the tidal wave propagation through the domain.This propagation can be inferred from the slope of the lines with equal flow magnitude in Figure 5c, albeit that is hardly visible since the scale of the domain (40 km) is small compared to the tidal wavelength (∼770 km).Additionally, results show that the flow adapts to the topography (Figure 5d), that is, it has a larger (lower) amplitude over a sand wave crest (trough).Contrarily, the simulation with a flat seabed and an increased effective roughness (Figure 5e) only shows tidal wave propagation.

Dependencies on Sand Wave Characteristics
Following the procedure in Section 2.4, both models allow for the study of the influences of sand wave characteristics on form roughness.For the second order perturbation approach analysis (Figure 6) this is done for sand wave height and wavelength.Form roughness, computed from the depth-averaged and spatially averaged flow, increases for higher sand waves and decreases for longer sand waves.Highest sand waves (h sw = 8 m) yield a form-induced slip parameter of s f = 0.7 • 10 3 1.8 • 10 3 m/s, for the amplitude-based and phase-based criteria, respectively.This corresponds to 7 18% of the grain roughness applied.Shortest sand waves yield s f = 3 • 10 4 5.5 • 10 4 m/s, corresponding to 3 5.5% of the grain roughness.The default sand wave characteristics (h sw = 4 m, λ sw = 350 m) yield s f = 2 • 10 4 4 • 10 4 m/s, which corresponds to 2 4% of the grain roughness.Both criteria (amplitude-and phase-based) yield the same qualitative dependencies on h sw and λ sw .Note that a quadratic dependency of form roughness on sand wave height is expected given that sand wave height only enters Equation 11 by means of the second order factor (U 2 ), which is effectively proportional to h 2 1 .The effects of sand wave asymmetry cannot be studied since the second order perturbation approach analysis does not distinguish between bed harmonics with a different relative phase.
In contrast, Delft3D allows for the computation of form roughness from the depth-averaged flow velocity and the sea surface elevation (black and green lines in Figure 7, respectively).The form-induced friction coefficient f f obtained is highly dependent on whether the criterion applied is phase-based or amplitude-based.Both phasebased criteria (depth-averaged flow and sea surface elevation) give a form-induced friction coefficient that is up to 20 times higher than the one obtained with the amplitude-based criteria.Differences between the two amplitude-based or the two phase-based criteria are always lower than 20%.All four criteria yield the same qualitative behavior as for the second order perturbation approach, that is, increasing form roughness for increasing sand wave height, and decreasing form roughness for increasing sand wave wavelength.In contrast to PA results, the phase-based form roughness obtained with Delft3D simulations is of the same order of magnitude as the grain roughness (thin horizontal lines in Figure 7).Furthermore, unlike PA, Delft3D allows for the inclusion of sand wave asymmetry in our study (Figure 7c).Notably, the obtained form roughness shows no dependency on the sand wave asymmetry.The default sand wave profile yields a form-induced friction coefficient of f f = 1.9 • 10 3 for the phase-based criterion, which corresponds to 82% of the grain-induced friction coefficient.

Discussion
A major innovation of our study is the extension of the concept of form roughness to a tidal setting.Previous form roughness studies only considered a steady and unidirectional flow, mainly characterized by a single parameter (e.g., flow magnitude).Differently, the tide is characterized by the amplitudes and phases of the different tidal constituents.In this study, we narrowed this down to an effective amplitude and the phase of the M2 constituent.As a result, we found four criteria to determine form roughness induced by tidal sand waves, amplitude-based and phase-based for either depth-averaged flow or sea surface elevation.The fact that these criteria all produce different form roughness values demonstrates that flow over a sand wave field cannot be fully mimicked by an increased roughness over a flat bed.The analysis of form roughness in a tidal setting instead of a fluvial setting thus brings in new considerations, which underlines the added complexity and novelty of our work.

Comparison With Existing Form Roughness Parametrizations
The study on river dunes of Lefebvre and Winter (2016) computes the associated form roughness from model simulations.Specifically, they compute the form friction factor from the average water level slope, thus yielding a unique criterion for form roughness quantification.The resulting dune roughness is between 10% and 80% of the total friction, and it increases with steeper and higher dunes, and decreases for longer dunes.The dependencies on sand wave height and length agree with our results.Furthermore, the order of magnitude of the obtained form roughness with respect to the grain roughness agrees with our phase-based Delft3D results (Figure 7), yet our second order perturbation approach results are one order of magnitude lower.
Alternatively, the ripple study by van Rijn (1993) formulates form roughness parametrizations for ripples and sand waves based on flume and field data.As for ripples, the resulting form roughness increases with ripple height in situations where currents dominate over wind waves.Ripple-induced form roughness decreases for longer ripples.These dependencies on ripple shape agree with our results.
Notably, form roughness studies generally state that form roughness only exists due to flow separation (Lefebvre et al., 2014;van Rijn, 1984), yet our results demonstrate that sand waves induce form roughness even under conditions where flow separation does not occur.Particularly, van Rijn (1993) argues that form roughness induced by symmetrical sand waves is zero because they can be seen as background topography, thus not inducing flow separation.The different approaches between his study and our study explain the (seemingly) contradictory outcomes.We focused on sand wave field scale, while van Rijn (1993) focused on a local scale, thus within the sand wave field.On smaller spatial scales, a sand wave profile can indeed be seen as a varying background water depth.This highlights the importance of selecting appropriate spatial scales depending on the specific knowledge gap and envisaged application of the results.

Lack of Dependency on Sand Wave Asymmetry
Our Delft3D results show no dependency of form roughness on the degree of sand wave asymmetry (Figure 7c), which disagrees with the increased form roughness obtained for increasing lee-side angles by Lefebvre and Winter (2016) in a fluvial setting.One of the possible reasons behind this discrepancy is the symmetric tidal forcing imposed in our model.The existence of asymmetric sand waves has been explained by the presence of residual currents (e.g., Besio et al., 2004).Thus, an asymmetric tidal forcing would better represent a situation with asymmetric sand waves, and could lead to asymmetry dependencies.However, it is unclear which residual current value is needed for which level of asymmetry, and allowing for a residual current would enlarge the parameter space, being thus out of the scope of this study.Future studies should examine the role of a residual current on the definition of form roughness in a tidal setting.
We examined the turbulent kinetic energy (TKE) over a symmetric and an asymmetric sand wave field at peak ebb, peak flood, and averaged over the tidal cycle (Figure 8).The TKE pattern over the symmetric sand wave field during peak flood (b1) is the mirror image of that during peak ebb (c1).This results in a symmetrical pattern of the tide-averaged TKE with respect to the sand wave crest (a1).In contrast, the TKE pattern over an asymmetric sand wave field shows a wake of high TKE on the lee side during peak flood (b2), and shows a weaker signal during peak ebb (c2).This results in an asymmetrical tide-averaged TKE pattern (a2).
Remarkably, the TKE averaged over the sand wave field during peak flood is higher for the asymmetric sand wave field than for the symmetric one.This is caused by the wake generated at the lee side of the asymmetric sand wave and results in a lower tidal flow during flood over the asymmetric sand wave field than over the symmetric one.Furthermore, the TKE averaged over the sand wave field during peak ebb is lower for the asymmetric sand wave field than for the symmetric one because of the milder stoss slope.This results in a higher tidal flow (in absolute value) during ebb over the asymmetric sand wave field than over the symmetric one.These two effects combined result in a depth-averaged flow of the same amplitude for both sand wave fields, thus yielding the same form roughness, even if the differences in sand wave asymmetry are felt by the flow.Lefebvre et al. (2014) obtained similar dependencies of the TKE to the tidal phase for an asymmetric sand wave field.Note that the hydrostatic assumption of Delft3D still holds for the most asymmetric sand wave shape considered, with lee-side angles of 3°for A sw = 0.7.

Second Order Perturbation Approach Versus Delft3D
In our study, we used two different models to quantify form roughness induced by sand waves on the tidal flow.Both approaches have different representations of the physical processes involved, the second order perturbation approach (PA) being strongly schematized, and Delft3D including more complex physical processes.For example, roughness is parameterized by means of a linear slip condition at the seabed, while Delft3D imposes a nonlinear friction law.Similarly, the eddy viscosity is assumed to be constant and uniform in the second order perturbation approach, while Delft3D uses the k ϵ turbulence closure.Despite these differences, both models provide qualitatively the same results, which serves as a verification of the dependencies of form roughness on the sand wave height and wavelength.Indeed, the second order perturbation approach, with a minimum of physical processes involved, is able to reproduce the dependencies obtained with the complex Delft3D model.The difference in order of magnitude between the phase-based results shows that it is the more complex processes, included in Delft3D, that cause higher roughness values.Furthermore, the second order perturbation approach is constrained to infinitesimally small amplitudes, while Delft3D allows for finiteamplitude sand waves.Yet, the second order perturbation approach complements Delft3D by being much more computational efficient, thus allowing for a broader parameter space inspection (Figures 6c and 6d).

Practical Relevance
The main motivation behind this study was to assess the importance of form roughness induced by sand waves at a sand wave field scale.Results show that form roughness can be of the same order of magnitude as the grain roughness.Thus, it should be incorporated into models that have too coarse grid sizes to explicitly resolve for sand waves.This way, the effects of sand waves would be implicitly accounted for.
Currently, basin-scale hydrodynamic models such as the DCSM for the North Sea (Zijl et al., 2023) and the MIKE model on the northern South China Sea (Zhou et al., 2018) do not use bedform-related information to assign a value to the bed roughness.The former computes the bed roughness based on a calibration of the sea surface elevation, which results in patches of different bed roughness along the domain, implying abrupt (and unrealistic) changes in the imposed bed roughness.The latter, being a 2DH model, imposes a bed roughness that depends solely on the mean water depth.Bed friction values imposed in the DCSM range within f = 0.7 • 10 3 4 • 10 3 (Zijl et al., 2023), while our results suggest a friction coefficient ranging within f = 2.3 • 10 3 10 • 10 3 , depending on the sand wave characteristics and on the criterion applied.Thus, our model results may serve to include bedform-related information in the calibration process of these models.This would result in roughness having a greater physical basis than in current practice.
The here highlighted existence of several criteria to quantify form roughness is a sign of the complexity of a tidal setting, but it also allows for the criterion selection that best fits the application.For example, the DCSM is calibrated against water level observations.One could thus use a sea surface elevation-based criterion to implement form roughness in the DCSM in order to match the already applied criterion.On the other hand, one could apply a linear combination of the roughness value obtained from the calibration against sea surface elevation observations and the form roughness obtained from a velocity-based criterion, such that the model settings account for both properties of the system.Future research should investigate how the different criteria affect larger-scale hydrodynamics and if the consideration of form roughness reduces the error in modeled flow velocities encountered in the DCSM.
The relative errors obtained on the other three criteria when applying a certain criterion to optimize form roughness are shown in Table 2.When using the amplitude-based criteria, the error in the other criteria is of the same order of magnitude (10 3 %).On the other hand, the phase-based criteria minimize the error in both the phase of U and ζ (10 4 %), while inducing a larger error for the amplitude based criteria (∼10 1 %).Future studies should focus on how these errors translate to larger-scale hydrodynamics.Note that different settings (e.g., mean water depth, tidal conditions) are likely to produce different form roughness values, while here we have only focused on sand wave-related dependencies.

Conclusion
The aims of this study were (a) to extend the concept of form roughness to a tidal setting, and (b) to assess how sand wave-induced form roughness depends on bedform characteristics.Form roughness had previously been considered in river settings only, the consideration of sand wave-induced form roughness in a tidal setting is an innovation of this work with respect to previous studies.
We developed a second order perturbation approach analysis to study the sand wave-averaged flow over a sand wave field.This was complemented by Delft3D simulations.Both models were run in a morphostatic setting.The combination of both modeling methodologies is another novelty of the present study.Furthermore, the use of two models serves as a verification of the obtained results, given that both yield qualitatively similar results.
Studying form roughness in a tidal setting requires a more complex analysis than in fluvial settings, where the flow is steady and unidirectional.The different characteristics of the tidal flow (amplitude and phase of each tidal constituent) are altered by the presence of a sand wave field, all in a different manner.This leads to several criteria defining form roughness induced by tidal sand waves.The distinct form roughness values obtained with each criterion demonstrate that the flow over a sand wave field cannot be fully mimicked by an increased roughness over a flat seabed.
Shorter and higher sand waves lead to an increased form roughness.Obtained results agree qualitatively with existing form roughness parametrizations for non-oscillatory settings (river dunes and ripples), except for the dependency on sand wave asymmetry.The obtained form roughness can be the same order of magnitude as the grain roughness, proving the importance of taking into account sand wave-induced roughness in larger-scale hydrodynamic models, where sand waves cannot be resolved due to grid resolution limitations.Future studies could apply these findings to specific regions, for example, the Netherlands Continental Shelf, to quantify the effects of sand wave-induced form roughness on larger-scale hydrodynamics.because in the basic state h 0 ≡ 0 and thus there is homogeneity in the x-direction.This implies that w 0 ≡ 0 for all (x, z).The momentum equation for the basic flow is then with the boundary conditions Here, the symbol {•} p denotes the pth Fourier component of the bracketed term, which involves a truncated convolution of the product inside.Such a convolution produces a cascading effect on the modes needed to solve the system of equations.Yet, we set a maximum P = 2, above which no further modes are considered.We solve truncated Fourier components (P = 2) using the Fourth Order Runge-Kutta method, and a central differences scheme for the z-derivatives of the forcing terms.

Figure 1 .
Figure 1.Bathymetric chart of the Netherlands Continental Shelf, highlighting three regions with different bed topography: (a, b) sand waves with different characteristics, (c) no sand waves.Note the different scale of the color bars.Data from the Netherlands Hydrographic Service.

Figure 2 .
Figure 2. Schematic representation of sand wave-induced form roughness: (a) the roughness that the flow experiences over a flat seabed is purely due to grain roughness.Sand waves add additional roughness to the seabed, which can be modeled either by (b) resolving the sand waves or (c) imposing an increased effective roughness over a flat seabed.

Figure 3 .
Figure 3. Main features of both models: (a) schematic of the second order perturbation approach, (b) two example bathymetries with different wavelengths used in Delft3D, (c) example of increased friction profile used in Delft3D simulations with a flat bed ( f eff = f g + f f , Section 2.4), (d) horizontal grid spacing Δx used in Delft3D.Note the different x axes in panels (b, c, d).Shaded areas in panels (b, c) denote the envelope regions of length L env where the amplitude of the bathymetry and the friction coefficient change gradually.The red dot in panel (b) is the x-position used in the Delft3D analysis to quantify form roughness (Section 2.4).

Figure 4 .
Figure 4. Examples of sand wave profiles with different degrees of asymmetry, with the default sand wave height (h sw = 4 m)and wavelength (λ sw = 350 m).Note that the second order perturbation approach only allows for the symmetric topography (A sw = 0).

Figure 5 .
Figure 5. Effects of sand waves on the depth-averaged flow (a, b: PA; c-g: D3D): (a) Results obtained with the second order perturbation approach, comparing the situation with a flat seabed and background roughness, sand wave field with background roughness, and flat seabed and effective roughness.The orange box indicates the zoomed-in region in panel (b).(c) Propagating tidal wave over a sand wave field obtained with Delft3D.The red box indicates the zoomed-in region in panel (d).(e) Same as panel (d), but over a flat bed with an increased effective roughness.Dashed vertical lines in panels (c, d, e) indicate the x-position for which results in panels (f, g) are shown (x = 20.75 km).(f) Same as panel (a), but performed with Delft3D.The blue box indicates the zoomed-in region in panel (g).The effective roughness values are those obtained with the amplitude-based criterion (s eff = 1.17 ⋅ 10 2 m/s, f eff = 3 • 10 3 , corresponding to 117% and 130% of the grain roughness, respectively).An extreme sand wave shape is used such that results are visible in the image (h sw = 8 m, λ sw = 100 m, symmetric).

Figure 6 .
Figure 6.Form roughness obtained from the second order perturbation approach (s f ): (a) in terms of h sw , (b) in terms of λ sw .Panels (c, d) show the dependencies on both h sw and λ sw with the amplitude-and the phase-based criteria, respectively.White lines in panels (c, d) mark the h sw and λ sw values selected for panels (a, b), h sw = 4 m, λ sw = 350 m.

Figure 7 .
Figure 7. Form roughness obtained from Delft3D ( f f ) in terms of (a) h sw , (b) λ sw , (c) A sw .Thin black line corresponds to the grain roughness ( f g ).Note the different vertical axes.

Figure 8 .
Figure 8. Turbulent kinetic energy (TKE) obtained with Delft3D simulations over a symmetric sand wave field (1) and over an asymmetric sand wave field with A sw = 0.7 (2).Panels (a) show the tide-averaged signal, while panels (b, c) show the peak flood or peak ebb signals, respectively.Black lines are meant for visualization purposes.In panels (a) they mark the level with TKE = 1.45 • 10 3 m 2 /s 2 , and in panels (b, c) they mark the level with TKE = 2.9 • 10 3 m 2 /s 2 (also marked on the color bars).Note the different scales of the color bars in panels (a-c).Arrows in panels (b, c) depict the direction of the flow.
The equations defining the basic flow are obtained by applying the expansion in Equation 5 into Equations 2-4 and keeping only the 0th order terms.The mass conservation equation and kinematic boundary conditions reduce to ∂w 0 ∂z = 0, w 0 = 0 at z = H, 0, (A1)

Table 2
Form Roughness Values and Relative Error (in %) Obtained for the Different Criteria (Rows) When Applying a Certain Criterion to Optimize Form Roughness (Columns)Note.Results obtained with the Delft3D simulations for the default sand wave profile.