Modelling the Past and Future Evolution of Tidal Sand Waves

This study focuses on the hindcasting and forecasting of observed offshore tidal sand waves by using a state-of-the-art numerical morphodynamic model. The sand waves, having heights of several meters, evolve on timescales of years. Following earlier work, the model has a 2DV configuration (one horizontal and one vertical direction). First, the skill of the model is assessed by performing hindcasts at four transects in the North Sea where sand wave data are available of multiple surveys that are at least 10 years apart. The first transect is used for calibration and this calibrated model is applied to the other three transects. It is found that the calibrated model performs well: the Brier Skill Score is ‘excellent’ at the first two transects and ‘good’ at the last two. The root mean square error of calculated bed levels is smaller than the uncertainty in the measurements, except at the last transect, where the M2 is more elliptical than at the other three transects. The calibrated model is subsequently used to make forecasts of the sand waves along the two transects with the best skill scores.


Introduction
On many continental shelves with tidal currents >0.5 m/s and a sandy bed, offshore tidal sand waves are found. Offshore tidal sand waves are rhythmic bed forms with typical crest-to-crest distances (wavelengths) of 100-1000 m, heights of 1-10 m and migration speeds of 1-10 m/year. They are found on the outer shelves of the North Sea [1], the Gulf of Cadiz [2], the Irish Sea [3] and the coastal waters of Japan [4] and Canada [5].
Due to their dynamic nature, the characteristics of sand waves (height, migration, spatial pattern) may interfere with offshore human activities. For example, migrating sand waves affect the water depth in navigation channels and they can uncover buried cables and pipelines, thereby risking damage due to, e.g., dragged fishing gear and anchors (see [6] and references therein). At the same time, cables should not be buried too deep as this is expensive and increases the risk of cable overheating [7]. Burial depth assessments are currently based on empirical methods where historical bed level trends are extrapolated into the future (see, e.g., [8,9]). However, lack of high resolution historical bathymetric data impedes the determination of trends from these data, which may result in unnecessarily large burial depths. Therefore, there is a strong need to apply more rigorous methods, based on process knowledge, that are able to optimise the burial depth of cables and pipelines.
Tidal sand waves form spontaneously when an oscillating tidal flow interacts with a wavy bottom [10]. In the vertical plane, residual circulation cells form due to a net (i.e., tidally averaged) production of vorticity. As a result, flow is accelerated towards the crest of bottom perturbations, causing a net convergence of sand at the crests. This is balanced by the slope of the perturbation, as sand moves more easily down the slope than up (gravitational effect). The balance between these forces results in a scale selection, with one wavelength growing faster than others (often referred to as the fastest growing mode). If the tidal signal is asymmetrical due to the presence of residual currents or overtides, the convergence of sand is out of phase with the bed perturbations, causing the sand waves to migrate [11]. The model by Hulscher [10] was extended in terms of the solution method and several physical processes; for an overview see, e.g., Besio et al. [6] and Leenders et al. [12] and references therein. The formation processes were identified with the use of linear stability analysis [13]. This approach identifies the fastest-growing mode and gives an indication of the migration rate, but is not able to predict the final shape of sand waves as its validity is restricted to perturbations with small amplitudes relative to the water depth. Modelling the dynamics of mature tidal sand waves requires a non-linear approach or, alternatively, can be approximated with empirical methods.
Currently, there are several alternative predictive methods available. One of them is the model of Fredsøe and Deigaard [14] which calculates a migration rate, sand wave height and length based on the amount of sand transport. However, this model is developed for sand waves formed in unidirectional currents and predictions require estimates of future sand transport. Knaapen and Hulscher [15] developed a sand wave amplitude model based on the Landau equation which is able to predict regeneration of sand waves after dredging, but this model does not provide any information on sand wave migration. Another example is the model SEDTUBE developed by Van Rijn [16], which is a 1D model that calculates a new bed level based on the divergence of sand transport. Here, the hydrodynamics are treated as model input, so there is no direct influence of bed level changes on the tidal flow. Knaapen [17] developed a predictor of sand wave migration based on the shape, but this does not take into account shape changes. Blondeaux and Vittori [18] proposed a combination of a process-based model and an empirical approach: the preferred wavelength and its migration rate are determined with a 3D linear stability model; the expected final waveheight is calculated empirically. All these models can give a good first estimate of sand wave evolution, but because there is no coupling between hydrodynamics and morphodynamics, the full evolution of sand waves over time cannot be resolved.
Németh et al. [19], Van den Berg et al. [20], Campmans et al. [21] and Van Gerwen et al. [22] all presented full, numerical process-based models to investigate sand wave height, shape and migration, and processes that affect this evolution. However, these studies were limited in the sense that they always considered highly idealised settings (i.e., initial sinusoidal sand waves with very small amplitudes, simplified tidal flow and turbulence), which makes their applicability to realistic sand wave fields not straightforward.
These considerations motivate the aims of the present study, which are twofold. The first aim is to assess the skill of a state-of-the-art numerical morphodynamic model in reproducing historical evolution of tidal sand waves (e.g., hindcast); the second aim is to examine this model's predictions of future tidal sand wave evolution (e.g., forecast). To this end, a calibrated model based on Delft3D [23] is applied to four different areas in the North Sea. Contrary to previous studies, this model uses realistic initial bathymetry, tidal flow and bed roughness as input, and model output is compared with observations done at least a decade after the initial survey.
In the next section, the study areas and the model are described, with a focus on the analysis of bathymetric and hydrodynamic data. Section 3 presents the results of model calibration, hindcasts and forecasts, all of which are discussed in Section 4. The final section contains the conclusions. Figure 1 shows an overview of the four areas considered in this study. In all of them, data of at least two bathymetric surveys are available, covering a minimum time span of 10 years. The areas are away from large-scale (order of 5-10 km size) sand banks in order to avoid interaction between these two types of bedforms [24]. Finally, environmental conditions (mean depth, tidal current, grain size) are different in these areas. Within each of these areas, a 5 km transect oriented perpendicular to the sand wave crests, is selected for detailed investigation ( Figure 1B,C). . Panel (B1-B4) are zoom-ins of these study areas, with the red line corresponding to the transects that are investigated. Panels (C1-C4) show the bed level z b (m) over distance x along the transect (km), with the red lines corresponding to surveys performed in 1994-2001 (used as initial bathymetry) and the black lines to surveys of 2011-2016. The panels (B1-B4) and (C1-C4) were drawn using survey data which are freely available from the Hydrographic Service of the Royal Netherlands Navy [25]. Panels (C1-C4) also show definitions of crest migration c crest , wavelength λ and waveheight h. The dominant tidal component is the semi-diurnal M 2 tide with a depth-averaged velocity amplitude of 0.6-0.7 m/s. This amplitude is constructed from output of ZUNO, which is a calibrated model for tides in the North Sea (for further details see, e.g., [26][27][28]). The tidal flow is characterised by a spring-neap cycle caused by the interaction of M 2 with S 2 (depth-averaged velocity amplitude of 0.2 m/s), resulting in a tidal range varying between 0.9 m (transects 1 and 2), 0.6 m (transect 3) and 1.5 m (transect 4) during neap tide and 1.5 m (transects 1 and 2), 1.0 m (transect 3) and 2.2 m (transect 4) during spring tide. Diurnal components are also present, such as K 1 and O 1 , both with depth-averaged velocity amplitudes of about 0.01 m/s. Net sand transport is mainly determined by the interaction between M 2 and its overtides M 4 and M 6 [29], the amplitudes of which are approximately 10% and 5% of the M 2 depth-averaged velocity amplitude, respectively.

