Do tidal sand waves always regenerate after dredging?

Tidal sand waves are rhythmic bedforms found on sandy continental shelves that pose a threat to offshore activities. While emphasis is placed on studying their natural morphodynamic evolution, little is known about if and how fast sand waves recover after dredging. This work presents an analysis of multibeam echosounder data collected at three former sand extraction sites on the Belgian continental shelf. At one of the sites, sand waves seemed to reappear approximately 5 years after dredging had stopped, which did not happen at the other two sites during the measurement period (5 and 9 years). The lack of recovery in those sites is likely the result of larger depths and smaller local sediment availability compared with the site where recovery occurred. Furthermore, these data reveal that in the latter site sand wave recovery was established mainly through local sediment redistribution.


Introduction
Tidal sand waves, sometimes referred to as large to very large dunes (Ashley, 1990), are rhythmic bedforms with wavelengths of hundreds of meters, heights of several meters and they migrate over the sandy bed of continental shelf seas at rates of 1-10 m per year (McCave, 1971;Damen et al., 2018).They form spontaneously as a free instability of an oscillating tidal current interacting with an erodible bed (see the review of Besio et al., 2008).Tidal sand waves are often observed together with larger and smaller bedforms, such as sandbanks with wavelengths of several kilometers and (mega-)ripples with wavelengths of decimeters to several tens of meters.This occurs, for example, on the Scotian Shelf (Li and King, 2007), in the southern North Sea (Bellec et al., 2010) and in the northern South China sea (Zhou et al., 2018).
Tidal sand waves are dredged when they become hazardous for navigation (Katoh et al., 1998;Dorst et al., 2011).Frequently, they are also removed as part of sand mining operations (Van Lancker et al., 2010) or installation of offshore structures, such as wind farms (Coates et al., 2015).Little is known about the morphodynamic response to dredging.This brings up the questions: do sand waves regenerate after a dredging event?And if so, on what time scales?And if not, why not?Answering these questions will help optimise dredging activities in navigation channels and wind farms and improve understanding of the impact of these activities on benthic communities that live among tidal sand waves (Cheng et al., 2021;Wyns et al., 2021).
Observations of sand waves after dredging are scarce, but analysis of field data by Katoh et al. (1998) in the Seto Inland Sea (Japan) suggests recovery of tidal sand waves within a period of approximately 10 years.Furthermore, field data in the Rotterdam-Eurogeul and Amsterdam-IJgeul, two navigation channels in the North Sea, point to regeneration of tidal sand waves within several years (Dorst et al., 2011).Campmans et al. (2021), using a highly idealised process-based model, found that sand waves tend to recover after dredging and proposed that the recovery timescale is a function of both the dredged volume and the dredging strategy.However, their model results have not yet been compared with observations, and the scarcity of such observations raises the question whether this is generally the case.In this study, bathymetric surveys, that have been conducted every few years using multibeam echosounders in areas previously subject to sand extraction in the Belgian part of the North Sea, are analysed to investigate whether tidal sand waves regenerated after they had been dredged.This analysis shows that this is not so evident.
Section 2 starts with a description of the study areas, followed by a description of the multibeam echosounder data and the methodology to analyse these in Section 3. Results will be presented in Section 4 and discussed in Section 5.The conclusions are given in Section 6.

Study areas
The area studied comprises the Kwinte Bank and Buiten Ratel, two of the large-scale Flemish banks on the Belgian continental shelf (Fig. 1).Both sandbanks are covered by sand waves and (mega-)ripples (small to large dunes).Here, the semi-diurnal M 2 -tide is dominant with depthaveraged velocity amplitudes of about 0.7-0.75m/s and an ellipticity of approximately 0.2 (obtained from the model ZUNO, for details see e. g.Zijl, 2013).The joint action with S 2 (depth-averaged amplitude of 0.2 m/s) causes spring-neap cycles, with maximum spring-flood currents of about 1 m/s.The spring-neap cycles are also visible in the tidal range which is about 5 m during spring and 3 m during neap tide.
Since the 1970s, sand has been extracted from multiple parts of these banks.This study considers three sites that have been closed to mining activities (Roche et al., 2017).Two areas (KBMA and KBMB) are located on the Kwinte Bank, where the largest amount of sand has been removed (Degrendele et al., 2010)  The third site (BRMC) is located on the Buiten Ratel sandbank.On BRMC, sand extraction took place mainly between 2008 and 2014, which has lowered the average depth to a level of 22-23 m LAT.The BRMC area was closed to extraction in early 2015 (Roche et al., 2017).This area will be called the 'western depression' from now on (Fig. 1 (e)).The three depressions are the focus of this study.
The sand in the central depression has a median grain size d 50 of 300-350 μm, is moderately-poorly sorted and contains shell fragments (De Backer et al., 2011).The sand in this depression is distinctively different from the rest of KBMA (Bellec et al., 2010).In the northern and western depression, the sand is coarser with median grain sizes of 350-400 μm (De Backer et al., 2011) and 350-450 μm, respectively (De Backer et al., 2017).In these areas, extraction caused near-complete depletion of the Holocene top layers (0-2 m left) (Degrendele et al., 2017).In KBMA, however, a thickness of 1-3 m remained and this area has been reopened to sand extraction since January 1st 2021 (Degrendele et al., 2021).The layer below the Holocene top layer consists of Quaternary deposits of varying origins: coastal, fluvial, lagoon and marsh-type deposits which are not all sandy in nature (Bellec et al., 2010).

