Next Article in Journal
Effects of Various Fuels on Combustion and Emission Characteristics of a Four-Stroke Dual-Fuel Marine Engine
Previous Article in Journal
A Practical Trajectory Tracking Scheme for a Twin-Propeller Twin-Hull Unmanned Surface Vehicle
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Modelling the Past and Future Evolution of Tidal Sand Waves

1
Institute for Marine and Atmospheric Research, Utrecht University, Princetonplein 5, 3584 CC Utrecht, The Netherlands
2
Faculty of Engineering Technology, Twente University, De Horst 2, 7522 LW Enschede, The Netherlands
3
WaterProof B.V., IJsselmeerdijk 2, 8221 RC Lelystad, The Netherlands
*
Author to whom correspondence should be addressed.
J. Mar. Sci. Eng. 2021, 9(10), 1071; https://doi.org/10.3390/jmse9101071
Submission received: 19 August 2021 / Revised: 17 September 2021 / Accepted: 27 September 2021 / Published: 30 September 2021
(This article belongs to the Section Coastal Engineering)

Abstract

:
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 M 2 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.

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

2. Materials and Methods

2.1. Study Sites

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).
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.
Depths along transects 1, 2 and 3 vary between 20 and 30 m and at transect 4 between 25 and 35 m (Figure 1(B1–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 1(C1)). At transect 2, λ varies between 100–550 m, with h = 1.0–4.5 m and c = 0–5.0 m/yr towards the northeast (Figure 1(C2)). The sand waves along transect 3 are the smallest in height with h = 0.5–2.0 m, but the largest in length with λ = 400–850 m. These sand waves migrate with speeds of 0–4.0 m/yr in a northeastern direction (Figure 1(C3)). At transect 4, the sand waves show the most variation perpendicular to the area of interest, with λ = 100–375 m, h = 1.5–5.0 m and c = 0–2.0 m/yr again towards the northeast (Figure 1(C4)). At all four locations, sand waves are superimposed by ripples and megaripples. These megaripples have heights of up to 0.5 m and wavelengths of up to 20 m, with average values of 0.2 m and 10 m, respectively ([8,9,33]).

2.2. 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
σ = z ζ H + ζ
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:
R ± = U ± g H ζ
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]:
k s = 3 d 90 k s , g r a i n + 20 l r d r ( d r γ r ) k s , r + 1.1 γ m r d m r ( 1 e 25 d m r / l m r ) k s , m r
with l r , d r , l m r and d m r the height and length of local ripples and megaripples, respectively, and d 90 the 90th percentile grain size. Form factors γ r and γ m r 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 α b s 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.

2.3. 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).
Along the transects, surveys from the Hydrographic Service of the Royal Netherlands Navy are used to construct the initial bed level [25]. These bathymetry data have a spatial resolution of 25 m and are linearly interpolated on the grids. These water depths have the vertical reference datum of Lowest Astronomical Tide (LAT) and are converted to the vertical reference datum of Normaal Amsterdams Peil (NAP) with the matrix NLLAT2018 provided by the Hydrographic Service in 2020. At transects 1 and 2, the survey used for the initial bathymetry was performed on 15 November 1999, at transect 3 on 15 November 2001 and at transect 4 on 16 May 1994. Outside of the area of interest, the bed level decreases gradually to a flat bed ( H 0 ). This flat bed corresponds to the average depth over length L, which is H 0 = 23.6 m at transect 1, H 0 = 26.3 m at transect 2, H 0 = 21.3 m at transect 3 and H 0 = 30.7 m at transect 4 (Figure 2a). The water depth in the 12 first cells after the open boundaries is kept at a constant value.
Transect 1 was surveyed again on 8 March 2012, transect 2 at 29 February 2012, transect 3 on 22 July 2011 and transect 4 on 10 October 2016 [25]. These bathymetric data are again interpolated linearly onto the grids and converted to NAP with NLLAT2018. This data is used for comparison with modelled bed levels for the years 2011, 2012 and 2016.
The Riemann astronomic components of M 2 , M 4 , M 6 and M 0 imposed at the open boundaries are constructed using model output from ZUNO. The boundary conditions for transect 1 are presented in Table 1, the ones for transects 2, 3 and 4 can be found in Supplementary Materials Table S1.
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 m r = 0.20 m, l m r = 10 m). Form factors γ r and γ m r 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 α b s and sand transport factor f. An overview of all model settings is presented in Table 2.

2.4. 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 α b s 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 α b s 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).

2.5. Analysis of Model Output