Study Sites
Depths along transects 1, 2 and 3 vary between 20 and 30 m and at transect 4 between 25 and 35 m ( Figure 1B1-B4). At transect 1, the dominant sediment type found at the bed is classified as fine to medium sand, with grain sizes varying between 200 and 300 µm ( [30,31]). The sand found at transects 2 and 3 is slightly coarser, with grain sizes in the range of 275-325 µm. The coarsest material is found along transect 4 with 350-450 µm [32].
Sand waves observed at transect 1 have wavelengths (λ, crest-to-crest distances) between 150 and 600 m, heights h varying between 0.5 and 3.0 m and migration speeds c of 1.5-5.0 m/yr in a northeastern direction ( Figure 1C1

Model Description
The present study uses the numerical model Delft3D [23] in a configuration based on Borsje et al. [34], Borsje et al. [35] and Van Gerwen et al. [22]; for a detailed model description see Section S1 of the Supplementary Materials (SM). The domain has one horizontal direction (along the transect) and a vertical direction; this configuration is called 2DV. The 2DV shallow-water equations are solved in (x, σ)-coordinates in the module FLOW, with x the coordinate along a transect of length L. The area of interest has a length and is located in the middle of the computational domain (x = 0). The vertical coordinate σ is defined as where ζ is the free surface elevation with respect to reference level z = 0 and H is the water depth below this reference level. The open boundaries are located at x = −L/2 and x = L/2 ( Figure 2a) and at these boundaries, the tidal flow is forced using Riemann invariants R, which are linear combinations of depth-averaged velocities (U, normal to the open boundary) and water levels: Here, R + is imposed at the inflow boundary (x = −L/2) and R − at the outflow boundary (x = L/2). The classification of inflow or outflow boundary depends on whether the tidal wave propagates in the positive or negative x-direction [36]. Riemann boundaries are weakly reflective ( [37,38]) and tidal waves can cross the open boundaries without reflection. Turbulence is computed using the k-model and roughness at the bed is imposed as a geometrical, Nikuradse roughness length k s . This roughness length is a function of the dimensions of the local ripples and megaripples and is calculated according to Van Rijn [39]: with l r , d r , l mr and d mr the height and length of local ripples and megaripples, respectively, and d 90 the 90th percentile grain size. Form factors γ r and γ mr depend on the dimensions of ripples and megaripples. As is explained in detail in Supplementary Materials Section S1, the output from the module FLOW is used in the module SED to solve the concentration equation using equilibrium profiles for suspended sediment at the open boundaries and to calculate bed load and suspended load transport using the equations of Van Rijn [39]. The bed load transport is corrected for longitudinal slope effects with the formulation of Bagnold [40]. This formulation contains a user-defined tuning parameter α bs which is used for calibration. The module MOR calculates the divergence of the sediment transport and updates the bed level via the sediment continuity equation. The calculated bed level changes are sped up using a morphological acceleration factor (MORFAC), which is applied under the assumption that morphological changes happen on much larger timescales than hydrodynamic changes. The new bed level is then used in FLOW to calculate hydrodynamics for the next timestep.

Model Configuration
Each transect has a length L = 50 km with in the middle (from x = − /2 to x = /2) an area of interest of length = 5 km (Figure 2a). The computational domain has a variable grid size ∆x that varies from 265 m at the boundaries to 5 m in the area of interest ( Figure 2b). Between the area of interest and the boundary there is an area with a maximum horizontal grid spacing of 50 m, so that sand waves migrating into the area of interest have a similar resolution as the bathymetry data. In the vertical direction, the domain is divided into 40 non-equidistant σ-layers. The percentage of the water depth in one layer increases from 0.05% at the bed to 6.5% at the surface ( Figure 2c).  Table 1, the ones for transects 2, 3 and 4 can be found in Supplementary Materials Table S1.

Loc 1
The k-model requires values for background horizontal eddy viscosity and diffusivity, which were set to 1 m 2 s −1 [22]. The imposed Nikuradse roughness length is k s = 0.0944 − 0.0949 m, calculated using Equation (3), with the dimensions of ripples d r = 0.034 m, l r = 0.2 m [41] and megaripples as reported by [8,33] and [9] (d mr = 0.20 m, l mr = 10 m). Form factors γ r and γ mr are set to 0.7 and 1, respectively. Median sand grain sizes are set to 275 µm at transect 1, 305 µm at transects 2 and 3 and 390 µm at transect 4 [42].
The hydrodynamic timestep ∆t was set to 6 s and the total hydrodynamic runtime varied between 30 and 40 days for the different transects. After sensitivity analysis, MORFAC was set to 148, so that one M 2 tidal cycle corresponds to a little over 2 months of morphological changes. Smaller values of MORFAC yielded almost the same results (differences of a few cm, Supplementary Materials Figure S2), but required longer computation times.
These model settings gave good results for both depth-averaged velocities and water levels at transect 1 when compared to ZUNO (Supplementary Materials Figures S5 and S6). Therefore, the model was calibrated on morphodynamics by adjusting slope parameter α bs and sand transport factor f . An overview of all model settings is presented in Table 2.

Design of Experiments
Transect 1 is chosen as the calibration location. The model, which has the initial bed level as measured in 1999 ran for 12.3 years and then the modelled bed level was compared to the measured bed level of 2012. Parameters α bs and sand transport factor f are varied such that the modelled bed level compares best to the observed bathymetry (runs H1-32, Table 3). The combination of f and α bs is then used to validate the simulated sand wave evolution along transects 2-4. Here, a survey of 1994-2001 is used as initial bathymetry and the skill is assessed by comparing the model results to a bathymetric survey of 2011-2016 (runs V2-4). Finally, the surveys measured in 2012 along transects 1 and 2 are used as initial bed level, to make a forecast for the year 2024 (runs F1-2). Hindcast runs transect 1 with f varying between 0.3 and 1 and α bs varying between 3 and 6 V2-4 Hindcast runs transects 2-4 with f = 0.45 and α bs = 5.5 F1-2 Forecast runs transects 1-2 with f = 0.45 and α bs = 5.5

Analysis of Model Output
The migration rate of the crests c crest in meters per year (troughs, c trough ) is defined as the mean (i.e., averaged over the area of interest) difference in postion of the crests (troughs) between the beginning and the end of the simulation, divided by the length of the simulation ( Figure 1C1). The wavelength λ is defined as the mean spacing between two subsequent crests ( Figure 1C2). The sand wave height h is defined as the mean vertical distance between crests and their adjacent troughs ( Figure 1C4). As each of these variables varies along the domain, the standard deviation σ is calculated as with N being the number of points in the area of interest, a-a certain variable (such as bed level, crest or trough position, sand wave length, sand wave height or migration rate) and a-the spatially averaged value of that variable. The model skill is assessed by calculating the root-mean-square error (RMSE) with respect to the observations: where a is the modelled variable and b-the observed variable. The model skill is also assessed by determining the Brier Skill Score (BSS) [43]: with z ini the initial bed level of 1994-2001, z obs the observed and z mod the modelled bed level for 2011-2016. The angular brackets . denote the mean over the area of interest. This BSS value can be decomposed [44], such that: where with r ∆z mod ,∆z obs = ∆z mod ∆z obs Here, α = 1 when sand is moved to the same position as observed (phase error), β = 0 when the same of amount sand is moved as observed (amplitude error), γ = 0 when the average modelled bed level is the same as observed and is a normalisation term. In atmospheric forecasting, a forecast with a value of BSS > 0.2 is considered useful [43].

Calibration
The model configuration described previously yielded good results for water levels and depth-averaged currents compared to ZUNO (Supplementary Materials Figures S5 and S6), therefore, calibration was done in terms of morphology. Figure 3 shows the RMSE (m) of the modelled with respect to the observed bed level of 2012 along transect 1 for different values of sand transport factor f and slope parameter α bs . The minimum value of RMSE (0.17 m) is obtained with f = 0.45 and α bs = 5.5. The BSS of this hindcast is 0.74, which is considered as excellent [43]. With α = 0.74 and β = 0.01, the BSS shows that the correct amount of sand is moved to the right location.    (Table 4). Only the depth of the troughs located at −2.45 km and −1.55 km is overpredicted slightly. Figure 4b shows the modelled bed level as a function of both distance x and time T in years. The position of crests and troughs is indicated with black dots and circles, respectively. The observed position in 2012 is indicated by the black crosses and squares. The modelled trough migration rate (c trough = 2.9 ± 0.7 m yr −1 ) is close to the observed trough migration (1.9 ± 0.9 m yr −1 ), resulting in the modelled trough position close to observations of 2012 (RMSE x,trough = 16.8 m-see also the red line in Figure 4a). The modelled crest migration rate is slightly lower than observed: c crest = 1.1 ± 1.0 m yr −1 versus 2.9 ± 1.3 m yr −1 , but the position of crests in 2012 is still close to observations (RMSE x,crest = 27.0 m). These migration rates are determined from the beginning and the end of simulation as explained in Section 2.5. The slope of the crest and trough positions in Figure 4b gives an indication of changes in the migration rate over time. In the first year, the slope of most crests and troughs deviates from the rest of the simulation, when the slope is fairly constant. Around 2.1-2.2 km, a discontinuity in the trough position is observed, which happens because of a shape change. Initially, this trough is wide with a minimum value on the right side of the trough. Throughout the simulation, this trough becomes narrower and the local minimum shifts to the left. The discontinuity happens when the left local minimum becomes deeper than the right one. A summary of the observations and model results along transect 1 is given in Table 4.    (Table 5). However, there is a tendency for both crests and troughs to be slightly lower than observed, resulting in RMSE h,crest = 0.44 m and RMSE h,trough = 0.42 m. The modelled trough migration rate (c trough = 1.3 ± 1.4 m yr −1 ) is very close to the observed rate (1.1 ± 1.1 m yr −1 ), yielding a small error in modelled trough positions (RMSE x,trough = 14.5 m, see also Figure 6a). The modelled crest migration rate is slightly lower than observed: c crest = 0.4 ± 0.6 m yr −1 versus 2.2 ± 1.2 m yr −1 and consequently RMSE x,crest = 28.7 m. Figure 6a shows that the sand wave located around 0.9 km dampens until in 2002 it completely merges with the sand wave located at 0.8 km. Figure 5b shows the observed and modelled bed levels for 2011 along transect 3. The RMSE is small (0.13 m), and the value of BSS is 0.24. BSS decomposition indicates that the lower BSS value is caused by a lower value of α, which is the phase error (Table 5) The modelled evolution of sand waves along transect 4 is depicted in Figures 5 and 6c. The bed level for 2016 has a RMSE of 0.72 m and a BSS of 0.26 compared with observations, this is again due to the phase error α (Table 5). Similar to transect 2, both crests and troughs tend to become lower than observed (RMSE h,crest = 1.04 m and RMSE h,trough = 0.55 m). Figure 5c also shows that the shape of the sand waves has changed: the crests have widened and the troughs sharpened. This affects the migration rate, which is lower than observed for both crests and troughs, resulting in RMSE x,crest = 38.7 m and RMSE x,trough = 14.3 m. In 2000, the sand wave located at 1.3 km is completely merged with the crest at 1.2 km. The sand wave originally located at 1.85 km has merged with the one at 1.7 km in the year 2010.      Figure 7 shows predicted bed levels for 2024 for transects 1 ((a)-(b)) and 2 ((c)-(d)) which have the best BSS scores. The black solid lines in Figure 7a,b : c : correspond to the initial bed level measured in 2012 and the red lines denote the prediction for 2024. Along transect 1, the predicted crest migration rate is c crest = 0.8 ± 1.0 m yr −1 and for the troughs this is c trough = 2.9 ± 0.8 m yr −1 . In this period, the crests remain at the same height or become slightly lower. Most of the troughs have deepened with an average value of 12 cm. Similarly to the hindcast, the migration rate deviates slightly in the beginning, but is constant in the rest of the simulation.

Forecast
Along transect 2, the crests are not really migrating (c crest = −0.2 ± 0.7 m yr −1 : ), but they are changing their shape towards a more rounded one. The troughs are migrating with a rate of c trough = 1.5 ± 1.0 m yr −1 . Both crests are becoming lower and troughs are getting deeper, with an average value of approximately 30 cm. The model predicts that the sand wave originally located at 0.9 km merges with the one at 0.8 km in the year 2018.

Discussion
In this paper, the skill of a 2DV morphodynamic sand wave model based on Borsje et al. [34], Borsje et al. [35] and Van Gerwen et al. [22] in reproducing observed sand wave behaviour was assessed. The effect of varying gridsize, timestep, number of σlayers and MORFAC in the calibrated model, is small. These results are presented in the Supplementary Materials in Section S3 ( Figures S1 and S2). There is a small effect on the modelled bed level when a different type of boundary conditions is imposed ( Figure S3), therefore choosing a different type of boundary conditions would result in slightly different calibration parameters.
The calibrated model yields bed levels at locations 1, 2 and 3 which have RMSE values which are of the same order as the vertical uncertainty in bathymetric measurements. At the crests and troughs along transect 4, the difference between modelled and observed bed level exceeds the vertical uncertainty range. This uncertainty depends on the water depth and is estimated to have a maximum value of ±0.57 m at transect 3 and ±0.64 m at transect 4 [45], which are the shallowest and deepest transects, respectively. The maximum horizontal uncertainty of the measurements is 6.1 m at transect 3 and 6.5 m at transect 4, which is very close to the grid size. In addition :: to the uncertainty of the measurements, the actual water depth might be different because mega ripples superimposed on sand waves are not taken into account. Moreover, due to the finite resolution of the bathymetric data, the actual crests and troughs can be missed.
The results presented in this study were calibrated using the slope parameter α bs and a scaling of the total sand transport with factor f . Increasing the slope parameter, increases the amount of sand transport from crest to trough, reducing the modelled sand wave height. This increased diffusion also influences the amount of sand waves present in the domain and shifts the positions of crests and troughs. Reducing the amount of sand transport with factor f is necessary to obtain migration rates comparable with observations. However, even after calibration, there was a tendency to underestimate crest heights and to overestimate trough depths, especially at transects 2 and 4. At transect 3, the modelled sand waves were growing which is not supported by observations. The model also simulates merging of small sand waves along transects 2 and 4, but observations reveal no merging. Both changes in shape and merging could potentially influence the migration rate. However, Figure 6 shows that there is little to no impact on the behaviour of the other sand waves, as the lines of crest and trough positions do not change after a merging event.
To further investigate why crests and troughs are under-and overestimated, runs were performed where bed levels were not updated. Further details are presented in the Supplementary Materials Section S6. The total sand transport consists of bed load (q b,tot ) and suspended load (q s ) transport. This bed load transport can be further divided into an advective part q b,adv and a diffusive part caused by slope effects q b,slope (Supplementary Materials Equations (S11), (S17) and (S18)). All the components of the transport were tidally averaged (tidal averages denoted by · ) and then the divergence of this residual transport ∂ q i ∂x was computed. Figure 8a-d shows in different colours the divergence ∂ q i ∂x (m/s) of each of the residual sand transport components together with bed level z b (m) in black over distance x (m) along transects 1-4. The sediment continuity equation (Equation (S10)) states that a convergence of sand corresponds to an increase of the bed level and vice versa. Exactly at most of the crests along transects 1, 2 and 4, the total bed load gradient is positive, indicating a lowering of the bed level. Just up-and downstream of these crests, the gradient is negative. This pattern is indicative of a lowering and widening of the crests and is most clearly visible at 0.55 km along transect 1 (Figure 8a), 0.05 km along transect 2 (Figure 8b), and at all crests along transect 4 (Figure 8d). What is striking, is that this pattern is not observed in the advective part of the bed load transport, but only in the slope part. The gradient of the advective bed load transport is negative above the crests and positive slightly to the right, which results in sand wave growth and migration. This is also what is observed along transect 3 (Figure 8c), even though all transports are significantly smaller here. The contribution of the suspended load transport is smaller than the bed load contribution at all of the transects, but shows clear negative peaks at transects 1-3. These peaks are slightly out of phase with the sand wave crests, which points to migration.
Above the troughs, both the total bed load and suspended load gradients are positive at transects 1 and 3, which means that both contribute to deepening. The total bed load gradient is positive here, because the advective part is larger than the slope part. It is a little more difficult to see what happens above the troughs at transects 2 and 4. To explain this, the time evolution of two crests and troughs is plotted in Supplementary  Materials Figures S7 and S8. Along transects 2 and 4, the crest heights change much more during the first 5 years of the simulation than during the rest. For the troughs, the opposite is the case, although they generally change much more gradually. Therefore, the bed levels after 5 years of morphological simulation time are also used to calculate the divergence of each of the residual transport components. The results for transects 2 and 4 are presented in Supplementary Materials Figure S9. Compared to the transports over the initial bed level, the transport over the crests has reduced a lot and the divergence above the crests and troughs is more of an equal magnitude. Above the troughs, both the divergence values of bed load and suspended load are positive and contribute to deepening.
The divergence can be translated into expected bed level changes via the sediment continuity equation (Supplementary Materials Equation (S28)). Further details of this approach can be found in Supplementary Materials Section S6. Two new metrics are introduced to identify the effect of each of the transport components: the global growth rate σ (yr −1 ) [46] and the global migration rate V (m yr −1 ) [47]; with the subscript i corresponding to the transport component. The root mean square bed level h rms is defined as h rms = h 2 tot 1/2 with h tot = z b − z b and the bar denoting a spatial mean over the area of interest (i.e., 1 /2 − /2 · dx). Table 6 shows the global growth rates σ in yr −1 and the global migration rates V in m yr −1 for each component of the sand transport. The total growth rate caused by bed load σ b is 0 along transect 1, as the advective and slope-related part balance each other. The suspended load only contributes very little to the growth, but it causes approximately 50% of the migration observed along transect 1. The rest of the migration is built up of almost equal parts of advective and slope-related bed load transport. At transect 2, the net growth due to bed load is negative because the slope effects are stronger than the advective transport. Here, all components add positively to the migration, with bed load being the largest driver. At transect 3, the advective part of the bed load transport causes the sand waves to grow. This is somewhat balanced by the slope-induced transport, but still the net effect is positive. This matches with the hindcast run, where sand wave growth was observed. The suspended load contributes only a little bit to the growth, but causes most of the migration along this transect. At transect 4, the total growth rate due to bed load is negative, as the slope-induced decay is stronger than the advective growth. The suspended load transport has a small positive contribution to the growth rate, but is not enough to balance this net decay. The global migration at transect 4 is small and almost entirely due to slope-related bed load transport. dx (m/s) over distance along transect x (km) for slope-related bed load transport (magenta), advective bed load transport (red), total bed load transport (green) and suspended load transport (blue) along transects 1(a), 2(b), 3(c) and 4(d). The black solid line corresponds to the bed level z b (m) along each transect. Note the different y-axes.
After calibration, excellent BSS scores are obtained at transects 1 and 2 and scores qualified as 'good' at transects 3 and 4 [43]. However, BSS scores should be considered carefully, which becomes clear at transect 4, where the RMSE value is still high. Previous studies showed that the BSS can increase in time, due to self-organisation of the morphological system [48]. This means that, initially, there may be some inconsistencies between hydrodynamic forcing and initial bed level. Upon starting the model, the bed level first adjusts to the imposed conditions (a "morphological spin-up") and then will start to evolve consistent with the forcing. This implies that errors may be larger at first, and decrease over time relative to the amount of morphological changes. Dam et al. [48] found that the BSS reaches a minimum value at 10-20 years after the start of the simulation and increases after this. The tidal sand wave evolution considered here happens on the same time scale as their "morphological spin-up time". To investigate whether morphological spin-up might be happening here, the evolution of BSS over time is investigated along transect 4. At transect 4, measurements were taken in 1994,2000,2003,2008,2012 and 2016. Figure 9 shows the value of the BSS of the modelled bed levels with respect to observations in all these years. The minimum BSS value is found approximately 10 years after the start of the simulation and afterwards increases with time. Table 6. Global growth rate σ (yr −1 ) and migration rate V (m yr −1 ) for each component of the sand transport. Another way to investigate morphological spin-up is to start the model from the initially measured bed level (measured in 2001 at transect 3), let the model run for a time τ so that the bed level can adjust to its forcings and take this new bed level after τ as the new initial bed level z b,2001+τ . This time τ is then the morphological spin-up time. The same method is applied to the last measured bed level (measured in 2011 at transect 3) to obtain the newly defined 2011 bed level z b,2011+τ , which is also adjusted to the tidal forcing. Then, the model is applied to this new initial state z b,2001+τ to obtain the prediction for 2011. This prediction for 2011 is then compared with z b,2011+τ and the BSS is computed. The resulting values of BSS as a function of morphological spin-up time τ is shown in Figure 10. The maximum value of BSS is obtained after a morphological spin-up time of approximately 10 years.
In this study, the calibration was performed at one transect and validated at three transects at different locations, i.e., a calibration in space, which yielded good results. However, in order to apply this model to offshore structures such as cable burials, the RMSE may need to be minimised further. This could, for example, be obtained by a location-specific calibration, which requires at least three measurements taken at different times: two to calibrate the model and one for validation; but this many surveys are often not available for one site. Another way to improve the model is to apply it to more transects, both inside and outside the North Sea. This is also impeded by the lack of data. Although the value of the BSS at transect 4 is qualified as 'good', the RMSE value is larger than the maximum uncertainty of the bathymetry measurement. One thing that sets this transect apart from the other three is the ellipticity of the M 2 tide, which is shown in Figure 11a. The ellipticity is defined as the ratio of the amplitudes of the minor and major axes of the tidal ellipse, and this value is much higher along transect 4 than at the other transects. For comparison, the ellipticity of the M 4 tide is depicted in Figure 11b. Recalibration of f and α bs did not result in a significant lowering of the RMSE at transect 4. Comparing Figure 1B4 with panels B1-B3, the bed at location 4 shows more variation perpendicular to the transect than at transects 1-3. This suggests that transect 4 is a good candidate for investigation with a 3D model. There are more degrees of freedom in a 3D model, such as the transverse slope terms and non-linear interactions between bed level perturbations in the longitudinal and transverse directions [49]. Leenders et al. [12] have studied the migration of sand waves over a sloping bathymetry using linear stability analysis, which resulted in a large variation in migration rates, depending on the location on the slope. However, finite amplitude sand waves have not been studied from a 3D point of view and this remains a topic of future research. This model has a few limitations in addition to the fact that the second horizontal dimension is not taken into account. The use of MORFAC speeds up the computations, but limits us in the sense that spring-neap cycles cannot be included in the current setting, as this would require much lower values of MORFAC [50]. Using linear stability analysis, Blondeaux and Vittori [51] found that including the spring-neap cycle significantly alters the fastest-growing wavelength, indicating that this might also have an effect on sand waves when looking at larger time scales, but this has not been studied within a finite amplitude regime yet.
Another simplification concerns the neglecting of wind-related phenomena (winddriven currents and waves). Campmans et al. [21] studied the effect of wind-driven currents and waves during storms on sand waves in the finite amplitude regime with the use of a highly idealised model. Besides assuming a constant vertical eddy viscosity, only bed load transport, no critical shear stress for erosion, monochromatic waves and no wave-current interactions, they considered the non-transient response of the sand waves to wind and waves. They concluded that wind-driven currents can alter the migration speed as obtained with forcing by tides only, depending on the angle of the wind with respect to the dominant tidal current direction. Furthermore, waves and wind waves tend to lower the sand wave heights obtained with only tides. However, when storm conditions were imposed only a fraction of the time (i.e., 1 week per year), both the modelled sand wave height and migration rates were very close to those of the tide-only case.
To explore the effect of wind-driven currents in our model, we conducted additional runs with no bed level updating for transect 1, but with a constant and uniform wind of U 10 = 8 m s −1 coming from the southwest (225 • in nautical convention) imposed at the surface. This value and direction are representative for the mean wind in that area [52]. A control experiment with a wind direction of 45 • was also performed. The resulting spatial variations residual bed load and suspended load along the transect have a slightly larger magnitude when southwestern wind is included (Supplementary Materials Figure S10a), but are slightly decreased in the control experiment (Supplementary Materials Figure S10b).
Including southwestern wind increases the depth-averaged flood velocity at x = 0 from 0.674 m s −1 to 0.677 m s −1 , while decreasing the depth-averaged ebb velocity from −0.593 m s −1 to −0.589 m s −1 . The flood duration becomes a few minutes longer and the depth-averaged residual velocity increases with 0.005 m s −1 . This causes no significant change in the global growth rates σ, but the global migration rates V are higher when southwestern wind is included: V b increases to 1.173 m yr −1 with wind and V s to 1.257 m yr −1 . In comparison, northeastern wind causes the migration to slow down (V b to 0.946 m yr −1 and V s to 0.886 m yr −1 ), resulting from a slight reduction of the flood velocity magnitude and duration, an increase of ebb velocity and a decrease of the residual velocity. Of course, both wind velocity and direction are highly variable in time. However, in this model, imposing time-varying wind forcing is not straightforward due to the morphological acceleration factor and remains a topic of future research.
Passchier and Kleinhans [53] observed no changes in tidal sand waves after storms of moderate intensity and little changes were found in the megaripples superimposed on sand waves at depths of 25-30 m. Only megaripples that were situated in shallower water (15-18 m depth) were completely obliterated during these storms, but recovered within one spring-neap cycle. This suggests that intense storms could indeed have a (temporary) effect on megaripples. This would result in a lower bed roughness, thereby affecting the tidal velocity and thus the sand transport. The effect of changes in bottom roughness on sand waves has been investigated and results (Supplementary Materials Figure S4) show that the resulting changes in the modelled bed levels are very small. Moreover, since storms only occur a few times per year, and megaripples recover relatively fast, the effect of storms on sand waves (which evolve on decadal time scales) is probably small. Confirming this by direct implementation of this feedback loop would require a smaller value of the morphological acceleration factor.
There is uncertainty regarding the imposed median grain size and the bed roughness. In this study, both are assumed to be constant along the whole transect, but this may not be the case. Several measurement campaigns found differences in median grain size between crests and troughs [54]. Damveld et al. [55] observed differences in the ripples superimposed on the sand wave troughs and crests. Besides the uncertainty in ripple and megaripple dimensions, there is an uncertainty in the roughness due to the form factors of Van Rijn [39] and due to the fact that there are different ways to impose bed roughness in Delft3D-FLOW [36]. Final contributors to the uncertainty are the fact that sorting effects and the interaction of morphology with biology are neglected. Damveld et al. [56] showed that sediment sorting effects can slightly lower the sand wave height, and the effect of benthos on sand wave heights also seems to be small [57].

Conclusions
The aim of this paper was to assess the skill of a 2DV numerical sand wave model in reproducing observed historical evolution along four transects in the North Sea. Along transects 1 and 2 the RMSE values of the modelled bed levels are in the order of a few centimetres to decimetres, which is smaller than the vertical uncertainty in bathymetry measurements and an 'excellent' Brier Skill Score is obtained at these transects. At transects 3 and 4, the obtained BSS score is 'good'. The RMSE values at transect 3 are also lower than the uncertainty, but at location 4, the RMSE values exceed the vertical uncertainty at the crests and troughs. At this transect, the tide is much more elliptical than at the others, so a further investigation of this transect with a 3D model is recommended.
Even though RMSE values are low and BSS values are high, modelled sand waves along transects 1, 2 and 4 tend to lower their crests and deepen their troughs. The first is caused by the slope effects in the bed load transport and the latter is the result of both advective bed load and suspended load transport. This is also visible in the predicted bed levels for 2024 at transects 1 and 2. Results also showed that the BSS values should be treated with care, as they may vary in time due to self-organisation of the system.