The importance of time-varying, non-tidal currents in modelling in-situ sand wave dynamics

Sand waves are found on shallow, sandy seabeds throughout the world and their dynamics may pose an imminent threat to offshore construction. Therefore, there is a pressing need to understand bed level dynamics in sand wave areas. These bed level dynamics lead to variations in sand wave shape and migration rate over time. However, these variations cannot be explained with the present-day process-based sand wave models, which all include a purely periodic tidal forcing. To explain these fluctuations a more intricate description of the hydrodynamics is necessary. The aim of this study is to explore the importance of time-varying, non-tidal currents for sand wave dynamics in the North Sea. We adopted the three-dimensional Delft3D-Flexible Mesh model, and were able to reconstruct time-varying, non-tidal currents on top of the periodic tidal forcing, while significantly reducing computation times. The simulated currents and water levels showed a good agreement with in-situ measurements. Compared to the situation with only tidal forcing, the simulated sedimentation and erosion rates were amplified up to 15 times due to time-varying, non-tidal currents. Additionally, periods of net erosion were found at locations in the sand wave transect where tidally forced models only showed net-sedimentation. It is therefore important to consider time-varying, non-tidal currents when predicting future sand wave dynamics in the field.


Introduction
Due to ambitious green energy goals, the pressure on the offshore area is increasing at an unprecedented pace.Through the Green Energy Deal, the European Union is planning to increase its offshore wind capacity to 300 GW by 2050.This entails a twelve folding of the current installed capacity of approximately 25 GW (European Commission, 2022).Especially the North Sea, with its suitable wind climate and large shallow areas, is to play a major role in this capacity increase.
At the same time, the bed of the North Sea is covered with large-scale rhythmic bed forms (Damen et al., 2018).Especially sand waves may interfere with offshore activities, due to their size and dynamic character (Morelissen et al., 2003).Sand waves can be found at many places around the world, such as the Taiwan Strait (Bao et al., 2014), San Francisco Bay (Sterlini et al., 2009), the Mid-Atlantic Bight (Bokuniewicz et al., 1977) and the Banks Strait near Australia (Auguste et al., 2021).In the Dutch North Sea these rhythmic bedforms have wavelengths of 100-1 000 m, wave heights of 1-10 m (Damen et al., 2018) and migrate with up to 20 m per year (Van der Meijden et al., 2023).
Sand wave migration or change in shape can expose subsea cables and pipelines and reduce navigation depths (Németh et al., 2003).Many activities in offshore areas, such as the construction and maintenance of windfarms, thus require bed level predictions on lifecycle timescales (i.e. decades).Using data-driven methods, historic sand wave migration rates can be determined, see for example Van der Meijden et al. (2023).By extrapolating these migration rates into the future, a range of possible bed levels can be estimated for engineering purposes.The use of process-based numerical models could improve the understanding of sand wave behavior and accuracy of predictions, and help quantifying uncertainties over time.
Idealized numerical modelling studies have led to a broad understanding of the processes underlying sand wave dynamics.Sand wave formation is caused by a complex interaction between (tidal) currents, bed topography and sediment transport (Hulscher, 1996).Initial perturbations of the sandy seabed grow into sand waves due to the net convergence of sediment transport at the crest of the perturbation.Due to dampening effects, such as gravitational forces, the sand wave height will eventually reach an equilibrium (Van Gerwen et al., 2018).Sand wave migration can be induced by a residual current (Németh et al., 2002;Besio et al., 2003) and the superposition of the M4 tidal component (Besio et al., 2004).Due to the latter, migration can occur in either flood or ebb direction, depending on the relative tidal phases of the M2 and M4 tidal current.Campmans et al. (2018) studied the effect of surface gravity waves and constant wind forcing and found that both processes cause a reduction of equilibrium sand wave height and an increase of sand wave migration.Although these idealized numerical models have provided valuable insights, their ability to predict in-situ sand wave dynamics has proven to be limited.Through calibration of an idealized numerical model, Campmans et al. (2022) were able to improve the match between modelled and observed sand wave lengths.However, even after calibration the equilibrium sand wave height was overestimated significantly.
The application of process-based models to in-situ sand wave cases is still in its infancy.Leenders et al. (2021) simulated 3D sediment transport patterns, using Delft3D-4, to investigate the influence of underlying sand banks on sand wave dynamics.However, due to computational limitations the model duration was limited to a couple of days.Krabbendam et al. (2021) applied this same model for predicting bed level developments in sand wave areas.Through calibration of the morphological parameters, a good match between the modelled and measured sand wave migration was found.However, the shape of the sand waves was altered during the simulation.The model simulation showed lowering of both the sand wave crests and troughs and rounding of steep bed gradients, in contrast to what was measured in the field.Consequently, the simulation results indicated that there are still processes which are either missing or need to be refined.
In previous modelling studies, it was assumed that the dynamics of sand waves are the result of the main tidal components (M2, S2, M4), combined with a constant residual current (Z0).The hydrodynamic forcing is thus purely periodic and repetitive.However, this type of forcing is unable to explain changes in sand wave shape and migration rate over time, which are often observed in the North Sea (see for example Fig. 1E and F and Menninga (2012)).In the shallow sea environment, where we find sand waves, several types of time-varying non-tidal influences on current velocity and water level are present.These include, but are not limited to, wind-driven currents, storm surges, currents driven by differences in atmospheric pressure, and density driven currents.Although the importance of time-varying non-tidal currents is recognized for shallower bed features such as sand ridges and sand banks (Van Rijn, 1998;Dibajnia et al., 2011), the exact effect of these currents on sand wave dynamics remains unexplored.These time-varying non-tidal currents may significantly influence sand wave dynamics and its full effect is not captured by the constant residual current applied in state-of-the-art sand wave models.The importance of these time-varying non-tidal currents is affirmed by observations of sand waves on the Taiwan Shoal by Bao et al. (2020).They discovered a significant influence of a passing tropical storm, which caused quick migration of the sand waves, combined with a significant reduction in sand wave height.We hypothesize, that these time-varying non-tidal currents could also explain changes in sand wave dynamics over time in other marine environments.Therefore, the aim of this paper is to assess the impact of realistic fluctuations of non-tidal currents on in-situ sand wave dynamics.To reproduce these hydrodynamic currents and their effect on sediment transport rates the process-based numerical model Delft3D FM is applied.