The migration rate of the crests c c r e s t in meters per year (troughs, c t r o u g h ) 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 1(C1)). The wavelength λ is defined as the mean spacing between two subsequent crests (Figure 1(C2)). The sand wave height h is defined as the mean vertical distance between crests and their adjacent troughs (Figure 1(C4)). As each of these variables varies along the domain, the standard deviation σ is calculated as
σ a = i = 1 N ( a i a ¯ ) 2 N
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:
R M S E a = i = 1 N ( a i b i ) 2 N
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]:
B S S = 1 ( z m o d z o b s ) 2 ( z i n i z o b s ) 2
with z i n i the initial bed level of 1994–2001, z o b s the observed and z m o d 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:
B S S = α β γ + ϵ 1 + ϵ
where
α = r Δ z m o d , Δ z o b s 2 , β = r Δ z m o d , Δ z o b s σ Δ z m o d σ Δ z o b s 2 ,
γ = Δ z m o d Δ z o b s σ Δ z o b s 2 , ϵ = Δ z o b s σ Δ z o b s 2
with r Δ z m o d , Δ z o b s = Δ z m o d Δ z o b s σ Δ z m o d σ Δ z o b s , Δ z m o d = z m o d z i n i , Δ z o b s = z o b s z i n i .
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].

3. Results

3.1. 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 α b s . The minimum value of RMSE (0.17 m) is obtained with f = 0.45 and α b s = 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.
Figure 4a compares the modelled bed level for 2012 (red line) with bathymetry measurements for this year along transect 1 (black solid line). The black dashed line indicates the initial bed level of 1999. Both the modelled crest heights and trough depths correspond well with observations, yielding a RMSE h , c r e s t = 0.13 m and a RMSE h , t r o u g h = 0.21 m, respectively (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 t r o u g h = 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 , t r o u g h = 16.8 m—see also the red line in Figure 4a). The modelled crest migration rate is slightly lower than observed: c c r e s t = 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 , c r e s t = 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.

3.2. Validation

Figure 5a compares the modelled bed level for 2012 (red line) with bathymetry measurements for this year along transect 2 (black solid line). The modelled bed level compares well with observations with BSS = 0.54 and RMSE = 0.35 m (Table 5). However, there is a tendency for both crests and troughs to be slightly lower than observed, resulting in RMSE h , c r e s t = 0.44 m and RMSE h , t r o u g h = 0.42 m. The modelled trough migration rate ( c t r o u g h = 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 , t r o u g h = 14.5 m, see also Figure 6a). The modelled crest migration rate is slightly lower than observed: c c r e s t = 0.4 ± 0.6 m yr 1 versus 2.2 ± 1.2 m yr 1 and consequently RMSE x , c r e s t = 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). According to the model, the sand waves are growing, i.e., the crests are higher and troughs are deeper than at the start of the simulation. This does not happen in the observations, resulting in RMSE h , c r e s t = 0.12 m and RMSE h , t r o u g h = 0.22 m. The modelled migration of both crest and troughs is underpredicted: observations give c c r e s t = 2.8 ± 1.2 m yr 1 and c t r o u g h = 2.3 ± 0.9 m yr 1 , respectively, while the model predicts c c r e s t = 1.0 ± 1.2 m yr 1 and c t r o u g h = 0.3 ± 0.5 m yr 1 . From Figure 6b it can be seen that the modelled crest and trough positions for 2011 are almost always to the left of the observed ones, resulting in RMSE x , c r e s t = 19.7 m and RMSE x , t r o u g h = 22.5 m.
The modelled evolution of sand waves along transect 4 is depicted in Figure 5 and Figure 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 , c r e s t = 1.04 m and RMSE h , t r o u g h = 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 , c r e s t = 38.7 m and RMSE x , t r o u g h = 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.

3.3. Forecast

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,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 c r e s t = 0.8 ± 1.0 m yr 1 and for the troughs this is c t r o u g h = 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.
Along transect 2, the crests are not really migrating ( c c r e s t = −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 t r o u g h = 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.

4. 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 α b s 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 , t o t ) and suspended load ( q s ) transport. This bed load transport can be further divided into an advective part q b , a d v and a diffusive part caused by slope effects q b , s l o p e (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];
σ i = 1 h r m s 2 h t o t h i t ¯
V i = 1 h t o t x 2 ¯ h t o t x h i t ¯
with the subscript i corresponding to the transport component. The root mean square bed level h r m s is defined as h r m s = h t o t 2 ¯ 1 / 2 with h t o t = z b z b ¯ and the bar denoting a spatial mean over the area of interest (i.e., 1 / 2 / 2 · d x ).
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.
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.
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 α b s did not result in a significant lowering of the RMSE at transect 4. Comparing Figure 1(B4) 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 (wind-driven 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].

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

Supplementary Materials

Author Contributions