Multibeam echosounder data
In KBMA, bathymetric data obtained during 31 surveys that took place in the period 2000-2019 were analysed, in KBMB 30 surveys from the period 2003-2019 and in BRMC 18 surveys from 2010 to 2019 (for further details see Supporting Information (SI), Table S1).Bathymetric data were obtained with two generations of the Kongsberg multibeam echosounders (MBES): the RV Belgica 100 kHz EM1002 in the period 2000-2008 and the RV Belgica 300 kHz EM3002 dual from 2009 to 2020 (Roche et al., 2017).Both generations of MBES are compliant with the IHO-S44 standard: the EM1002 has standard Order 1 with an uncertainty of ±0.6 m at depths of 20 m, and the EM3002 has the standard Special order, with an uncertainty of ±0.3 m at depths of 20 m.Roll, heave and heading corrections happened real time using a Seatex MRU5 motion sensor and an Anschutz Standard 20 gyrocompas.Data were corrected for the tide with a tidal reduction method using tide gauge data obtained from multiple stations along the Belgian Coast (Degrendele et al., 2010).Different GPS positioning systems were used: Sercel NR103 (uncertainty <5 m) for the EM1002 time series and a Thales Aquarius 02 (decimetric uncertainty) for the EM3002 dual dataset.
Multibeam echosounder measurements also provide an acoustic backscatter signal, which measures the intensity of the acoustic energy scattered back to the receiving sonar antenna.For a given angle of incidence, the backscatter strength varies with the nature of the seabed: hard, rough bottoms made up of coarse sediments return much more energy than smoother, soft bottoms made up of fine sediments.For this reason, over the past two decades, the backscatter has been used more and more as a proxy to characterise the seabed nature (Lurton et al., 2015).It is not straightforward to compare backscatter signals obtained with different uncalibrated MBES and therefore only backscatter from the EM3002 dual dataset are included in the analysis.These signals were corrected for level attenuation and instantaneous insonified area.The acoustic signal within the angular interval ±[30 • , 50 • ] that best reflects the nature of the sediment, is interpolated on a grid.These signals are hereafter averaged over a geographical area, so that time variations of these averages reflect a change of the seafloor.This standardised method of processing backscatter signals is described in detail in Roche et al. (2018).

Analysis
The bathymetric models derived from the multibeam echo soundings have a resolution of 1 × 1 m and highlight bedforms of different dimensions (Continental Shelf Service, Federal Public Service Economy, 2022).In order to isolate the tidal sand waves, these data were processed in a way similar to Van Dijk et al. (2008) and Cazenave et al. (2013).For both KBMA, KBMB and BRMC, a fixed grid was defined and data from every survey were interpolated on this grid, to ensure a constant domain size.Data gaps were filled using 2D linear interpolation to obtain bed levels z b in m at each grid point.A 2D tapered cosine window with lobes of 10% of the domain length in each direction was applied before taking the 2D Fourier transform to suppress the effect of the edges.Hereafter, a fourth order, high-pass Butterworth filter was applied in the wave number space, so that all bedforms with wavelengths λ smaller than 50 m are retained (λ cutoff = 50 m).These bedforms are assumed to be (mega-)ripples h mr in m.This high-pass bed level signal was subsequently subtracted from z b , so that the result contains only the tidal sand waves and sandbanks (h sbsw in m).Hereafter, the h mr topography was no longer considered.This resulting h sbsw was interpolated on three transects, drawn in the depressions (Fig. 1); additional transects are considered in the SI (Fig. S1).These transects are chosen in such a way, that they are oriented perpendicular to the sand wave crests.Along these transects, sandbanks (h sb ) were removed using a Gaussian moving-average window with a size of 400 m, which is twice the average sand wavelength.This resulted in sand wave topography h sw (m).An overview of the methodology applied to obtain h sw from the bathymetric models is provided in Fig. 2. The root-mean-square wave height h rms (m) was calculated from h sw as: where the bar indicates averaging, i.e. 1