North Sea cases
To allow for a wider understanding of the importance of timevarying, non-tidal currents on sand wave dynamics two contrasting locations in the North Sea are selected.These sites are chosen such that they are representative of the range of conditions within dynamic North Sea sand wave areas both in terms of hydrodynamics and sand wave characteristics and dynamics (see Damen et al. (2018) and Van der Meijden et al. (2023)).In this section first the study sites are discussed, followed by the data availability for each site.

Study sites
To be able to validate the hydrodynamics within the model, the first site was chosen where Acoustic Doppler Current Profiler (ADCP) measurements are available.A site within the Hollandse Kust Zuid (HKZ) wind farm zone was found (see Fig. 1A and C).At this site sand waves are observed with a height of around 3 m and a wavelength between and 700 m, at an average depth of approximately 23 m.The sand waves migrate in flood direction with approximately 1-2 m per year (Van der Meijden et al., 2023).
Since the migration rate at the HKZ site is relatively low, for the second site a sand wave field is chosen where higher migration rates are observed.Based on Van der Meijden et al. (2023) the area west of Texel is selected (see Fig. 1A and B), where migration rates of over 10 m per year are observed.The sand waves in this area have a wave height of approximately 2 m and a wavelength between 200 and 300 m.The average depth is approximately 23 m at this transect.At this location no hydrodynamic field measurements are available.