Conceptualization, A.N., H.d.S., B.B. and L.P.; methodology, J.K., A.N., H.d.S., B.B. and L.P.; software, J.K.; validation, J.K.; formal analysis, J.K.; investigation, J.K.; resources, A.N., B.B. and L.P.; data curation, J.K.; writing—original draft preparation, J.K.; writing—review and editing, A.N., H.d.S. and B.B.; visualization, J.K.; supervision, A.N. and H.d.S.; project administration, H.d.S.; funding acquisition, J.K., A.N., H.d.S., B.B. and L.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research is part of the Industrial Doctorates programme with project number NWA.ID. 17.038 funded by Netherlands Organisation for Scientific Research (NWO).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Acknowledgments

We thank K. Koudstaal (WaterProof B.V.) for providing technical support.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. McCave, I.N. Sand waves in the North Sea off the coast of Holland. Mar. Geol. 1971, 10, 199–225. [Google Scholar] [CrossRef]
  2. Lobo, F.; Hernández-Molina, F.; Somoza, L.; Rodero, J.; Maldonado, A.; Barnolas, A. Patterns of bottom current flow deduced from dune asymmetries over the Gulf of Cadiz shelf (southwest Spain). Mar. Geol. 2000, 164, 91–117. [Google Scholar] [CrossRef]
  3. Van Landeghem, K.J.; Uehara, K.; Wheeler, A.J.; Mitchell, N.C.; Scourse, J.D. Post-glacial sediment dynamics in the Irish Sea and sediment wave morphology: Data–model comparisons. Cont. Shelf Res. 2009, 29, 1723–1736. [Google Scholar] [CrossRef]
  4. Katoh, K.; Kume, H.; Kuroki, K.; Hasegawa, J. The Development of Sand Waves and the Maintenance of Navigation Channels in the Bisanseto Sea. Coast. Eng. 1998, 3490–3502. [Google Scholar] [CrossRef]
  5. Duffy, G.P.; Hughes-Clarke, J.E. Application of spatial cross correlation to detection of migration of submarine sand dunes. J. Geophys. Res. Earth Surf. 2005, 110. [Google Scholar] [CrossRef] [Green Version]
  6. Besio, G.; Blondeaux, P.; Brocchini, M.; Hulscher, S.J.; Idier, D.; Knaapen, M.A.; Németh, A.A.; Roos, P.C.; Vittori, G. The morphodynamics of tidal sand waves: A model overview. Coast. Eng. 2008, 55, 657–670. [Google Scholar] [CrossRef] [Green Version]
  7. Det Norske Veritas. Subsea Power Cables in Shallow Water Renewable Energy Applications, Report DNV-RP-J301; Det Norske Veritas: Høvik, Norway, 2014. [Google Scholar]
  8. Paulsen, B.T.; Roetert, T.; Raaijmakers, T.; Forzoni, A.; Hoekstra, R.; Van Steijn, P. Morphodynamics of Hollandse Kust (zuid) Wind Farm Zone, Report 1230851-000-HYE-0003; Deltares: Delft, The Netherlands, 2016. [Google Scholar]
  9. Raaijmakers, T.; Roetert, T.; Bruinsma, N.; Riezebos, H.J.; Van Dijk, T.; Forzoni, A.; Vergouwen, S.; Grasmeijer, B. Morphodynamics and Scour Mitigation for Hollandse Kust (noord) Wind Farm Zone, Report 11202796-000-HYE-0002; Deltares: Delft, The Netherlands, 2019. [Google Scholar]
  10. Hulscher, S.J.M.H. Tidal-induced large-scale regular form patterns in a three-dimensional shallow water model. J. Geophys. Res. 1996, 101, 20727–20744. [Google Scholar] [CrossRef]
  11. Besio, G.; Blondeaux, P.; Brocchini, M.; Vittori, G. On the modeling of sand wave migration. J. Geophys. Res. Oceans 2004, 109. [Google Scholar] [CrossRef] [Green Version]
  12. Leenders, S.; Damveld, J.; Schouten, J.; Hoekstra, R.; Roetert, T.; Borsje, B. Numerical modelling of the migration direction of tidal sand waves over sand banks. Coast. Eng. 2021, 163, 103790. [Google Scholar] [CrossRef]
  13. Dodd, N.; Blondeaux, P.; Calvete, D.; de Swart, H.E.; Falqués, A.; Hulscher, S.J.M.H.; Rózyński, G.; Vittori, G. Understanding Coastal Morphodynamics Using Stability Methods. J. Coastal Res. 2003, 19, 849–865. [Google Scholar]
  14. Fredsøe, J.; Deigaard, R. Mechanics of Coastal Sediment Transport; World Scientific: Singapore, 1992; Volume 3, pp. 265–279. [Google Scholar] [CrossRef]
  15. Knaapen, M.A.; Hulscher, S.J. Regeneration of sand waves after dredging. Coast. Eng. 2002, 46, 277–289. [Google Scholar] [CrossRef] [Green Version]
  16. Van Rijn, L.C. Principles of Sedimentation and Erosion Engineering in Rivers, Estuaries and Coastal Seas Including Mathematical Modelling Package (Toolkit on CD-ROM); Aqua Publications: Amsterdam, The Netherlands, 2005. [Google Scholar]
  17. Knaapen, M.A.F. Sandwave migration predictor based on shape information. J. Geophys. Res. Earth Surf. 2005, 110, 1–9. [Google Scholar] [CrossRef] [Green Version]
  18. Blondeaux, P.; Vittori, G. A model to predict the migration of sand waves in shallow tidal seas. Cont. Shelf Res. 2016, 112, 31–45. [Google Scholar] [CrossRef]
  19. Németh, A.A.; Hulscher, S.J.; Van Damme, R.M. Modelling offshore sand wave evolution. Cont. Shelf Res. 2007, 27, 713–728. [Google Scholar] [CrossRef]
  20. Van den Berg, J.; Sterlini, F.; Hulscher, S.J.; Van Damme, R. Non-linear process based modelling of offshore sand waves. Cont. Shelf Res. 2012, 37, 26–35. [Google Scholar] [CrossRef]
  21. Campmans, G.H.; Roos, P.C.; de Vriend, H.J.; Hulscher, S.J. The Influence of Storms on Sand Wave Evolution: A Nonlinear Idealized Modeling Approach. J. Geophys. Res. Earth Surf. 2018, 123, 2070–2086. [Google Scholar] [CrossRef] [Green Version]
  22. Van Gerwen, W.; Borsje, B.W.; Damveld, J.H.; Hulscher, S.J. Modelling the effect of suspended load transport and tidal asymmetry on the equilibrium tidal sand wave height. Coast. Eng. 2018, 136, 56–64. [Google Scholar] [CrossRef] [Green Version]
  23. Lesser, G.R.; Roelvink, J.A.; van Kester, J.A.; Stelling, G.S. Development and validation of a three-dimensional morphological model. Coast. Eng. 2004, 51, 883–915. [Google Scholar] [CrossRef]
  24. Komarova, N.L.; Newell, A.C. Nonlinear dynamics of sand banks and sand waves. J. Fluid Mech. 2000, 415, 285–321. [Google Scholar] [CrossRef] [Green Version]
  25. Hydrographic Service of the Royal Netherlands Navy. Hydrographic Surveys. 2017. Available online: http://opendap.deltares.nl/thredds/catalog/opendap/hydrografie/surveys/catalog.html (accessed on 17 November 2018).
  26. Dastgheib, A.; Roelvink, J.; Wang, Z. Long-term process-based morphological modeling of the Marsdiep Tidal Basin. Mar. Geol. 2008, 256, 90–100. [Google Scholar] [CrossRef]
  27. Zijl, F. Development of the next generation Dutch Continental Shelf Flood Forecasting Models, Memo, 1205989-003-ZKS-0002; Deltares: Delft, The Netherlands, 2013. [Google Scholar]
  28. Zijl, F.; Verlaan, M.; Gerritsen, H. Improved water-level forecasting for the Northwest European Shelf and North Sea through direct modelling of tide, surge and non-linear interaction. Ocean Dyn. 2013, 63, 823–847. [Google Scholar] [CrossRef]
  29. Van De Kreeke, J.; Robaczewska, K. Tide-induced residual transport of coarse sediment; Application To the Ems Estuary. Neth. J. Sea Res. 1993, 31, 209–220. [Google Scholar] [CrossRef]
  30. Fugro. Geophysical Site Investigation Survey/Hollandse Kust (Zuid) Wind Farm Development Zone/Wind Farm Site I to IV/Report GH176-R1-4; Fugro: Nootdorp, The Netherlands, 2016; Volume B. [Google Scholar]
  31. Fugro. Geophysical Site Investigation Survey/Hollandse Kust (noord) Wind Farm Survey 2017/Report GH216-R3; Fugro: Nootdorp, The Netherlands, 2018; Volume 3. [Google Scholar]
  32. Damen, J.M.; van Dijk, T.A.; Hulscher, S.J. Spatially Varying Environmental Properties Controlling Observed Sand Wave Morphology. J. Geophys. Res. Earth Surf. 2018, 123, 262–280. [Google Scholar] [CrossRef]
  33. Van Dijk, T.A.; Kleinhans, M.G. Processes controlling the dynamics of compound sand waves in the North Sea, Netherlands. J. Geophys. Res. Earth Surf. 2005, 110, 1–15. [Google Scholar] [CrossRef]
  34. Borsje, B.W.; Roos, P.C.; Kranenburg, W.M.; Hulscher, S.J. Modeling tidal sand wave formation in a numerical shallow water model: The role of turbulence formulation. Cont. Shelf Res. 2013, 60, 17–27. [Google Scholar] [CrossRef]
  35. Borsje, B.W.; Kranenburg, W.M.; Roos, P.C.; Matthieu, J.; Hulscher, S.J. The role of suspended load transport in the occurence of tidal sand waves. J. Geophys. Res. Earth Surf. 2014, 119, 701–716. [Google Scholar] [CrossRef] [Green Version]
  36. Deltares. Delft3D 3D-FLOW User Manual, Version: 3.15.34158; Deltares: Delft, The Netherlands, 2014. [Google Scholar]
  37. Engquist, B.; Majda, A. Absorbing boundary conditions for the numerical simulation of waves. Math. Comput. 1977, 31, 629–651. [Google Scholar] [CrossRef]
  38. Verboom, G.K.; Slob, A. Weakly-reflective boundary conditions for two-dimensional shallow water flow problems. Adv. Water Resour. 1984, 7, 192–197. [Google Scholar] [CrossRef]
  39. Van Rijn, L.C. Principles of Sediment Transport in Rivers, Estuaries and Coastal Seas; Aqua Publications: Amsterdam, The Netherlands, 1993. [Google Scholar] [CrossRef]
  40. Bagnold, R.A. An Approach to the Sediment Transport Problem From General Physics, in US Geological Survey Professional Paper 422-I; U.S. Government Printing Office: Washington, DC, USA, 1966; pp. 1–37. [CrossRef] [Green Version]
  41. Sleath, J.F. Sea Bed Mechanics; McCormick, M.E., Ed.; Wiley: Hoboken, NJ, USA, 1984. [Google Scholar]
  42. Maljers, D.; Gunnink, J. Interpolation of Measured Grain-Size Fractions. 2007. Available online: http://www.emodnet-seabedhabitats.eu/pdf/Imares_Interpolation%20of%20measured%20grain-size%20fractions.pdf (accessed on 17 May 2020).
  43. Sutherland, J.; Peet, A.; Soulsby, R. Evaluating the performance of morphological models. Coast. Eng. 2004, 51, 917–939. [Google Scholar] [CrossRef]
  44. Murphy, A.H.; Epstein, E.S. Skill scores and correlation coefficients in model verification. Mon. Weather Rev. 1989, 117, 572–582. [Google Scholar] [CrossRef] [Green Version]
  45. International Hydrographic Organization. IHO Standards for Hydrographic Surveys—IHO Publication No. 44; International Hydrographic Bureau: Monte-Carlo, Monaco, 2008; Volume 6. [Google Scholar]
  46. Garnier, R.; Calvete, D.; Falqués, A.; Caballeria, M. Generation and nonlinear evolution of shore-oblique/transverse sand bars. J. Fluid Mech. 2006, 567, 327–360. [Google Scholar] [CrossRef] [Green Version]
  47. Vis-Star, N.C.; De Swart, H.; Calvete, D. Patch behaviour and predictability properties of modelled finite-amplitude sand ridges on the inner shelf. Nonlin. Processes Geophys. 2008, 15, 943–955. [Google Scholar] [CrossRef] [Green Version]
  48. Dam, G.; Van der Wegen, M.; Labeur, R.; Roelvink, D. Modeling centuries of estuarine morphodynamics in the Western Scheldt estuary. Geophys. Res. Lett. 2016, 43, 3839–3847. [Google Scholar] [CrossRef] [Green Version]
  49. Blondeaux, P.; Vittori, G. Three-dimensional tidal sand waves. J. Fluid Mech. 2009, 618, 1–11. [Google Scholar] [CrossRef]
  50. Ranasinghe, R.; Swinkels, C.; Luijendijk, A.; Roelvink, D.; Bosboom, J.; Stive, M.; Walstra, D. Morphodynamic upscaling with the MORFAC approach: Dependencies and sensitivities. Coast. Eng. 2011, 58, 806–811. [Google Scholar] [CrossRef]
  51. Blondeaux, P.; Vittori, G. Formation of tidal sand waves: Effects of the spring - neap cycle. J. Geophys. Res. Oceans 2010, 115. [Google Scholar] [CrossRef]
  52. Reistad, M.; Breivik, Ø.; Haakenstad, H.; Aarnes, O.J.; Furevik, B.R.; Bidlot, J.R. A high-resolution hindcast of wind and waves for the North Sea, the Norwegian Sea, and the Barents Sea. J. Geophys. Res. Oceans 2011, 116. [Google Scholar] [CrossRef] [Green Version]
  53. Passchier, S.; Kleinhans, M. Observations of sand waves, megaripples, and hummocks in the Dutch coastal area and their relation to currents and combined flow conditions. J. Geophys. Res. Earth Surf. 2005, 110. [Google Scholar] [CrossRef]
  54. Van Oyen, T.; Blondeaux, P.; Van Den Eynde, D. Sediment sorting along tidal sand waves: A comparison between field observations and theoretical predictions. Cont. Shelf Res. 2013, 63, 23–33. [Google Scholar] [CrossRef]
  55. Damveld, J.H.; van der Reijden, K.J.; Cheng, C.; Koop, L.; Haaksma, L.R.; Walsh, C.A.; Soetaert, K.; Borsje, B.W.; Govers, L.L.; Roos, P.C.; et al. Video transects reveal that tidal sand waves affect the spatial distribution of benthic organisms and sand ripples. Geophys. Res. Lett. 2018, 45, 11837–11846. [Google Scholar] [CrossRef] [Green Version]
  56. Damveld, J.; Borsje, B.; Roos, P.; Hulscher, S. Horizontal and Vertical Sediment Sorting in Tidal Sand Waves: Modeling the Finite-Amplitude Stage. J. Geophys. Res. Earth Surf. 2020, 125. [Google Scholar] [CrossRef]
  57. Damveld, J.H.; Borsje, B.W.; Roos, P.C.; Hulscher, S.J. Biogeomorphology in the marine landscape: Modelling the feedbacks between patches of the polychaete worm Lanice conchilega and tidal sand waves. Earth Surf. Processes Landforms 2020, 45, 2572–2587. [Google Scholar] [CrossRef]