Lx
∫ ⋅dx, along a transect with length L x and x the along-transect coordinate.The crest (V c ) and trough volumes (V t ) are as follows; where h crest = (h sw > 0) and h trough = (h sw < 0).These volumes are interpreted as bed level volumes per unit width (m 3 /m, see also the shaded areas in Fig. 3).Along each transect, unfiltered data are used to compute the average bed level z b (m -LAT) to investigate trends in background bathymetry.The unfiltered bed level variance σ 2 (m 2 ) is computed to consider temporal changes in tidal sand waves and megaripples combined.Unlike the sand waves in the central depression, those in the northern depression do not recover during the measurement period (Fig. 4).There are no clear sand waves left after the extraction stopped in 2010; and no increase in h sw is observed.The profile obtained in 2019 shows even less variation than in 2010.The small peak located at x = 1300 m in 2003 migrates with a nearly constant velocity of 6 m/yr.

Fig
Similarly, the sand wave height in the western depression decreased rapidly in the years prior to the closure of BRMC and no changes are visible in the period 2015-2019 (Fig. 5).The profiles obtained in 2014 and 2019 are merely shifted horizontally along the transect.The average wavelength along this transect is approximately 155 m and the migration rate is fairly constant at 0-2 m/yr.
Figs. 6 (a) shows the root-mean square wave height h rms along each transect as a function of time.In all three cases, h rms decreased prior to the closure of the sites (dashed, coloured lines).In the northern and western depression, h rms was constant after closure at 0.2 m and 0.7 m, respectively (red and green scatters).In the central depression, h rms increased from 0.4 m in 2008 to 0.65 m in 2019 (blue scatters).The mean water depth z b (m -LAT) along the transects in the depressions decreased up to the moment of closure to extraction, but did not change after (Fig. 6 (b)), suggesting that the depressions caused by sand extraction are not being filled (yet).This is supported by the evolution of the bed level variance σ 2 (m 2 ) over time (Fig. 6 (c)).In the northern and western depressions, σ 2 remained constant since closure, while there was a slight increase in the central depression, reflecting the sand wave growth.The sand waves just to the southeast and northwest of the depressions, show the same trends: both the average bed level and the variance did not change (SI, Figs.S2-S7).The crest and trough volumes V c and V t per m width (m 3 /m) evolved in the same way: both decreased prior to closure (Fig. 6 (d)).Only in the central depression did these volumes increase again after 2008, indicating that troughs were deepening and crests were growing at the same rates.In the northern and western depression, both volumes remained constant after closure.Fig. 7 shows that there is a persistent difference in backscatter signal between the central depression (blue shaded area) and that in the sand wave fields just to the north (green shaded area) and to the south (yellow shaded area) in KBMA.In the area to the north, the sand is coarser and to the south, it is finer than in the central depression.This difference is stable in time, as is visible in Fig. 7(b).Here, the backscatter signals averaged over the central depression (blue) and areas to the north (green) and south (yellow), respectively, are plotted as a function of time t (yr).The shaded areas correspond to ± one standard deviation.Although there is some overlap in the standard deviations, no significant trend is observed, suggesting a general stability of the nature of the surface sands at the time scale considered.

Discussion
Analysis of field data by Katoh et al. (1998) and Dorst et al. (2011) suggested that sand waves will regenerate within a few years after dredging.However, results from the Belgian continental shelf indicate that this is not always so evident.Tidal sand waves in the central depression (KBMA) seemed to reappear approximately 5 years after sand extraction had stopped, but in the northern and western depressions no sign of regrowth was detected during the observational period, which spans, respectively, 9 and 5 years after closure.However, regrowth may become visible in these areas in future surveys.The conclusions are independent of the type of filter, cut-off frequency or the moving average window applied (SI, Figs.S10-S12).The filter affects peaks in the bottom topography by a few decimeters, which results from the removal of the small bedforms.
To explain the differences in sand wave behaviour at the three sites, we made use of available model studies that describe the initial formation of sand waves using linear stability analysis (Besio et al., 2003;Van Santen et al., 2011;Blondeaux and Vittori, 2011;Campmans et al., 2021).This theory assumes that initially there are only small bedforms, whose amplitudes grow or decay exponentially in time, with rates that are independent of their initial amplitude.In the cases studied here, application of this theory is reasonable, because at all sites, sand waves that remained after dredging had heights on the order of 5% of the water depth panels (a)).These studies revealed that differences in growth rates are attributed to differences in amplitude and ellipticity of the depth-averaged tidal current, water depth, sediment characteristics and availability, and dredging strategy.We present an overview of these characteristics for each site in Table 1.It would be interesting to know the differences in sand waves characteristics prior to dredging, but unfortunately, these data do not exist for the Kwinte Bank.Tidal characteristics at the three sites are very similar according to the ZUNO model, and it is therefore not likely that they explain the different responses to dredging.However, in these areas, the ZUNO grid size is about 1000 m (Zijl, 2013).Observations obtained by Garel ( 2010) indicate that the current in the central depression is canalised in the

Table 1
Overview of environmental and sand wave characteristics for each of the depression areas.sense that the main axis of the tidal ellipse is oriented along the depression.Although these observations were taken over only a couple of days, they seem to indicate that the tidal currents might be slightly enhanced compared to measurements obtained just outside of the depression.To test what effect this could have, we performed simulations with the model of Besio et al. (2003), which is a local, linear model describing the initial growth of tidal sand waves on a horizontal flat seabed.Using the measured water depth, grain size and tidal flow characteristics for the central depression, the model predicts an efolding time scale (the time necessary to increase the amplitude with a factor e) of a couple of years.A 5% increase of the depth-averaged current amplitude while keeping the rest of the parameters constant, leads to a decrease of the e-folding time scale with 25%.This may be one of the reasons why growth rates are larger in KBMA.Hydrodynamic measurements were not obtained in the other two areas.Campmans et al. (2021) have shown that the recovery time may be a function of the dredged volume and the time between successive dredging events.In the years prior to closure, the yearly sand extraction rate standardised by surface area is highest for BRMC (Degrendele et al., 2010;Roche et al., 2017).However, the KBMA and KBMB exact extraction rates before 1999 are unknown and by this point the central and northern depression had already formed.It is therefore difficult to draw conclusions from the extraction rates.
The mean depths and grain sizes at the three sites are quite different, so they might explain the differences in response to dredging.Experiments with the model of Besio et al. (2003) show that a change in water depth has a much larger influence on the initial growth than a change in grain size.The e-folding time scale doubled and tripled if the water depth changed from the one from the central depression to the one of the northern (+4 m) or western depression (+7 m), respectively, while changes in grain size yielded changes in the e-folding time scale of approximately 10%.Besides the enhanced growth in shallower water due to stronger tidal bottom stress, the bed is also more prone to windand wave action in shallower water.This has the potential to mobilise more sand, as has been demonstrated by Campmans et al. (2017).
The local sediment availability is quite different between the three areas: the central depression is the only site where part of the sandy Holocene top layer remains.In the northern and western depression, the Pleistocene deposits have reached the surface (Degrendele et al., 2021).These deposits have not been studied in detail for these areas, but they have been for the central depression and the nearby Middelkerke Bank.Trentesaux et al. (1999) and Bellec et al. (2010) showed that these deposits may consist of very different materials.Besides, these sediments are assumed to be so compact that their critical shear stress for initiation of motion is larger than the maximum shear stress induced by the local water motion.Thus, there is probably no sufficient local source of erodible sand at the bed in the northern and the western depression.
In principle, sand could be transported into the depressions from surrounding areas and result in sand wave growth.Terseleer et al. (2016) indicated that this happened during the extraction period in BRMC.However, results from the KBMA area indicate that this does not happen, as there is no change in average bed level both in-and outside the central depression (Fig. S3), nor in the sand wave topography of the neighboring sand wave fields (Fig. S2).The width of the sand wave field to the northwest of the central depression did decrease somewhat (Fig. S8), but the backscatter signal in the central depression did not change over time (Fig. 7).This signal is expected to change if the sand from the northwest ended up in here as this sand has a significantly larger median grain size (see also Bellec et al., 2010).Therefore, we conclude that the growth of the sand waves in the central depression is probably the result of the local redistribution of sand: it is eroded in the troughs and deposited at the crests, which is fully in line with the free instability concept.
Several studies have looked at the effects of limited sediment availability on the formation of bedforms in a unidirectional current (Porcile et al., 2020;Jarvis et al., 2022) and in oscillatory currents (Blondeaux et al., 2016).They show that, initially, there are no differences in the formation.Once the height of the bedforms approaches the thickness of the sediment layer, however, their length tends to increase, their height to decrease and their shape becomes more three dimensional compared to the case with unlimited availability of sand.The monitoring of both the northern and western depression is still ongoing, so in the future it might be possible to observe regeneration in these areas as well.Then, the influence of local sediment availability may be visible from sand wave shape changes.These observations could be used to derive a relation between recovery time and environmental parameters such as water depth, characteristics of the bed and of the tide, which could be helpful for optimizing dredging and sand extraction strategies in the future.

Conclusions
Little is known about the response of tidal sand waves after dredging, but the few available studies assume that they will recover within a few years.To further explore this, multibeam data of three areas on the Flemish banks on the Belgian continental shelf were analysed.These areas were subject to sand extraction that has since been stopped.Using Fourier analysis and filtering techniques, the tidal sand waves were separated from the underlying sandbanks and superimposed (mega-) ripples, so that the sand wave topography along transects could be followed in time.In the central depression (KBMA), tidal sand waves seemed to form approximately 5 years after extraction activities had stopped.In the northern and western depression (KBMB and BRMC), no growth occurred during the observation period, respectively 9 and 5 years after closure.Possible explanations for this difference in response time are differences in mean water depth and local sand availability, both of which slow down the recovery and are significantly different from the central depression.The local availability of erodible sand may play a role, as bed level changes and backscatter timeseries indicate that the regeneration of sand waves in the central depression originates from local reworking processes: the sand is eroded in the troughs and deposited on the crests.
resulting in the formation of two depressions: the 'central depression' in KBMA and the 'northern depression' in KBMB (Fig. 1 (c)-(d)).In the central depression, the average depth is about 15-16 m with respect to Lowest Astronomical Tide (LAT), in the northern depression the average depth is 19-20 m LAT.Both depressions are approximately 5 m deeper than its surroundings.The KBMA and KBMB areas were closed to extraction in 2003 and 2010, respectively.

Fig. 1 .
Fig. 1.(a) Location of study areas, bathymetry from EMODnet 2020.(b) zoom-in of dashed box in panel (a), reference bed level from 2000 to 2002 (source: Continental Shelf Service, Federal Public Service Economy, 2022).(c)-(e) zoom-ins of KBMA, KBMB and BRMC (obtained 19 September -25 November 2019, source: Continental Shelf Service, Federal Public Service Economy, 2022) with location of the transects of interest (black dashed lines) and depressions (the areas most affected by extraction; black solid lines).Note the different colorbars.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 3 (a) shows the sand wave topography h sw as a function of distance x (m) along the transect in the central depression at March 3rd 2003 and November 25th 2019; Fig. 3 (b) shows h sw in colors as a

Fig. 3 .
Fig. 3. (a) Sand wave topography h sw (m) as a function of distance x (m) along a transect located in the central depression in KBMA obtained on March 3rd 2003 (dashed line) and November 25th 2019 (solid line).The red and blue shaded areas denote crest-and trough volumes V c and V t , respectively (see Section 3.2).(b) Color plots of sand wave topographies h sw (m) as a function of time t (yr) and distance x (m).The dashed line corresponds to the closure of KBMA to sand extraction, and the black triangles at x = 1600 m indicate when surveys took place.The color plot is interrupted due to a data gap between 2015 and 2018.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .
Fig. 5. (a) Sand wave topography h sw (m) as a function of distance x (m) along a transect located in the western depression in BRMC obtained on September 24th 2014 (dashed line) and September 19th 2019 (solid line).(b) Color plots of sand wave topography h sw (m) as a function of time t (yr) and distance x (m).The dashed line corresponds to the closure of BRMC to sand extraction, and the black triangles at x = 2300 m indicate when surveys took place.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .
Fig.6.The root-mean-square bed level h rms (m) (a), the average bed level z b (m -LAT) (b), the bed level variance σ 2 (m 2 ) (c) and crest and trough volume V c and V t (m 3 /m, asterisks and circles, respectively) (d) as a function of time t (yr) along transects in the central, northern and western depression (blue, red and green scatters, respectively).The coloured dashed lines in both panels corresponds to the closure to sand extraction.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. (a) Backscatter signals in angular range ±[30 • , 50 • ] obtained in KBMA are averaged over the central depression (blue) and the areas to the north (green) and south (yellow).(b) Averaged signals as a function of time t (yr) (colored dots) ± one standard deviation (shaded areas).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)