Measurement data
At both sites Multibeam Echosounder (MBES) bathymetry data is made available by the Netherlands Hydrographic Office of the Royal Dutch Navy (NLHO).At the HKZ study site this consists of highresolution data for the years 2000 and 2010, composed by Deltares (2016Deltares ( ) from different datasets, spanning 1999Deltares ( -2001Deltares ( and 2009Deltares ( -2012 respectively.Additionally, a 2016 MBES survey, which was carried out in relation to wind farm development at the site (Deltares (2016) and Fugro (2016)), is available via RVO (2016).At the Texel site a MBES dataset is available for 2002 (only part of the area shown in Fig. 1D) and 2008 from the NLHO database, which is supplemented with 2010 MBES data from Rijkswaterstaat.The NLHO data is available via Deltares (2017).
Prior to the development of a wind farm in the HKZ area, two ADCP measurements buoys, HKZA and HKZB, were deployed from June until June 2018 (see Fugro and Deltares, 2018).By deploying two ADCP buoys, approximately 2 km apart, the measurement system was made more robust in case of technical difficulties with the equipment.For this study data from the HKZB ADCP buoy, which has the highest data availability for current measurements, is used.The location of the HKZB ADCP buoy is shown by the cross in Fig. 1C.The ADCP buoys measured, among others, the current profile over depth, at intervals of 2 m between 4 and 20 m below the surface, and the local water level, at a constant interval of 10 min.The measurement data is available at RVO (2018).
Unfortunately, the water level data is scarce due to problems with the measurement equipment and only covers 26% of the measurement period (~6 months).Moreover, both ADCP buoys showed an offset of the measured water level, which varied between the ADCP buoys and the measured periods.Comparing the results to a large-scale model the offset of the HKZB buoy was found to be in the order of a few decimeters, with a larger offset in the second measurement period after interception (see Fugro and Deltares (2018)).For this reason, only the first period of water level measurements (approximately 4.5 months) from the HKZB ADCP buoy, is used in further analysis.At the Texel site no hydrodynamic measurements are available.

Delft3D flexible mesh
In this study the newly developed Delft3D Flexible Mesh (FM) model is used to simulate hydrodynamics, sediment transport and morphology  (Deltares, 2023).The Delft3D FM model is the successor of the Delft3D-4 (or Delft3D-FLOW) model (Lesser et al., 2004), which has extensively been used for sand wave modelling (see Borsje et al. (2013);van Gerwen et al. (2018); Wang et al. (2019); Damveld et al., 2020aDamveld et al., , 2020b;;Leenders et al. (2021); Krabbendam et al. (2021) and Liang et al. (2022)).One of the main differences between the Delft3D-4 and Delft3D FM model is the ability to use unstructured grids (flexible meshes) in the latter (Kernkamp et al., 2011).Another key advantage of the Delft3D FM model for this study is the ability to run morphological simulations in parallel (i.e. on multiple nodes).This allows for optimal use of the computational power available and a significant decrease in runtime.Furthermore, the Delft3D FM model has a time-varying timestep, which is set automatically based on a user defined maximum flow Courant number.

Model description
The basic model equations as used in this study are listed below.For a more intricate description of the model reference is made to the user manual (Deltares, 2023).

Hydrodynamics.
The hydrodynamics within the Delft3D FM model are solved using the Shallow Water Equations (SWE), consisting of a continuity equation and momentum equations in three directions.Here, the horizontal scales are assumed to be much larger than the vertical scales, reducing the vertical momentum balance to the hydrostatic pressure.The time-integration is done implicitly, except for part of the advection term, which is solved explicitly, resulting in a Courant limitation for the timestep.Since in this study a 2DV model is applied and sigma coordinates are used for the vertical grid, the SWE reduce to: In these equations u represents the horizontal velocity in the along transect direction (x).ω is the vertical velocity in the σ-direction (with respect to the moving sigma plane).ζ represents the water level with respect to the reference level and H represents the total water depth (water depth (d) + water level (ζ)).ρ w is the density of water, P x is the hydrostatic pressure gradient in x-direction and F x represents the horizontal Reynold's stress.ν V is the vertical eddy viscosity.Since the k− ε turbulence closure model is chosen for this study, after Borsje et al. (2013), the vertical eddy viscosity is solved using the transport of turbulent kinetic energy (k) and energy dissipation (ε).The vertical eddy viscosity is then computed as follows: Where c μ is a calibration coefficient, set as 0.09 (Deltares, 2023).
The above equations are solved using boundary conditions at the free surface and bed level.Assuming an impermeable bed and free surface, the kinematic vertical boundary conditions are found: ω| σ=− 1 = 0 and ω| o=0 = 0 The dynamic boundary conditions, representing the stresses at bed level and at the free surface are shown in the equations below.Where τ b and τ s represent the shear stress at the bed and free surface level respectively.
Due to the exclusion of wind forcing in this study, the shear stress at the water surface is zero.The shear stress at the bed can be computed using the following equation: Where u b represents the near bed velocity (first cell above the bed) and C is the user defined Chézy coefficient.

Sediment transport and bed update.
To compute sediment transport the Van Rijn 2004 equation is chosen in this study.A more detailed description of this sediment transport equation and its implementation in Delft3D FM can be found in Van Rijn et al. (2004).Since suspended sediment transport is excluded from the model runs only the equations for bed load transport are included here.
The bed load transport (S b ) is computed using the following formula: Where ρ s and d 50 are the density and the mean diameter of the considered sediment, respectively.τ b,cr is the critical bed shear stress, according to the Shields curve.D * is the dimensionless particle size and τ ′ b,c is the grain related shear stress due to currents.For the calculation of these parameters, see the extended equation scheme in the Supplementary Materials.
The adjusted bed load transport (S ′ b ), which takes the effect of bed slope into account is computed using the bed slope parameter (α bs ): Where β represents the slope of the bed in the flow direction and φ is the angle of internal friction of the sediment, which is assumed to be 30 • .Since a 2DV set-up is used there are no longitudinal bed slopes which effect the bed load transport.
After the bed load transport is computed, the bed level change is calculated using the Exner equation (continuity of sediment), which in 2DV reduces to: Where z b is the bed level and ε p is the bed porosity, which is assumed to be 0.4.

Model set-upcomputational mesh and bathymetry
In this study a 2DV model setup is applied, adapted from earlier studies by a.o.Borsje et al. (2013) and Krabbendam et al. (2021).Within the original 2DV model set-up used in Borsje et al. (2013), a domain length of 50 km is used.However, significant variations in hydrodynamics can occur over these length scales in reality, which become apparent when moving from idealized to in-situ modelling cases.Therefore, the length of the domain is reduced to make the set-up more suitable for in-situ sand wave modelling.In this way the hydrodynamics at the model boundary better match what is happening in the sand wave area.The new domain has a length of 16 km.The direction of the domain is chosen such that it matches the main tidal axis, which is a good approximation for the direction of sand wave migration (Deltares, 2016).
The horizontal and vertical grid of the Delft3D FM models are based on Van Gerwen et al. (2018).Since steep slopes are observed in the sand wave bathymetry (see Fig. 1), a horizontal grid size of 2 m is used in the P.H.P. Overes et al. sand wave area, following Van Gerwen et al. (2018).Towards the boundaries the grid size increases to almost 400 m, using a growth factor of 1.1, as advised in the user manual (Deltares, 2023).The vertical grid consists of 40 sigma layers, with increasing size from 0.05 % of the water column at the bed to 14 % at the surface.The horizontal and vertical grid are shown in Fig. 2.
In the middle of the domain, over a length of 7 km sand waves are included in the initial bathymetry of the model.To compose the initial bed level, the measured bathymetry is split into a dynamic and a static part, see Deltares (2016).The dynamic part of the bathymetry includes both sand waves and megaripples, of which the latter are dampened out through a moving average filter.The sand waves are gradually dampened in the outermost km on both sides of the 7 km long sand wave domain, after Borsje et al. (2013).This new sand wave bathymetry is superimposed upon the static bathymetry from the Dutch Continental Shelf Model (DCSM), see Fig. 3.In this way the bed level matches that of the DCSM model at the edges of the model domain, while the effect of sand waves on the hydrodynamics gradually builds up around the area of interest.

Model set-uphydrodynamics
At the far sides of the model domain two open boundaries are present.In the original set-up by Borsje et al. (2013) and the revised version by Krabbendam et al. (2021) a Riemann type forcing is applied at these boundaries.This type of forcing boundary is well suited for tidal conditions since it links the velocity and water levels through relations found in propagating tidal waves.Moreover, its weakly reflective properties dampen shorter waves, which leads to a shortened hydrodynamic spin-up.However, this connection between the water level and velocities makes it inadequate for replicating the non-tidal currents included in this study.In this study a combination of a velocity boundary at the side of the dominant tidal direction and a water level boundary on the opposite side is applied.At both sites in this study the flood current was dominant, leading to a velocity boundary in the SW and water level at the NE edge.The validation of this set-up with the ADCP measurements is shown in the next section.It was observed that even without the weakly reflective properties of the Riemann boundary the tidal wave could exit the domain smoothly and was not reflected at the boundaries.
To define boundary conditions for the 2DV sand wave models, the regional 3D DCSM is used (see Deltares, 2018).This model includes tidal water motions as well as meteorological influences and density driven flows (due to temperature and salinity).The model has been validated using water level gauges and current profile measurements, at both onand offshore locations (Deltares, 2018).Since meteorological input is available from 1995 onwards, the model may be used to hindcast a specific period.Through such a hindcast the velocity and water level at the boundaries of the sand wave models are generated.

Model set-upparameter settings
The k-epsilon turbulence model is used to determine the vertical turbulent eddy viscosity, following Borsje et al. (2013).The background horizontal eddy viscosity is set at 0.1 m 2 s -1 .The Chézy bed roughness is set at 70 m 1/2 s − 1 , after Van Gerwen et al. (2018).The maximum flow Courant number, upon which the timestep is determined, is set at 0.5.Sediment transport is simulated using the Van Rijn 2004 equations (Van Rijn et al., 2004).Since the bed load is expected to be the main mode of sediment transport in the chosen areas, suspended sediment transport is not included in the model simulations.For simplicity we modelled the sediment to be of a single fraction, which is log-uniformly distributed.
We chose a median grainsize (d n,50 ) of 350 μm at both study sites, based on Damen et al. (2018).A value for the bed slope parameter (α bs ) of 3 is applied following Van Gerwen et al. (2018).No morphological scaling is used (i.e. a morphological scale factor of 1).The morphological parameters used are not calibrated or validated.An overview of the used model parameters is given in Table 1.

Model cases
To be able to isolate the effects time-varying, non-tidal currents, the complexity of the forcing is increased in a stepwise manner.The tidal components which are included are: M2 (main tidal component), S2 (causing spring neap variation) and M4 (identified as causing sand wave migration (Besio et al., 2004)).Additionally, a time-averaged residual current is included, which also induces sand wave migration (Németh et al., 2002;Besio et al., 2003).In case only tidal forcing is included, filtering is applied to the hydrodynamic timeseries from the DCSM, to define the tidal current velocities and water levels at the boundaries for specific tidal components.The different cases are summarized in Table 2. Case I resembles the set-up applied in most previous sand wave model studies, where the M2 and M4 tidal components are included, combined with a residual current.The constant residual current is taken as the time-averaged current from the DCSM model, which includes tidal and density driven currents as well as meteorological influences.In Case II this is supplemented with the S2 tide, to test a more complete tidal forcing.The Case III model includes the full signal.At the HKZ site the model is run for the 2 years that ADCP measurements are available (June 2016-June 2018) to validate the results.For the Texel site no hydrodynamic measurements are available.There the model is run for the 2 years following the most recent bathymetry survey (February 2010-February 2012).The first two days of the model period are used as hydrodynamic spin-up.In this period the morphological changes are set to zero.The frequency of occurrence of the residual current strength over the model periods, only applied in the Case III model, is shown in Fig. 4. At the Texel site the residual current is generally stronger and shows more variation over time.However, relative to the time-averaged residual current (μ) the variation (σ) at the HKZ site is higher (i.e. a higher coefficient of variation).

Results
Using the model set-up as described above the hydrodynamics, sediment transport and resulting morphological change are simulated and analyzed within the sand wave transects.First the modelled hydrodynamics at the HKZ site are validated using field measurements.This indicates the ability of the model to reproduce time-varying, nontidal velocities and water levels.The validated model is then used to simulate sediment transport and morphological changes and compared with tidal models to indicate the relative importance of these timevarying, non-tidal hydrodynamics for sand wave dynamics.

Hydrodynamic validation
To evaluate the performance of the Delft3D FM model on reproducing the local water levels and currents, the HKZ Case III (full forcing) model is validated using ADCP measurements from the HKZB buoy (see Fugro and Deltares, 2018).The comparison between the measured and modelled velocity is made at a depth of 12 m below the surface, approximately halfway the water column.The velocity in the sand wave model is taken from the sigma layer closest to this level.The velocity parallel to the transect is extracted for the comparison.The transect is aligned with the main tidal direction.The water levels can directly be compared between the measurements and the model results.
A good match is found between the modelled and measured water levels, as shown in  2021) is included in the Supplementary Material.Compared to previous studies the accuracy of the modelled hydrodynamics has improved significantly.In the new set up the water levels and current velocities are decoupled, unlike with the Riemann boundaries, allowing for natural phase differences between tidal water levels and velocities.Moreover, these alterations allow for inclusion of non-tidal velocity and water level fluctuations.Since the hydrodynamics are well represented in the model, the next step is to include sediment dynamics.

Morphological behaviour
To test the influence of time-varying non-tidal hydrodynamics on sediment transport and morphological behavior of sand waves, transport of sediment is added to the models.First the HKZ site is discussed, where the migration rates are low, followed by the Texel site, where migrations rates are higher.The three cases which are considered are summarized in Table 2.

HKZ sitelow migration rates
To investigate the relevance of time-varying, non-tidal currents the 3 cases as presented in Table 1 are simulated for two years.To compare the spatial bathymetry changes over time between the different types of forcing, the sedimentation and erosion patterns are shown in Fig. 6.The  figures show a representative period of 50 days.Within this 50-day period the morphological changes are in the order of centimeters.To allow for intercomparison, the colour scale is kept constant between the different models.
First, the widely used sand wave model set-up, which includes the M2 and the M4 tide and a constant residual current, is tested in Case I.The results of this model show a constant sedimentation erosion pattern for each tidal cycle throughout the two-year run, see Fig. 6.Each tidal cycle, sediment is deposited on (the lower part of) the steep slope of the sand wave, which results in a constant migration of the sand wave.When we add the S2 tidal component (Case II), the sedimentation and erosion rates rise and fall with the spring and neap tide respectively (see Fig. 6).In this run periods of larger dynamics, during spring tide, are interrupted by periods where the bed level changes are close to zero.
The sedimentation-erosion pattern from the Case III model, with the full hydrodynamic signal showed more variation compared to Case I and Case II (see Fig. 6).In short periods of time, even a few tidal cycles, relatively large bed level changes took place.In between these events, the bed level was more stable.In this sedimentation-erosion pattern the spring-neap tidal cycle can still be distinguished, although it is clear that between consecutive spring tides large differences occurred due to nontidal currents.In some periods, such as around day 556 and day 582, also erosion of the steep slope of the sand waves occurred, while in the purely tidal models only net sedimentation at the steep slope was visible over the tidal cycle.This indicates that in these tidal cycles the ebb tide was dominant, while over the long term a flood dominance was found.
The steep slope of the sand wave is clearly the most dynamic part of the bathymetry in this area.The sedimentation and erosion at this slope can be seen as a proxy for sand wave migration.When we look at the volume change in this area (control box from sand wave crest to trough), we can easily compare the results of the different forcing types over time.This comparison is shown in Fig. 7 for the slope around 7600 m.The Case I model shows a constant sedimentation of the steep slope.In Case II the sedimentation volume changes with the spring-neap tidal cycle, but no erosion is found.The Case III model shows a highly variable sedimentation and erosion pattern at the steep slope.From day 572 up to day 580 we see high sedimentation rates, indicating temporarily enforced migration.The sedimentation volumes per tidal cycle reached over 15 times the sedimentation volumes from Case I, with only tidal forcing.After this event we see erosion of the steep slope around day 582.Averaged over time the Case III model shows a sedimentation rate which is a factor 2.2 higher than the sedimentation rate from Case I.It is clear form these results that the time-varying, non-tidal currents have a considerable effect on the sedimentation and erosion of sand waves at this location.

Texel sitehigh migration rates
In the HKZ case the model results indicate that the non-tidal variable currents are a key driver for the sand wave behavior.To verify whether this also holds for more dynamic sites, two-year simulations have been conducted for the Texel site.Given the high migration rates (~10 m/yr), the amount of sediment transport and morphological change is much higher in this area.The morphological change in this area within the 50day period shown in Fig. 8 is in the order of centimeters, up to just over 10 cm locally.In the Case I and Case II models, with only tidal forcing, already significant sedimentation of the steep slope is found.In the Case I model we again see a constant net sedimentation and erosion, in terms of rate and pattern, over time (see Fig. 8).With the addition of the S2 tidal component in Case II, the effect of the spring-neap tidal cycle is visible in the magnitude of the sedimentation and erosion rates as shown in Fig. 8.During neap tide still quite some sedimentation is found on the steep slope of the sand wave.When time-varying, non-tidal currents are included, we again see a highly variable sedimentationerosion pattern Fig. 4. Frequency of occurrence of residual current strength in model period (HKZ: 2016-2018, Texel: 2010-2012), from numerical simulations with a 10min output interval.The time-averaged residual current strength, used in the Case I and Case II models, is indicated by a solid line and has a strength +0.024 m/s and +0.073 m/s for the HKZ and Texel case respectively.The respective standard deviations are 0.076 m/s and 0.134 m/s.Values are extracted from regional model.Positive residual current indicates flood direction.sedimentation rates of the Case II model oscillate around this same value, due to the spring-neap tidal cycle.In the Case II model only net sedimentation of the steep slope was found.In the Case III model, the mentioned periods of amplified sedimentation are clearly visible in the results.In these periods the sedimentation rate of the Case III model is 3-4 times as high as what was found in Case I. On the other hand, around day 273 and 296 significant erosion of the steep slope was found.Averaged over time the Case III model shows a sedimentation rate which is approximately one third higher than the Case I model.

Discussion
This study shows, for the first time, the importance of time-varying, non-tidal currents for sand wave dynamics in the North Sea.The results of the models including these currents show significant time variation in bed level changes and thus shape and migration rate of sand waves.Periods of high bed mobility are interrupted by calm periods, where bed level changes are close to zero.The time-varying, non-tidal currents even cause net-erosion at locations where due to tidal forcing only netsedimentation occurs over the tidal cycle.
Through various model alterations and the application of the more efficient Delft3D FM model long timeseries including time-varying, nontidal currents could be simulated accurately.Relative to previous model studies the addition of non-tidal forcing has improved the representation of hydrodynamics in the model significantly.In the Supplementary Materials hydrodynamic comparisons using the tidal set-up from Borsje et al. (2013) and Krabbendam et al. (2021) are included.In the studies by Borsje et al. (2013) and Krabbendam et al. (2021) Riemann boundary conditions are applied, which are especially made for representing propagating tidal waves, but inadequate for non-tidal hydrodynamics.Since the time-varying, non-tidal currents are non-repetitive, the widely used morphological scale factor (Ranasinghe et al., 2011) could not be applied to speed up the morphological changes.However, using the Delft3D FM model software the runtimes of the morphological models decreased significantly relative to the previously used Delft3D-4 model.Without parallelization a relative decrease in runtime of 25-50 % was observed, which can be attributed to, among others, an increase in code efficiency.Moreover, through parallel model runs the available computational power can be utilized more effectively, especially when computational nodes have multiple cores.Using 8 nodes the runtimes of the morphological models were decreased to 90-100 h for a two-year model run, approximately 1/13 of the computational time needed using Delft3D-4 (i.e. the model suite adopted in a.o.Borsje et al. (2013) and Krabbendam et al. (2021)).
The two considered locations in the Dutch North Sea showed qualitatively similar effects of the inclusion of time-varying, non-tidal hydrodynamics on sand wave dynamics.In both cases a highly variable sedimentation erosion pattern was visible.The absolute difference between the Case I and the Case III sedimentation volumes was largest for the Texel site.Here tidal asymmetry already causes significant sedimentation, which is regularly amplified by over 5 times due to timevarying, non-tidal currents.On the other hand, the relative effect of these currents is larger at the HKZ site, where an amplification of netsedimentation by a factor 10 was found on a regular basis.Moreover, when comparing the time-averaged sedimentation rate from the Case III model to that of the Case I model, this rate more than doubled at the HKZ site, whereas at the Texel site an increase of approximately one third was found.This suggests that although time-varying, non-tidal currents are important in both areas, especially in areas with little tidal asymmetry, and thus relatively low sand wave migration rates, these currents can have a relatively large influence on sand wave dynamics.From the residual current strength over time (see Fig. 4) it was observed that the absolute size of the standard deviation was larger at the Texel site, while the size of the standard deviation relative to the time-average (the "coefficient of variation"), is larger at the HKZ size.These values could thus give an indication of the absolute and relative impact of timevarying, non-tidal currents.Since the sites were carefully selected to be representative of the range of both hydrodynamic conditions and sand wave characteristics and dynamics in active sand wave fields in the Dutch North Sea, these current variations are likely to be important throughout the Dutch North Sea.By excluding this part of the forcing, sand wave migration rates could be underestimated by a factor two.
The Van Rijn et al., 2004 bed load transport formula used in this study is based upon the principle of a threshold current velocity for sediment transport (Van Rijn et al., 2004).The results are thus highly sensitive to the used threshold, which is, in this case, based upon the Shields curve.When the maximum tidal currents are close to this threshold, the non-tidal currents will play a significant role in sand wave dynamics.In both studied sites the threshold for sediment transport was reached more often when time-varying, non-tidal currents were included.At the HKZ site in the model including these current (Case III) the threshold was reached approximately 55 % of the time versus 45-50 % for the tidal models.At the Texel site, the ability of the ebb tidal current to transport sediment was limited to the spring tides due to this threshold.This resulted in an exceedance of the threshold around 40 % of the time in the Case III model versus 25-35 % of the time in the tidal models.A good approximation of the real threshold (or transition area) is thus vital to estimate the importance of non-tidal currents.This also means that the sediment particle size used in the model significantly influences the calculated bed load transport.Since geotechnical data is  to be accurately represented in the used numerical models.
The influence of surface waves is not included in the models used in this study.However, it is expected that the inclusion of surface waves will only increase the contrast between calm and dynamic periods.This hypothesis is supported by two conclusions from a study into the effect of surface waves on sand wave dynamics by Campmans et al. (2018).Firstly, the effect of surface waves was found to mainly increase the present migration rate.This means that in absence of sand wave migration, i.e. in the calm periods, the effect of surface waves will be small.Secondly, Campmans et al. (2018) found that extreme waves contributed to sand wave migration in a disproportionate manner.Even if these conditions have a small probability of occurrence, they have a larger absolute contribution to sand wave migration than the commonly occurring smaller waves.Since periods of extreme wave action are likely combined with strong wind-driven currents, these effects will likely strengthen each other.Together this leads to even more polarized and stochastic sand wave dynamics.
Since a 2DV model set-up was applied, the along crest component of the current is not included.The main axis of the main tidal currents is well aligned with the model transect, but other currents may be present in all directions.The along crest component of the current has a typical magnitude of approximately 10 % of the current perpendicular to the crest.However, at times the along crest component is amplified and can reach magnitudes close to the tidal currents perpendicular to the crest.Although currents along the crest of the sand wave do not directly lead to sand wave migration, they do count towards reaching the threshold for sediment transport.This means that including this direction in future models will likely increase the effect of non-tidal currents.Moreover, in many areas, significant variations in the sand wave bathymetry are observed in the along crest direction.This suggests that more 3D processes might play a role.The step towards 3D sand wave models may thus further improve model results.
By including more realistic hydrodynamics in sand wave models a step is made towards engineering applications.Since the hydrodynamics in the model showed a good agreement with field measurements, future research can now focus on morphodynamics.As stated in the method section, the morphological parameters used in these model runs are not calibrated and are taken from literature.Moreover, the addition of suspended sediment transport will lead to slightly more diffused sediment transport and thereby influence the final shape and height of the sand waves (Borsje et al., 2014).These morphological models can then be used for predictions of sand wave migration as well as trench infill rates in sand wave areas.If only tidal forcing is included, sediment transport rates are easily underestimated.A constant residual current can be a poor substitute for the real, time-varying, non-tidal currents since sediment transport is not linearly dependent on currents speed.Moreover, the timing of residual currents with respect to the tide is of importance.However, tidal models do allow for straight forward upscaling through for example the morphological scale factors.A combination of both types could be a solution.The long-term tidal influence can then be combined with for example the impact of a design storm on sand wave dynamics.This opens the door for probabilistic assessments of possible future seabed levels and trench infill rates in sand wave areas.In this way sand wave models can contribute to reducing the uncertainty in and improving the design of offshore engineering solutions.

Conclusions
Based on repeated sand wave measurements in the field, we observed changes in sand wave shape during their evolution.These temporal shape variations may be explained by time-varying, non-tidal currents, which were previously not included in studying sand wave dynamics.The Delft3D FM model offers new possibilities to efficiently and accurately reproduce long timeseries of offshore tidal and non-tidal water levels and currents and their effect on sediment dynamics.Consequently, a computation time reduction of 13 times is reached, allowing for one-on-one coupling between hydro-and morphodynamics, and bringing the long simulation periods, desired in industry projects, within reach.First, the simulated hydrodynamics at an offshore wind farm site in the North Sea are validated using in-situ measurements, which shows a good match.Second, at two studied sites the time-varying, non-tidal currents were found to have a significant influence on the temporal sedimentation-erosion rates at the considered sand wave transects.Sedimentation of the steep slope was periodically amplified up to 15 times and even erosion of the steep slope was observed, where tidal models only showed a slow, but steady sedimentation.The largest amplifications were found at the site with the lowest long-term migration rate, indicating a larger relative contribution of time-varying, non-tidal currents to local sand wave dynamics there.Since the chosen sites are representative for the range of Dutch North Sea conditions, in terms of hydrodynamics and sand wave characteristics, these influences will be significant throughout the Dutch North Sea.Incorporating these types of currents in future studies will lead to better predictions of time-varying hydrodynamics and associated sand wave dynamics and trench infill rates and will allow for improved designs of offshore constructions.

Fig. 1 .
Fig. 1. (A) Locations of the two study sites on the Dutch Continental Shelf, (B) 2010 bathymetry measurement at the Texel site, including chosen transect, (C) 2016 bathymetry measurement at the HKZ site, including chosen transect (cross indicates location of HKZB buoy deployment), (D) sand wave bathymetry profile along the Texel transect, (E) zoom of Texel transect showing subsequent decay and growth of the sand waves (F) zoom of HKZ transect showing steepening of the sand waves (F) sand wave bathymetry profile along the HKZ transect.
Fig. 5.The ebb tidal amplitude is slightly overestimated in the model, by approximately 0.1 m, but overall the water levels are highly consistent between the model and the measurements.Moreover, the sand wave model is well able to reproduce the measured currents along the transect, as shown in Fig. 5.The current velocities in the sand wave model are slightly overestimated compared to the measurements.The flood tidal currents show a large spread in peak current velocities, due to non-tidal effects, which is well represented in the model.Outliers can be attributed to measurement inaccuracies.A comparison between the hydrodynamics in the sand wave model and the DCSM model is included in the Supplementary Material and indicates the quality of the nesting.Also a comparison with the previous model set-up from Borsje et al. (2013) and Krabbendam et al. (

Fig. 5 .
Fig. 5. Comparison of water levels (left) and (transect parallel) current velocities (right) between the HKZB buoy measurements and the HKZ Case III 2DV sand wave model.A linear fit through the data is included as well as the values for the Mean Average Error (MAE), Root Mean Squared Error (RMSE) and the Pearson correlation coefficient (r).The colors indicate the point density, where yellow colors show the highest density.

Fig. 6 .
Fig. 6.Time-stack plot of computed sedimentation and erosion per tidal cycle for the Case I (upper), Case II (middle) and Case III (lower) sand wave model at the HKZ site.Dotted lines added at crest and trough locations.Bed level profile based on computed levels at day 550 Case I.

Fig. 7 .
Fig. 7. Tide-averaged sedimentation and erosion volumes for all cases at the steep slope of the sand wave (between crest and trough) around 7 600 m in the HKZ transect.Red area indicates control box (from crest to trough).Dotted line shows time-averaged sedimentation of the Case III model over the full model period.

Fig. 8 .
Fig. 8. Time-stack plot of computed sedimentation and erosion per tidal cycle for the Case I (upper), Case II (middle) and Case III (lower) sand wave model at the Texel site.Dotted lines added at crest and trough locations.Bed level profile based on computed levels at day 550 Case I.

Fig. 9 .
Fig. 9. Tide-averaged sedimentation and erosion volumes at the steep slope of the sand wave (between crest and trough) around 9 900 m in the Texel transect.Red area indicates control box (from crest to trough).Dotted line shows time-averaged sedimentation of the Case III model over the full model period.

Table 1
Overview of default parameter settings.

Table 2
Considered model cases for both study locations.