Figure 1. Overview of the four study sites, with panel (A) showing their locations in the North Sea (EMODnet bathymetry 2020). Panel (B1B4) are zoom-ins of these study areas, with the red line corresponding to the transects that are investigated. Panels (C1C4) 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 (B1B4) and (C1C4) were drawn using survey data which are freely available from the Hydrographic Service of the Royal Netherlands Navy [25]. Panels (C1C4) also show definitions of crest migration c c r e s t , wavelength λ and waveheight h. Depths are relative to Normaal Amsterdams Peil (NAP) and coordinates are in World Geodetic System 84/Universal Transverse Mercator 31 N. Note the different colour bars in all plots.
Figure 1. Overview of the four study sites, with panel (A) showing their locations in the North Sea (EMODnet bathymetry 2020). Panel (B1B4) are zoom-ins of these study areas, with the red line corresponding to the transects that are investigated. Panels (C1C4) 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 (B1B4) and (C1C4) were drawn using survey data which are freely available from the Hydrographic Service of the Royal Netherlands Navy [25]. Panels (C1C4) also show definitions of crest migration c c r e s t , wavelength λ and waveheight h. Depths are relative to Normaal Amsterdams Peil (NAP) and coordinates are in World Geodetic System 84/Universal Transverse Mercator 31 N. Note the different colour bars in all plots.
Jmse 09 01071 g001
Figure 2. (a) Initial bed level z b (m) over distance x (km) along transect 1. (b) Grid spacing Δ x (m) in the horizontal direction and (c) grid spacing Δ z in the vertical direction.
Figure 2. (a) Initial bed level z b (m) over distance x (km) along transect 1. (b) Grid spacing Δ x (m) in the horizontal direction and (c) grid spacing Δ z in the vertical direction.
Jmse 09 01071 g002
Figure 3. RMSE (m) of modelled 2012 bed level along transect 1 for different values of sand transport factor f and slope parameter α b s .
Figure 3. RMSE (m) of modelled 2012 bed level along transect 1 for different values of sand transport factor f and slope parameter α b s .
Jmse 09 01071 g003
Figure 4. (a) Modelled bed level z b (m) for 2012 (red) over distance x (km) along transect 1. The black dotted line corresponds to initial bed level (1999) and the black solid line to the measured bed level of 2012. (b) Bed level z b (m) as a function of distance x (km) and time T (yr). The black dots and circles mark the position of crests and troughs, respectively, over time. The black crosses and squares correspond to the measured positions in 2012.
Figure 4. (a) Modelled bed level z b (m) for 2012 (red) over distance x (km) along transect 1. The black dotted line corresponds to initial bed level (1999) and the black solid line to the measured bed level of 2012. (b) Bed level z b (m) as a function of distance x (km) and time T (yr). The black dots and circles mark the position of crests and troughs, respectively, over time. The black crosses and squares correspond to the measured positions in 2012.
Jmse 09 01071 g004
Figure 5. Modelled bed level z b (m) for 2011, 2012 or 2016 (red) over distance x (km) along transects 2–4 (panels (ac)). Black dotted lines correspond to initial bed levels measured in 1999, 2001 or 1994 and black solid lines to bed levels observed in 2011, 2012 or 2016 used for comparison.
Figure 5. Modelled bed level z b (m) for 2011, 2012 or 2016 (red) over distance x (km) along transects 2–4 (panels (ac)). Black dotted lines correspond to initial bed levels measured in 1999, 2001 or 1994 and black solid lines to bed levels observed in 2011, 2012 or 2016 used for comparison.
Jmse 09 01071 g005
Figure 6. Bed level z b (m) as a function of distance x (km) and time T (yr) for transects 2–4 (ac). The black dots and circles mark the position of crests and troughs, respectively, over time. The black crosses and squares correspond to the measured positions in 2011, 2012 or 2016. Note the different colour bars.
Figure 6. Bed level z b (m) as a function of distance x (km) and time T (yr) for transects 2–4 (ac). The black dots and circles mark the position of crests and troughs, respectively, over time. The black crosses and squares correspond to the measured positions in 2011, 2012 or 2016. Note the different colour bars.
Jmse 09 01071 g006
Figure 7. Predicted bed levels z b (m) for 2024 (red lines) over distance x (km) along transects 1 (a) and 2 (c) and as a function of distance x (km) and time T (yr) (panels (b,d) for transects 1 and 2, respectively). The black dots and circles mark the position of crests and troughs, respectively, over time. Note the different colour bars.
Figure 7. Predicted bed levels z b (m) for 2024 (red lines) over distance x (km) along transects 1 (a) and 2 (c) and as a function of distance x (km) and time T (yr) (panels (b,d) for transects 1 and 2, respectively). The black dots and circles mark the position of crests and troughs, respectively, over time. Note the different colour bars.
Jmse 09 01071 g007
Figure 8. Horizontal gradient of tidally averaged sand transport q d x (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.
Figure 8. Horizontal gradient of tidally averaged sand transport q d x (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.
Jmse 09 01071 g008
Figure 9. Evolution of the Brier Skill Score (BSS) of modelled bed levels with respect to observations taken at time T (yr) along transect 4.
Figure 9. Evolution of the Brier Skill Score (BSS) of modelled bed levels with respect to observations taken at time T (yr) along transect 4.
Jmse 09 01071 g009
Figure 10. Evolution of the Brier Skill Score (BSS) of modelled bed levels as a function of morphological spin-up time τ (yr) along transect 3.
Figure 10. Evolution of the Brier Skill Score (BSS) of modelled bed levels as a function of morphological spin-up time τ (yr) along transect 3.
Jmse 09 01071 g010
Figure 11. Ellipticity (ratio of amplitudes of minor and major axis) of the depth-averaged M 2 (a) and M 4 tidal velocity (b) along the four transects. Ellipticity is calculated from ZUNO output.
Figure 11. Ellipticity (ratio of amplitudes of minor and major axis) of the depth-averaged M 2 (a) and M 4 tidal velocity (b) along the four transects. Ellipticity is calculated from ZUNO output.
Jmse 09 01071 g011
Table 1. Boundary conditions imposed at the open boundaries x = L / 2 and x = L / 2 .
Table 1. Boundary conditions imposed at the open boundaries x = L / 2 and x = L / 2 .
Loc 1 R + , x = L / 2 R , x = L / 2
AmplitudePhaseAmplitudePhase
M 2 1.02 ms 1 73 0.38 ms 1 73
M 4 0.14 ms 1 108 0.06 ms 1 296
M 6 0.05 ms 1 92 0.04 ms 1 89
M 0 0.02 ms 1 −0.001 ms 1
Table 2. Overview of default model settings for hydrodynamics, sand transport and numerics.
Table 2. Overview of default model settings for hydrodynamics, sand transport and numerics.
ParameterSymbolTransect 1Transect 2Transect 3Transect 4Dimension
Domain lengthL50505050km
Length of area of interest5555km
Undisturbed water depth H 0 23.626.321.330.7m
Roughness length k s 0.09440.09450.09450.0949m
Median sand grain size d 50 0.2750.3050.3050.390mm
Bed slope parameter α b s 5.55.55.55.5-
Sand transport scale factorf0.450.450.450.45-
Hydrodynamic time step Δ t 6666s
Horizontal grid spacing Δ x 5555m
No of σ -layers-40404040-
Morphological acceleration factorMORFAC148148148148-
Morphological simulation time T m 12.312.39.722.4years
Table 3. Overview of experiments.
Table 3. Overview of experiments.
RunsDescription
H1–32Hindcast runs transect 1 with f varying between
0.3 and 1 and α b s varying between 3 and 6
V2–4Hindcast runs transects 2–4 with f = 0.45 and α b s = 5.5
F1–2Forecast runs transects 1–2 with f = 0.45 and α b s = 5.5
Table 4. Observed and modelled sand wave evolution along transect 1.
Table 4. Observed and modelled sand wave evolution along transect 1.
Transect 1
ObservationsModel
c c r e s t 2.9 m yr 1 1.1 m yr 1
σ c , c r e s t 1.3 m yr 1 1.0 m yr 1
c t r o u g h 1.9 m yr 1 2.9 m yr 1
σ c , t r o u g h 0.9 m yr 1 0.7 m yr 1
RMSE h , c r e s t 0.13 m
RMSE x , c r e s t 27.0 m
RMSE h , t r o u g h 0.21 m
RMSE x , t r o u g h 16.8 m
RMSE 0.17 m
α 0.74
β 0.01
γ 0.01
ϵ 0.01
BSS 0.74
Table 5. Observed and modelled sand wave evolution along transects 2, 3 and 4.
Table 5. Observed and modelled sand wave evolution along transects 2, 3 and 4.
Transect 2Transect 3Transect 4
ObservationsModelObservationsModelObservationsModel
c c r e s t 2.2 m yr 1 0.4 m yr 1 2.8 m yr 1 1.0 m yr 1 1.3 m yr 1 −0.4 m yr 1
σ c , c r e s t 1.2 m yr 1 0.6 m yr 1 1.2 m yr 1 1.2 m yr 1 0.2 m yr 1 0.5 m yr 1
c t r o u g h 1.1 m yr 1 1.3 m yr 1 2.3 m yr 1 0.3 m yr 1 1.0 m yr 1 0.6 m yr 1
σ c , t r o u g h 1.1 m yr 1 1.4 m yr 1 0.9 m yr 1 0.5 m yr 1 0.5 m yr 1 0.5 m yr 1
RMSE h , c r e s t 0.44 m 0.12 m 1.04 m
RMSE x , c r e s t 28.7 m 19.7 m 38.7 m
RMSE h , t r o u g h 0.42 m 0.22 m 0.55 m
RMSE x , t r o u g h 14.5 m 22.5 m 14.3 m
RMSE 0.35 m 0.13 m 0.72 m
α 0.57 0.26 0.27
β 0.03 0.02 0.01
γ 0.01 0.02 0.01
ϵ 0.01 0.02 0.01
BSS 0.54 0.24 0.26
Table 6. Global growth rate σ (yr 1 ) and migration rate V (m yr 1 ) for each component of the sand transport.
Table 6. Global growth rate σ (yr 1 ) and migration rate V (m yr 1 ) for each component of the sand transport.
Transect 1Transect 2Transect 3Transect 4
σ b , a d v 0.027 yr 1 0.031 yr 1 0.0226 yr 1 0.054 yr 1
σ b , s l o p e −0.027 yr 1 −0.033 yr 1 −0.017 yr 1 −0.076 yr 1
σ b 0.000 yr 1 −0.002 yr 1 0.006 yr 1 −0.022 yr 1
σ s 0.002 yr 1 0.004 yr 1 0.003 yr 1 0.006 yr 1
V b , a d v 0.461 m yr 1 0.387 m yr 1 −0.198 m yr 1 −0.135 m yr 1
V b , s l o p e 0.598 m yr 1 0.755 m yr 1 0.447 m yr 1 0.844 m yr 1
V b 1.059 m yr 1 1.142 m yr 1 0.249 m yr 1 0.710 m yr 1
V s 1.070 m yr 1 0.779 m yr 1 0.485 m yr 1 0.067 m yr 1
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Krabbendam, J.; Nnafie, A.; de Swart, H.; Borsje, B.; Perk, L. Modelling the Past and Future Evolution of Tidal Sand Waves. J. Mar. Sci. Eng. 2021, 9, 1071. https://doi.org/10.3390/jmse9101071

AMA Style

Krabbendam J, Nnafie A, de Swart H, Borsje B, Perk L. Modelling the Past and Future Evolution of Tidal Sand Waves. Journal of Marine Science and Engineering. 2021; 9(10):1071. https://doi.org/10.3390/jmse9101071

Chicago/Turabian Style

Krabbendam, Janneke, Abdel Nnafie, Huib de Swart, Bas Borsje, and Luitze Perk. 2021. "Modelling the Past and Future Evolution of Tidal Sand Waves" Journal of Marine Science and Engineering 9, no. 10: 1071. https://doi.org/10.3390/jmse9101071

APA Style

Krabbendam, J., Nnafie, A., de Swart, H., Borsje, B., & Perk, L. (2021). Modelling the Past and Future Evolution of Tidal Sand Waves. Journal of Marine Science and Engineering, 9(10), 1071. https://doi.org/10.3390/jmse9101071

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop