Monitoring underwater nourishments using multibeam bathymetric and backscatter time series

Natural and man-induced coastal erosion endanger life and environment in coastal areas worldwide. For sedimentary barrier coasts, beach and underwater nourishments are an efficient coastal protection strategy. To optimize nourishments and to understand their impact on the marine environment, monitoring strategies are required. In this study, we investigate the potential of multibeam echosounder (MBES) data, providing both bathymetry and backscatter (BS), for monitoring the evolution of the nourished sediment and morphodynamics over time. A time series of seven MBES measurements, as well as two sets of box cores, vibrocores and seismic data were acquired of a channel-side nourishment near the Wadden Sea island Ameland (The Netherlands), between April 2017 and May 2019. In general, a high confidence of the acoustic reliability of the BS time series measurements is demonstrated. The unsupervised Bayesian classification method, supported by groundtruthing, is employed to produce a time series of sediment maps, revealing sediments ranging from sandy mud to sand with varying amounts of shell fragments. Based on the sediment maps, the nourished sediment could be distinguished from the natural sediment. Within one year, the shell-rich pre-nourishment seabed is recreated by washing out finer sediments, which are deposited towards the main tidal channel. Using the seismic data and vibrocores, the shell-rich pre-nourishment seabed could be identified in the subsurface after being buried by the nourishments, supporting the general findings. Furthermore, a rapid development of steep bedforms with increasing sediment sorting is observed in parts of the nourished areas. This study shows that high-resolution sediment maps obtained from a time series of MBES BS together with bathymetry reveal morphodynamic and sedimentary processes of nourishment evolution and can advance underwater nourishment strategies.


Introduction
Coastal erosion is a hazard to humans and infrastructure in coastal zones worldwide. Erosion is mainly driven by storm waves, crossand long-shore currents but also sea level rise leads to loss of coastal land (Woolfe, 2017). For barrier coasts like the Wadden Sea islands, beach and underwater nourishments are an efficient coastal protection strategy to further mitigate erosion and to counteract sea-level rise. In the Netherlands, yearly nourishments exceed 12 million m 3 (Pit et al., 2017). Without such human intervention it is unlikely that the natural processes can supply sufficient sediment to hold the Wadden Sea islands in a dynamic equilibrium to relative sea-level rise and further erosions (Elias et al., 2012). However, consequences of coastal engineering on the long-term response of the Wadden Sea environment is in general relatively unexplored (Elias et al., 2012). Coastal engineering, like nourishments or seawalls, can induce changes of the sediment T.C. Gaida et al. and geometrical corrections to the received signal at the MBES, a measure of the backscatter (BS) strength can be retrieved (Lurton, 2010). This is an intrinsic property of the seabed and therefore applicable to characterize seabed substrate types. As shown by several scientific studies, collecting BS data over large areas allows to extend in-situ measurements from a single location (e.g., box core samples or video footage) to a large-scale and high-resolution sediment map (Anderson et al., 2002;Eleftherakis et al., 2004;Ierodiaconou et al., 2018;Gaida et al., 2019).
In order to compare repetitive BS measurements over time for monitoring the seabed environment, external sources of variation in the BS measurement need to be investigated. External sources are either related to the acoustic system or the environment. The former includes varying sonar settings, aging antennas, varying transducer sensitivity or biofouling of the transducer. Environmental sources are in particular related to the water column properties due to the occurrence of sediment suspension, gas bubbles or varying salinity and temperature (Lurton and Lamarche, 2015). In addition, the orientation of small-scale bedforms (i.e., organized seabed roughness) relative to the navigation direction causes an azimuthal dependency of the BS strength ). An ideal design and implementation of a BS monitoring program requires full control of the acoustic system and the environment (Anderson et al., 2008). The calibration of a MBES for BS measurements, aiming for a full control of the sources related to the acoustic system, has just recently attracted more attention (Lurton and Lamarche, 2015). Currently, two approaches exist to calibrate the MBES for BS measurements: (1) an absolute calibration via a crosscorrelation with a fully calibrated reference system (e.g., singlebeam echosounder) Ladroit et al., 2018;Weber et al., 2018;Wendelboe, 2018) or (2) a relative calibration using natural reference areas with a known and temporally stable BS response (Roche et al., 2018). As long as the MBES is not calibrated, the measured and processed BS data represent a relative measurement and the correct terminology is either relative BS strength or BS level. To control environmental sources, Montereale-Gavazzi et al. (2019) used acoustic and optical BS sensors to measure the sediment suspension concentration (SSC) between 0.2 and 2.4 m above the seabed. However, the control of all parameters is often unfeasible due to time or budget limitations, environmental settings or weather conditions and alternative approaches are required as well.
To date, only a limited number of studies employed repetitive MBES BS measurements, either focusing on the temporal reliability of the BS levels (Madricardo et al., 2017;Roche et al., 2018;Montereale-Gavazzi et al., 2019) and the acoustic classification results  or using the BS measurements for seabed change detection of natural seabeds via manual interpretation (Urgeles et al., 2002) and automatic classification methods (Rattray et al., 2013;Montereale-Gavazzi et al., 2017). Roche et al. (2018) showed that the use of natural reference areas allows for comparing BS measurements over time, based on an investigation of uncalibrated BS surveys in a time span of 12 years at three shallow-water areas. Montereale-Gavazzi et al. (2019) conducted a field experiment of 15 to 47 repetitive tracks in three areas to investigate the effect of short-term hydrodynamic variations and SSC on the BS level over a tidal cycle (∼13 h) . Snellen et al. (2019) studied the repeatability of five MBES surveys acquired over 15 months in a very stable environment, and demonstrated that sediment classification based on MBES BS achieves repeatable results in this stable environment. Montereale-Gavazzi et al. (2017) used a natural reference area and a consistent data acquisition to assure the temporal stability of the measured BS data, which allowed for the quantification of temporal changes of sediment distributions, and identified the expansion of sandy areas over gravel beds. Despite the few studies, there is still a need to develop strategies to quantify the acoustic reliability of BS measurements, in particular related to the environmental sources but also addressing the system sources when natural reference areas are not at hand. In addition, previous studies have not used repetitive BS measurements to monitor the evolution or impact of human interventions, such as nourishments, on a natural system.
This study aims to exploit the potential of MBES data, providing both bathymetry and BS, to monitor underwater nourishments. To this end, we acquired a time series of seven MBES BS datasets in the tidal inlet of the Wadden Sea island Ameland (The Netherlands) within the period from April 2017 to May 2019. We address the acoustic reliability of the MBES measurements and investigate to what extent we are able to acquire and process a MBES BS time series in such a way that the resulting maps represent the actual seabed properties. At the same time we explore the value of MBES BS time series, in addition to MBES bathymetry, for studying the morphodynamics and distribution of the nourished and native sediments.

Geological setting
Ameland is one of the back-barrier islands in the north of the Netherlands (Fig. 1a). The islands are separated by tidal inlets, connecting the North Sea and the Wadden Sea (back-barrier tidal basin), with flood and ebb tidal deltas on the land and seaward side, respectively. Tidal channels incise the flood tidal delta and are connected to the tidal inlet and the ebb tidal delta. The Ameland inlet extends 6 km into the North Sea with a maximum water depth of 28.5 m and a local tidal range around 2 m. The long-term development of the Ameland inlet results in a period of sediment loss and coastal retreat (Israel and Dunsbergen, 1999). The study area is located around the tidal channel, called Borndiep, in the eastern part of the tidal inlet and is about 80 m south-west from the coastline of Ameland (Fig. 1c) with a size of 3.8 km × 0.5 km (Fig. 1d). The investigated area comprises a part of the tidal channel and a very steep retrograding channel wall. The main sediment types in the area range from sandy mud to medium sand, with a varying amount of shells and gravel. Since 1947 revetments are being placed close to the shore to prevent the tidal channel from migrating landward (Vermaas et al., 2019). The nourishments were conducted in the central and south-eastern part of the study area, close to the channel wall (Fig. 1d).

Acoustic data acquisition
The MBES data acquisition was carried out with a Kongsberg maritime EM 2040C dual-head, operated with the data acquisition software SIS and installed on the survey vessels Siege and Amasus. A Trimble SPS851 GNSS receiver provided real-time kinematic positioning with a horizontal and vertical accuracy of 1 and 2 cm, respectively. The motion sensor was a IXSea Octans gyro. Up to five CTD (conductivity, temperature and depth) measurements per survey were taken with a Valeport sound velocity probe. Table 1 reports the sonar settings used to operate the MBES. The MBES datasets are temporally separated by a minimum of 10 weeks and a maximum of 26 weeks (Table 2). A total of 17 track lines were acquired in NW-SE direction with an overlap of about 120% for adjacent track lines. The navigation equipment, sonar settings as well as the location, number and direction of survey lines were consistent among the different surveys, except for the first survey (April 2017) ( Table 1). The survey vessel (Siege and Amasus) and the MBES were changed after the first 5 surveys (April 2017 to September 2018: Siege; February to May 2019: Amasus). For data control additional survey lines with different azimuthal directions (∼45 • interval) covering the same reference area (called ''compass rose'' pattern in Lurton et al., 2018) as well as calibration lines were acquired. During three surveys in April 2017, September 2018 and May 2019 a rough sea state (wave heights up to 1 m) yielded to a vessel roll of up to 5 • , 12 • and 5 • and a vessel pitch of 8 • , 3 • and 6 • in the north-western part of the survey area, respectively.
The seismic data were acquired in November 2017 with an Edgetech X-star 3200XS sub-bottom system consisting of the SB512i tow fish. A chirp signal with frequencies from 0.5 to 7.2 kHz and a pulse length of 30 ms were used. After filtering, the emitted signal was reduced to a near-Klauder wavelet with a center frequency of 3.5 kHz and an effective bandwidth (between the −6 dB points) of 2.5 to 5.5 kHz. The resulting vertical resolution is about 25 cm. We ran seismic lines mainly in NW-SE direction with a distance of 20 m between adjacent lines. In Fig. 1 only the seismic line presented in this study is shown. The data was post-processed with the Delph software suite and the seismic reflectors were manually interpreted with the support of the vibrocores (Vermaas et al., 2019). The time-depth conversion was performed assuming a sound velocity of 1500 m/s for the water column and 1600 m/s for the subsurface.

Ground truth data
Bed sample locations were selected after collecting the MBES data to sample the variations in the BS level observed in the preliminary BS mosaics. A cylindrical box-core sampler was used, with a diameter of 30 cm, which retrieves to a certain extent the undisturbed seabed. We conducted two sampling campaigns in May 2017 (24 samples) and May 2018 (22 samples) ( Table 2). Each box core was photographed and surface samples (0 to 2 cm) were taken on board for grain-size analysis in the lab. Since the majority of the particles larger than 2 mm were composed of shells and shell fragments, particles larger than 2 mm were sieved out and dry weighted. Grain-size distributions of the mud and sand fractions were determined by laser diffraction (Malvern 2000), a technique that measures the scattering of light. Samples were neither prepared nor peptized and measured without ultrasone conditions.
In the interpretation of BS data in this study, the median grain size, 50 , of the mud and sand fractions (volume percentage) is used and we consider the fraction larger than 2 mm (weight percentage) separately. The 50 value ranges between 220 to 280 μm, with one sample having a value of 20 μm. The amount of particles larger than 2 mm (shell and gravel) varies between 0 and 18%. Sediment cores were taken with a vibrocore in May 2018 and described, photographed and analyzed in the lab. Sediment layers and grain-size distribution were determined similar to the box core samples.

Nourishments
Underwater nourishments were performed between June 2017 and February 2018. The sediments for the nourishments were extracted from an area to the north-west of Ameland (Fig. 1c). Box core and vibrocore samples were taken at several locations in the extraction area before the nourishments started. In general, the surficial and subsurface sediment consists of fine to medium sand with a moderate amount of silt and shell fragments in the extraction area. In the subsurface (50 to 100 cm) also traces of clay were found. At the different sample locations the 50 value ranges from 205 to 217 μm except for one location where a 50 value of 143 μm was measured. The nourished area is located in the south-east of the study area (Fig. 1d). T.C. Gaida et al.

Multibeam data processing
A combination of the software Qimera and developed Matlab scripts are used for bathymetric data processing. Soundings are cleaned with a spline-filter, vertically referenced to Normal Amsterdams Peil (NAP) and gridded into a 0.5 m by 0.5 m grid using the GIS software packages ARCmap 10.5 and SAGA. The BS processing algorithm, written in Matlab, is implemented as follows. In the first step the BS ''beam amplitude'' data and all necessary sonar and environmental parameters are extracted from the Kongsberg sonar datagrams. The second step is to retrieve the BS strength (in dB per 1 m 2 at 1 m) from the received echo level (in dB re 1 V) as accurately as possible. The external sources related to environmental properties, hardware and sonar settings, contributing to the received echo level, are described by the sonar equation (modified from Augustin and Lurton, 2005;Schimel et al., 2018) where is the source level (in dB re 1 μPa at 1 m), modulated by the transmission directivity pattern (in dB) as a function of the frequency and the transmission angle with respect to the sonar axis.
(in dB re 1 V/ μPa) is the transducer sensitivity, (in dB) is the receiver gain applied by the receiver electronics and (in dB) is the directivity pattern at reception expressed as a function of and the receiving angle with respect to the sonar axis. These hardwareand software-related terms are considered and accounted for during the real-time processing within the Kongsberg MBES (Hammerstad, 2000). How well the sonar manufacturer accounts for and and to what extent the conversion of the from acoustic pressure (in Pa) to analogue values (in dB re 1 V) and finally to digital values (in dB re 1 μPa) is calibrated is not fully known. The following environmental terms require additional post-processing steps for a Kongsberg system. The ensonified area is not only affected by sonar characteristics (i.e., beam aperture, pulse length) and properties of the study area (i.e., water depth and signal travel distance ) but also by the seabed morphology. However, the real-time processing of a Kongsberg MBES assumes a flat seabed and calculates based on measurements from previous pings (Hammerstad, 2000). In environments with a rough seabed morphology, this simplification has a significant effect on the ensonified footprint and the true incident angle and, consequently, on the BS level. To account for both issues, the seabed morphology correction procedure described in Gaida et al. (2018) is applied.
The transmission loss (in dB) depends on the water conditions and the signal travel distance and can be written as where the first term accounts for the sound attenuation in the water column and the second term for the energy loss of the signal due to geometrical spreading. Sound attenuation in seawater, = + , comes from relaxation of dissolved substances and pure water viscosity (Lurton, 2010), (in dB/km) (Francois and Garrison, 1982a,b), and can also be a result of attenuation caused by suspended material, (in dB/km) (Urick, 1948). The coefficient is dependent on temperature, salinity, acidity, pressure and . The coefficient is dependent on the volume and the geoacoustic properties of the suspended material and . The Kongsberg real-time processing considers a constant for the entire water column. To calculate a variable with depth the CTD measurements and the expression for , given in Francois and Garrison (1982a,b), are used. The approach to assess the effect of on the BS level due to additional suspended sediment will be described in Section 3.2. BS mosaics are produced by calculating an averaged angular response curve (ARC) over a moving window of 100 pings and then normalizing the BS level at every incident angle with the averaged BS level between 30 • and 60 • for starboard and port side separately.
To improve the quality of the BS mosaics, the average BS level is interpolated between −30 • (starboard) and +30 • (port). For three surveys (04/17, 09/18 and 05/19) air bubbles near the transducer head generated along-track artifacts in the north-western part of the survey area, which are removed by using an along-track stripe filter. Strong BS artifacts are observed in the northern-middle part of the survey area for the September 2018 measurements. The causes could not be found out and therefore this part is removed from the analysis.

Sediment suspension modeling
Sound attenuation due to suspended material results from the combined effect of viscous absorption and scattering (Urick, 1948). Viscous absorption defines a frictional energy loss caused by the interaction of the propagating sound field with the viscous fluid and suspended solid particles (Richards et al., 1997). Sound scattering represents an energy redistribution by the small rigid particles. Urick (1948) demonstrated that sound attenuation caused by suspended material can be estimated by ] 20 ln (10) (3) , is the particle radius [m], and are the sediment and water densities [kg m −3 ] and is the kinematic viscosity [m 2 s −1 ]. The first term in the squared brackets corresponds to the particle scattering and the second term to the viscous absorption. This equation is valid as long as is small compared to unity and 3 (Hoitink and Hoekstra, 2005). To assess the sound attenuation due to suspended sediments in our study area, we use measurements from a tidal inlet in the German Wadden Sea as an analogue (Bartholomä et al., 2009).

Acoustic sediment classification
In this study, the unsupervised Bayesian method for acoustic sediment classification (ASC) is applied to the processed BS data. The method was first developed by Simons and Snellen (2009), where a detailed description is given, and further developed by Amiri-Simkooei et al. (2009) and Gaida et al. (2018). The main concept of the method is based on the central limit theorem stating that the BS strength from a specific sediment type follows a Gaussian distribution for a sufficient number of measurements. The method employs an optimization technique to fit an increasing number of Gaussian distributions to the measured BS data per beam angle. A Chi-square test is used to define the number of Gaussian distributions yielding to the optimal data fit. Here we apply, following the approach described in Gaida et al. (2018), the fitting procedure to incident angles between 40 and 65 • for starboard and port side to retrieve a statistical robust estimate of the optimal number of Gaussian distributions. When the number of Gaussian distributions is estimated, the boundaries between the Gaussian distributions are determined for three reference angles (in this study between 56 and 63 • ). The boundaries define the acoustic classes (ACs) and the number of ACs define the number of sediment types which can be distinguished based on the measured BS data in the survey area. Based on the percentage distribution of the ACs at the reference angles, the ACs are assigned to the BS data at all considered angles (in this study 15 to 70 • ). The application to each survey and each beam angle separately provides a classification approach independent of system dependent sources affecting the absolute BS level, like different transducer sensitivity per sonar head, remaining beam pattern artifacts or changing transducer sensitivity between different surveys.
As long as the relative BS level per beam angle (e.g., due to sediment suspension or insufficient CTD measurements) is not affected, the BS can be used for classification. By using the resulting classes from each survey, the classification results can be compared between the different surveys. For a detailed mathematical description of the method, we refer to Simons and Snellen (2009), Amiri-Simkooei et al. (2009) and Gaida et al. (2018).
The box core samples are used to assign sediment types to the ACs. Therefore, the AC with the highest occurrence (mode) within a radius of 5 m around the sample is calculated and the correlation between sediment types and ACs is evaluated. Based on the correlation, a sediment map is produced from the classified BS data. The calculation of the BS per sample location uses the median of the BS values measured within a radius of 5 m.

External sources affecting the backscatter level
In order to establish the acoustic reliability of the data, we relate modeled and measured uncertainties of the Ameland BS time series to the external sources affecting the BS level. In this contribution, uncertainty is not strictly considered as a statistical quantity (i.e., standard deviation), it describes rather the variation in BS level caused by the influence of the external source. Table 3 summarizes the potential external sources and shows the action which is executed to evaluate or control the particular external source. In addition, the estimated uncertainty in decibels, and whether or not these uncertainties influence the BS mosaics or the ASC, is listed. In the following subsections, we elaborate on the effect of the water column properties, transducer sensitivity and survey azimuth on the BS level in more detail.

Water column: Variation in salinity and temperature
The variation of the water absorption coefficient, , with depths is relatively small, but shows a strong difference among the surveys (e.g., between March and Sept 2018 of up to 27 dB/km, Fig. 2a). This is mainly caused by large seasonal temperature differences. The variation within an individual survey (between 5 to 6 h) varies for the different surveys. This is indicated by the error bars displaying the lowest and highest absorption coefficient per survey. The causes are probably seasonal and tidal differences in the water mass dynamics, which indicates the importance of regular CTD measurements. Fig. 2b displays the difference in sound attenuation per beam angle considering the lowest and highest absorption for a water depth of 25 m. It illustrates the maximum uncertainty of the BS level related to an insufficient sampling of the water column properties (i.e., CTD measurements) in the study area. Solely the October 2017 and May 2018 surveys indicate a noticeable variation in sound attenuation during the survey. The average values are 0.4 and 0.2 dB and the maximum values are 0.9 and 0.4 dB for the outermost beams, respectively. The processing algorithm accounts for the absorption variation in time and depth. However, it needs to be considered that these results are based on 2 to 5 CTD measurements within a time of 5 to 6 h per survey and the maximal variations of the water column properties might not be captured. Therefore, Fig. 2b indicates a rough estimate of a potential uncertainty of the BS level due to variation in water column properties.

Water column: Sediment suspension
In order to estimate the maximum expected sound attenuation caused by suspended sediment, we use the maximum SSCs ( in Eq. (3)) measured during a tidal cycle (Scenario 1, 35 mg/l) and a yearly cycle (Scenario 2, 120 mg/l) in a tidal inlet between the German Wadden Sea islands Langeoog and Spiekeroog (Bartholomä et al., 2009); an area with similar sedimentary and morphodynamic processes to the Ameland inlet. Bartholomä et al. (2009) measured SSC values from October 2006 to June 2007 using an in-situ laser particle sizer and acoustic doppler current profiler (ADCP) as well as taking water samples from the water column (near water surface to 0.75 m above seabed). We extrapolate the peak values from 0.75 m above the seabed to the entire water column with a slight decrease with increasing distance from the seabed (Fig. 3a). This provides a worst-case estimate of the effect of sediment suspension on the BS measurements. Furthermore, we use a maximum water depth of 25 m, an angular range from nadir to 70 • and a grain size ( in Eq. (3)) of 20 μm (Scenario 1a and 2a) and 250 μm (Scenario 1b and 2b). The grain size of 20 μm displays the lowest 50 value obtained from the box core samples in our study area (i.e., lab measurements). The second value results from the reference area where a high rate of flocculation of the fine-grained sediments was observed with an averaged floc size of 250 μm. Fig. 3b represents the sound attenuation over the beam angles for Scenario 1a, b and 2a, b. A maximum sound attenuation of about 2.3 dB is calculated for the highest SSC and the smallest grain size for the outermost beam angle. The decrease in grain size from 250 to 20 μm results in a significant increase in sound attenuation which even exceeds the increase in SSC from 35 to 120 mg/l. The smaller grain size yields to a decrease in sound scattering but it results in an even larger increase in viscous absorption. The sound attenuation of the most extreme case (Scenario 2a) would induce a significant uncertainty for the BS data. However, assuming that in our study area flocculation of the fine grained suspended sediments exists and that the MBES measurements were not conducted shortly after a storm, it is more likely that peak values for sound attenuation in our study area are similar as in Scenario 1b where the sound attenuation does not exceed 0.1 dB. These values would be significantly lower than the uncertainty of 1 dB related to the intrinsic uncertainty of the MBES (Hammerstad, 2000) and therefore the effect is assumed negligible for the datasets considered in this contribution. However, the modeling results also indicate that SSC can lead to a significant effect on the BS level, and therefore surveying during or shortly after storms but also shortly after the nourishments should be avoided.

Dual-head transducers
In order to identify a discrepancy in BS levels between the two transducer heads of the Kongsberg EM 2040C MBES, we calculate the averaged ARC for both heads using the MBES BS data from the entire survey area. Due to the large overlap (120 %) and high number of track lines (17), we assume that both transducers ensonify, on average, approximately the same seabed.
The averaged BS values over the angular range from 15 to 70 • (Fig. 4), show that (1) the transducer heads on the Siege (five surveys with a deviation between sonar heads from 0.5 to 1.4 dB) are better calibrated to each other than on the Amasus (two surveys with a deviation from 3.0 and 3.5 dB) and (2) a high influence of a new installation and a different MBES on the BS level. The estimated difference in the BS levels per incident angle (residual BS curve) between both sonar heads are used to adjust the BS levels relatively to each other (transducer head cross-calibration). This adjustment mitigates the influence on the BS mosaic whereas for the ASC an adjustment is not necessary because the method is applied to each transducer and beam angle separately.

Survey azimuth
To investigate the influence of the survey line direction (survey azimuth) on the BS level due to organized seabed roughness (e.g., small-scale ripples), we follow the approach developed in Lurton et al. (2018). The approach involves a systematic coverage of a reference area with different survey azimuths and comparing the resulting ARC per azimuthal direction. The BS data from both sonar heads are crosscalibrated using the BS residual curves. Comparing the ARCs of six sailing directions from the October 2017 and May 2018 surveys (Fig. 5), we observe the largest deviation for the angular range between 20 to 40 • with maximum values of 3.4 dB and 4.2 dB in the October 2017 and May 2018 survey, respectively. The maximum deviation averaged over T.C. Gaida et al.

Table 3
Potential external sources affecting the BS level. Uncertainty represents the averaged value over the angular range from 15 to 70 • and the value in the brackets represents the maximum uncertainty within the same angular range. The last two columns indicate whether or not the results are expected to be affected by the specific external source.   Bartholomä et al. (2009), as an analogue to the Ameland inlet. Scenario 1 represents the general peak SSC during a tidal cycle and Scenario 2 represents the maximum SSC measured during a yearly cycle. (b) Calculated sound attenuation due to sediment suspension with a 50 value of (a) 20 μm and (b) 250 μm for a water depth of 25 m and a frequency of 300 kHz based on Urick (1948) and Hoitink and Hoekstra (2005).  The main survey directions were approximately 120 and 310 • for all surveys (i.e., parallel to the shoreline). For these azimuthal directions, we observe an averaged deviation of 0.7 dB and 1.1 dB with maximum values reaching 1.9 dB and 2.2 dB for the October 2017 and May 2019 survey, respectively. That means the effect is lower compared to other azimuthal angles but still for some incident angles the effect of small-scale ripples might mask the actual seabed properties. Using these values as an estimate, one should consider that the orientation of the small-scale ripples might change during the acquisition, which would have the same effect as varying survey azimuths. In general, the averaged values are similar to the intrinsic uncertainty of the MBES (1 dB) (Hammerstad, 2000) and therefore have a minor effect on the BS data and the monitoring purpose in this study. Nevertheless, the results indicate that the survey line direction should be constant during the survey otherwise the effect of small-scale ripples on the BS data becomes more dominant.

Backscatter and bathymetric time series
The MBES bathymetric time series is used to calculate the differences in bed elevation. The intra-nourishment (10/17) and postnourishment (03/18 to 05/19) bathymetries are referenced to the pre-nourishment (04/17) bathymetry to indicate bed aggradation and degradation. The April 2017 bathymetry is referenced to a prenourishment measurement conducted in May 2017 to indicate the natural seabed dynamics. The bathymetric difference maps represent an aggradation from April 2017 to March 2018 of up to ∼10 m vertically (Fig. 6). The post-nourishment maps show a degradation in the nourished area and a further aggradation north-west from the nourished area, indicating erosion, sediment transport and deposition towards the tidal inlet. In addition, large bedforms are extensively formed in the northern part of the nourished area and to the north-west of this.
The BS mosaics show, in general, varying patterns over the entire time series, indicating changing seabed properties (Fig. 6). High BS levels are observed on the landward side (shoreface) whereas lower BS levels are observed on the tidal-inlet side before the nourishments were conducted (04/17). The intra-nourishment measurement (10/17) shows a lower BS for the areas that were nourished just before the acoustic measurements (blue polygon, Fig. 6). This indicates that the rough and hard natural seabed is replaced with smoother and softer nourished material. However, the areas which were nourished three months before the acoustic measurement show high BS values again (green polygon, Fig. 6). Similar observations can be made for the postnourishment measurements, which show higher BS levels for the entire nourished area over time. This indicates that the natural dynamics of the area result in an increase in seabed roughness and hardness of the nourished area. In addition, strong BS variations are observed along the bedforms for all surveys.

Estimation of the number of acoustic classes
The results of the Chi-square test for the seven processed MBES BS datasets (Fig. 7) show that the averaged Chi-square values, including the standard deviation, reach the criterion for an optimal data fit (black dashed line) using five Gaussian distributions nearly for all datasets. The February 2019 dataset shows high standard deviations, indicating varying fitting results for the different angles and sides resulting in an ambiguous estimate between four and seven classes. To comply with the other datasets, five classes are chosen for the classification of the February 2019 BS data as well.

Correspondence between acoustics and ground truth
To label each AC and eventually retrieve a sediment map, the acoustics (05/18) are compared with the ground truth data (05/18). Fig. 8a indicates a positive correlation of the BS with the weight percentage of grains larger than 2 mm. AC 5 corresponds exclusively to a bottom type containing more than 1% shell fragments and gravel. AC 4 shows a spread between very little (0.1%) and a large amount (18%) of shell fragments and gravel. The sample with the very large amount (18%) amount of shell fragments contains also 13% of mud, which is a likely cause for the lower BS value and the classification as AC 4 instead of AC 5 . AC 3 is evenly distributed between no shells and traces of shells and AC 2 and AC 1 correspond mainly to samples with no shell content. Here, traces define scattered shell-fragments with a quantity of less than <0.1%. Considering the 50 value of the mud and sand fraction, Fig. 8b shows that, in general, the sediment is very homogeneous with 50 values ranging from 220 to 280 μm in the study area. The BS does not show a correlation to such a small range of 50 values. Although AC 1 only contains 1 sample, the 50 value of 20 μm (classified as sandy mud (sM) after Folk, 1954) discriminates it from the samples in AC 2 with a 3 dB lower BS value. Here, one has to consider that this sample corresponds to a thin sandy mud layer (i.e., a veneer with a thickness of ∼1 cm) deposited on a sand layer ( 50 = 240 μm). Following theoretical BS calculations, one would expect around 8 dB difference between sandy mud and medium sand (APL, 1994), which might indicate that the acoustic signal is also affected by the sand layer underneath.
A qualitative comparison between the acoustics and box core samples is given in Fig. 9. The box core samples (BC007 and BC018), extensively covered with shell fragments, are located in areas with T.C. Gaida et al.    high BS values and high ACs (AC 4 and AC 5 ). Box core samples (BC005 and BC006), indicating well-sorted fine to medium sand with no shell fragments, are located in areas of low BS and low ACs (AC 2 to AC 3 ). Box core sample BC015, indicating well-sorted fine to medium sand with traces of shell fragments and gravel, are located in areas of intermediate ACs (AC 3 to AC 4 ). The quantitative and qualitative comparisons show that the acoustics are mainly driven by the amount of shell fragments in the very homogeneous environment of sand. The final assignment of bottom type to AC, including the accuracy measurement, is listed in Table 4. This assignment leads to an overall accuracy of 64 %.
The first box core sampling (09/05/2017) was conducted 14 days after the pre-nourishment survey (25/04/2017) due to logistical problems and weather conditions. Since we noticed an agreement between 50 and AC for dynamically more stable areas (coarse sediments and not in the tidal channel) and a disagreement for areas with higher morphodynamics (finer sediments and tidal channel), we decided not to include the box core samples in the present study because of the time difference. Therefore, solely the box core samples from 22/05/2018 to 24/05/2018 in combination with MBES BS data from 23/05/2018 are used for the assignment of sediment type to AC. The resulting assignment (Table 4) is propagated to the entire time series. Although the absolute range in decibels varies between 25 and 30 dB among the MBES BS datasets, a consistent number of ACs is estimated by the Bayesian method for the entire time series. This indicates that the width of the Gaussian distributions and consequently the AC boundaries account for the wider decibel range.

Temporal evolution of surficial sediments
The ASC maps, displayed in Fig. 10, reveal the temporal evolution of the surficial sediments. The pre-nourishment map (04/17) shows extensively distributed coarse sediments (shell fragments and gravel) towards the coastline of Ameland. The highest amount of shells and gravel is found from the middle to the south-east part of the study area (complies with the to-be nourished area) indicating an extensive shell layer. The intra-nourishment map (10/17) shows that the nourished sand ( 50 = ∼220 μm) covers the shell layer almost completely (see area in blue and black polygon in Fig. 10). Within the blue polygon ( Fig. 10), the only small purple areas (AC 5 ) that remain, coincide to some extent with the location of part of the revetments, indicating that the highest AC also corresponds to the revetments. Therefore, it is more difficult to identify the revetments in the pre-nourishment sediment map where an extensive shell layer is located towards the coastline. Further to the south-east, the intra-nourishment map (10/17) indicates a high amount of surficial shell fragments and gravel within the nourished area (green and black polygon in Fig. 10). Considering that the nourishments have taken place two months (green polygon) up to only a few days (south-east of black polygon) before the acoustic survey, we hypothesize that a washing out of the fine sediments and an increase of the coarse fraction (shells and gravel) occurred very rapidly. Similar observations are made for the post-nourishment maps (03/18 to 05/19) within the nourished areas. The extent of the shell layer increases with time (AC 4 and AC 5 ; pink and purple area), indicating a continuous washing out of the fine fraction. Moreover, the northwestern part of the study area seems to have an increase in finer sediments (sandy mud to sand) from 10/17 to 05/18, which coincides with the bed aggradation (Fig. 6). Therefore, we hypothesize that the fine sediments of the nourished material are redeposited north-west of the nourished area. The ASC maps visualize that, after the nourishments have stopped in February 2018, the surficial sediment distribution has approached the pre-nourishment state, indicating a recovery (natural sedimentary equilibrium state) of the nourished area (Fig. 10).
Furthermore, the development of bedforms is observed, most strongly in October 2017 to March 2018 and in the February 2019 datasets (Fig. 6). In Fig. 11, showing the bathymetry and sediment classification from the intra-(10/17) and post-nourishment (03/18) surveys, two distinct zones of bedforms can be identified, with higher amplitudes and longer wavelengths on the landward side and lower amplitudes and shorter wavelengths on the tidal-inlet side. Fig. 12 shows two cross-profiles visualizing the bathymetry and ASC along different types of bedforms. Profile 1 is located in the nourished area and Profile 2 about 50 m south-west of the nourished area. As a reference, the pre-nourishment bathymetry displays megaripples with a wavelength of 5 to 10 m with amplitudes of up to 1.5 m along Profile 1 (not shown here). Within three months (June to October 2017) very high (∼2.5 m) and steep (∼25 • ) megaripples (Fig. 12a) with wavelengths of about 40 m are rapidly developed on the nourishments. From October 2017 to March 2018 the megaripples grow into sandwaves with an increase in height to 3 m and an enormous increase of the wavelengths up to 120 m in relatively shallow waters (10 to 14 m). In Profile 2 the smaller megaripples grow into larger megaripples (Fig. 12b). This indicates that T.C. Gaida et al. the artificially supplied sediment, with the local hydrodynamics, causes a high degree of morphodynamic activity within the nourished area. The second observation, as revealed by the ASC, is the different sorting patterns over bedforms in both areas and time periods. The clearest pattern is observed in Profile 1 (nourished area) for the measurement in March 2018 (Fig. 12a), where coarse sediments occur on the stoss sides (AC 5 : gravel and shell-containing sands), finer sediments towards the crests (AC 2 and AC 3 : sand) and even finer sediments (AC 1 : sandy mud) in the troughs. A similar pattern is observed for Profile 2 of the measurement in March 2018, except that the magnitude in ACs only varies between AC 3 and AC 1 . In Profile 1 of the measurement in October 2017, we observe a contrasting pattern of finer sediments on the stoss and coarser sediments on the lee sides and in the troughs. In Profile 2 of the measurement in October 2018 we do not identify a clear sediment sorting pattern which might be related to an early stage of the development of bedforms or an insufficient resolution of the MBES. The BS data is corrected for the effect of the seabed morphology as described in Section 3.1. To what extent MBES BS can be used to resolve sediment sorting along bedforms depends on the relationship between the bedform wavelength and the acoustic resolution of the MBES (see Section 5.3).

Integration of MBES and sub-bottom profiling
In this contribution, we use the seismic data to investigate to what extent the classified, pre-nourishment seabed sediments, as revealed by the ASC, can be acoustically identified in the subsurface. In the processed sub-bottom profiler (SBP) transect (11/17) within the nourished area (Fig. 13), the acoustic characteristics of the reflector reveal three different zones: (1) in the center, the discontinuous, high-amplitude reflector indicates a shell layer, as corroborated by vibrocores (05/18), (2) in the south-east, a continuous and high-amplitude reflector correlates with a consolidated clay layer, as identified in the vibrocores (Fig. 15), and (3) in the north-west a high-amplitude reflector indicates the possible presence of buried revetments.
When plotting the seismic profiles (Fig. 14), the detected shell layer in the subsurface corresponds to the pre-nourishment seabed (04/17) between 200 and 1600 m with a precision of 0 to 30 cm, which is classified as AC 4 and AC 5 (due to the shell-rich seabed sediments at the time, see also Fig. 10). The agreement shows that the layer thickness of the pre-nourishment shell layer and the impedance contrast with the surrounding sediment (i.e., nourished sand) is sufficiently high to be acoustically resolved by the SBP. Between 1600 and 2100 m the clay layer is mapped by the SBP. Between 1600 and 1800 m the SBP shell layer overlies the SBP clay and the MBES shell layer, while between 1800 and 2100 m no shell layer is detected in the SBP. According to the box core samples taken in May 2017 this area consists of sand with varying amount of shells, confirming the ASC. Although the identification of the clay layer in the SBP is certain based on the validation of the reflector with vibrocores, the MBES data indicate that the shell layer is located at the same depth in this area. The accurate detection of the shell layer with the SBP in that area might be hampered T.C. Gaida et al.  due to the directly underlying consolidated clay, as can be seen in VC11. This might also influence the vertical precision where the clay is interpreted in the SBP data. In addition, the vertical position of the SBP (∼decimeters) is less accurate than of the MBES (∼centimeters) and, furthermore, the reflectors observed might be from a lower part of the shell layer with a higher shell content. It shows that the combination of MBES (BS and bathymetry) time series measurements and SBP data can be used for improved interpretation of the data and are a good addition to each other.
The vibrocores (22/05 to 23/05/2018) confirm the assignment of sediment type to AC, in addition to the box core samples (22/05 to 24/05/2018). In VC08 (Fig. 15), the top of the core comprises sand with just traces of shell fragments and in core VC11 sand with a high amount of shell fragments is found, which matches with AC 3 and AC 5 , respectively (MBES measurement in 24/05/18). The horizontal bar indicates the location of the pre-nourishment seabed where sand with a relatively large amount of shell fragments can be found in VC08. For VC11 the mean depth is located in the consolidated clay layer but considering the positioning uncertainty of the vibrocore, it could also Fig. 13. Seismic section acquired with the SBP. Shell layer (red), consolidated clay layer (blue) and buried revetments (green) are identified in the seismic section. The location of the seismic section is displayed Fig. 1d. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) be located in the shell layer. In that case the pre-nourishment ASC, indicating sand with a high amount of shells (AC 5 ), would coincide with the vibrocore. It shows that the pre-nourishment ASC results are representative for the newly formed subsurface and that they can support the interpretation of the seismic reflectors.

Monitoring seabed evolution using MBES backscatter
Environmental monitoring requires temporally stable measurements to assess changes of the environment over time (Lurton and Lamarche, 2015). Using acoustic BS to study the evolution of underwater nourishments means that the uncertainty of the BS needs to be smaller than the necessary resolution to characterize the status of the environment. For the environmental sources, the maximum averaged uncertainty of 1.1 dB caused by the possible presence of small-scale ripples, is lower than the magnitude of the required resolution of MBES BS to discriminate between, for example, fine, medium and coarse sand (2 dB) or fine sand and silt (5 dB) (APL, 1994). The change of the vessel and the MBES after five surveys induced the largest uncertainty related to the MBES. The difference in BS levels between the sonar heads, due to the change of the vessel and MBES, increased from 1.4 to 3.5 dB. Firstly, this indicates an imperfect cross-calibration between the sonar heads and secondly it shows a varying transducer sensitivity between the different MBES mounted on different vessels. A relative calibration of the MBES using a reference area with a temporally stable seabed as recommended in Lurton and Lamarche (2015) and Roche et al. (2018) and performed in Montereale-Gavazzi et al. (2017) could be used to quantify the transducer sensitivity. However, the Wadden Sea is a highly dynamic environment and therefore it was not feasible to perform for every survey a relative MBES calibration.
Although the strategies for assessing MBES BS uncertainties that we presented in this study, provide a good estimate of possible variation in BS by the system and environmental sources, we did not quantify the effect of air bubbles, for example due to rough weather conditions, or the effect of the varying frequencies (270-330 kHz) caused by the dual-head system and the ping modus (dual swath). A promising approach would be the use of the water column data (entire signal) to approximately estimate the sound attenuation caused by air bubbles. Due to the small frequency separation, we expect a negligible effect on the BS level according to the research on multi-frequency BS (a frequency separation of at least 1 octave is required) (Hughes Clark, 2015; Gaida et al., 2020). Malik et al. (2018) presented modeled uncertainties of BS measurements. They related the major uncertainty sources of the BS level to the seabed ensonified area (1 to 3 dB), absorption coefficient (up to 6 dB), random fluctuation (5.5 dB) and sonar calibration (system dependent). In this contribution, we minimized the uncertainty related to the first source by calculating the along-and across-track seabed slope with a very high-resolution of 1.5 m. The absorption coefficient was corrected in post-processing for the entire water column using up to five CTD measurements. We compensated for the statistical fluctuation of the BS level by averaging the data over a number of pings and beam angles. In addition, the Bayesian method considers the remaining statistical fluctuation by generating Gaussian distributions with a standard deviation representing the statistical fluctuations. Finally, the method is independent of a sonar calibration and applicable to relative BS measurements .
However, the ASC maps still showed in some areas along-track stripes (Fig. 10). An explanation could be the influence of air bubbles or varying survey azimuths due to a rough sea state. In case of the varying survey azimuths the observed BS variations of up to 4.1 dB for some incident angles and azimuthal angles in combination with the large overlap of about 120 % (i.e., each location is ensonified from three different tracks lines) might cause these artificial patterns. Another reason could be that due to the very dynamic environment the sediment composition or the seabed morphology have changed during the acquisition of three different track lines (∼1 h). Variation in tidal flow direction could alter the orientation of small-scale ripples and variation in the tidal flow speed could favor the settling and remobilization of very small sediment particles changing the watersediment interface. Montereale-Gavazzi et al. (2019) showed variation of >2 dB between ebb and flood caused by surficial substrate changes due settling of the mud fraction in a highly dynamic area close to a river mouth. Hence, preferring certain survey times with similar tidal flow conditions, decreasing the length of the track lines and avoiding a rough sea state (here, defined as a vessel roll and pitch >5 • ) might help to reduce the along-track stripes due to BS variations.

Ground-truthing in a dynamic environment for acoustic seabed monitoring
Sediment mapping based on MBES BS requires ground-truthing to assign sediment properties to ACs. The ground-truthing has to be timed as close as possible to the acoustic survey to be still representative for the actual seabed at the time of the acoustic measurements, especially for such a dynamic environment like a tidal inlet of the Wadden Sea. Therefore, we discarded the first box core sampling set and only used the second box core sampling set and propagated the ground truth information from this one assigned AC set to the entire time series via the unsupervised Bayesian classification method (see Section 4.3.2). A similar approach was presented in Montereale-Gavazzi et al. (2017) except that they used a reference area with a stable seabed to relatively calibrate the BS level for each survey. Here, we followed the approach of Snellen et al. (2019) which demonstrated good repeatability of the classification results obtained from different surveys using the same unsupervised method. They showed that the method can be applied to relative BS levels and allows to compare the resulting ACs from different surveys given that the number of ACs is consistent among the time series. In case the number of estimated ACs differs due to a new sediment supply or sorting processes, the comparison of the ACs among the time series would be hampered and an additional ground-truthing followed by a new assignment of sediment type to AC is recommended.
Furthermore, the assignment of AC 1 was based on only one sample. First of all, one sample is statistically not representative, even though the contrasting BS validate a separate AC. Secondly, despite the welltimed sampling with respect to the acoustic measurements, there are still a few tidal cycles in between and the ∼1 cm thick mud veneer may indicate settlement of muddy sediments, which in combination with high tidal dynamics can yield to large BS variabilities (>2 dB) as shown by Montereale-Gavazzi et al. (2019). They interpreted the variabilities as short-term sediment changes due to the deposition and remobilization of the muddy particles. To some extent these observations could cause the moderate disagreement (expressed as an overall accuracy of 64%, Table 4) between the ground truth and acoustic data in our area. Future monitoring campaigns need to further reduce the time difference (and/or acquire more information on the SSC) to allow an accurate comparison between ground truth and acoustic data.
Further reasons for the moderate disagreement between box core samples and the acoustic data might be related to (1) the positional inaccuracies between MBES and sampling device, (2) noisiness of the BS data, (3) and other challenges in estimating in-situ properties from (small) samples in the lab (e.g., underrepresentation of coarse particles) and (4) the uncertainty whether a single sample is representative of a larger acoustically ensonified region (Goff et al., 2000). In this study, we observed that the BS level was highly dominated by the shell fragment content (≤2 mm). One reason is that in our study area the sediment consists of mainly homogeneous sand (220-280 μm) and therefore the shell content is acoustically the main difference. As demonstrated in several studies, small variations in the shell content (or larger particles in general) can highly influence the BS level (e.g., Goff et al., 2000;Ivakin, 2012;Gaida et al., 2019). In particular for shell fragments it is difficult to obtain a spatially representative measurement because samples contain only a few shell fragments (see point 3 and 4) and consequently the correlation of seabed properties (e.g., mean grain size or shell content) to BS data is biased (Goff et al., 2000). A combination of box core samples with video footage (e.g., drop cameras, video tows or ROV/AUV platforms) would extent the detailed physical information from a very refined location (i.e., 30 cm) to a larger area Gaida et al., 2019). However, a high energy environment as the Ameland inlet might yield to inaccurate positioning and low visibility and hinder the application of video footage.
When applicable, video footage would help to further investigate why the revetments were not clearly distinguishable from the shell fragments (both assigned to the highest AC) even though the revetments represent a much harder material. Either the revetments are covered with a thin veneer of sediments decreasing the BS or the acoustic response at 300 kHz from the revetment and shell fragments is very similar. The BS is not only a result of the hardness contrast (i.e., acoustic impedance) but also depend highly on the seabed roughness (Lurton, 2010). It is likely that the shell fragments appear acoustically much rougher then the revetments, which might compensate for the hardness contrast.

Spatial resolution of MBES backscatter
The grain-size sorting patterns along bedforms that we acoustically identified at different locations (nourished and natural areas) and time periods (intra-and post-nourishments) are comparable to other studies using MBES BS over bedforms (Lamarche et al., 2011;Van Dijk et al., 2012;Koop et al., 2019). Lamarche et al. (2011) observed differences of up to 10 dB between troughs (high BS level) and crests (low BS level) of sand waves with wavelengths between 100 and 250 m and related the differences to varying sediment types. Koop et al. (2019) acoustically identified coarser sediments in the troughs and on the crests and finer sediments on the steepest part of the stoss side of megaripples (wavelengths between 10 and 25 m). The latter discussed the trade-off between spatial resolution (ability to spatially distinguish between different sediment patterns) and geoacoustic resolution (ability to acoustically distinguish between different sediment types). However, we also need to consider the relationship between the spatial/geoacoustical resolution and the wavelength of bedforms, ranging from few centimeters (ripples) to tens of kilometers (sand banks) in coastal regions (Ashley, 1990;Blondeaux, 2012).
In order to assess the capability of MBES BS to identify grain-size sorting along these periodic features, we need to distinguish between three scenarios displayed in Fig. 16: Scenario 1: Bedform wavelength ( bedform ) is smaller than the ensonifed footprint area ( ) and therefore the bedform contributes to the seabed roughness. Consequently, the BS level is affected by wavelength, amplitude and orientation of the bedforms . Scenario 2: Bedform wavelength is larger than the ensonified footprint but smaller than the spatial resolution of the bathymetry slope calculation ( slope ) and therefore the true incident angle and the correct ensonifed footprint area cannot be calculated. The BS level is affected by the slope of the bedform. Scenario 3: Bedform wavelength is significantly larger than the spatial resolution of the slope calculation (i.e., slope correction reveals variation of slope over bedform). Hence, the BS level can be corrected for the true incident angle and the correct ensonified footprint area. Only in the third scenario the BS level represents the actual sediment properties.
For the MBES and the settings used in this study (Table 1) following values are obtained: = 4 to 16 cm (for incident angle = 70 to 15 • ), = 12 to 480 cm (for water depths between 5 to 25 m and for incident angles between 15 to 70 • ) and slope = 150 cm (i.e., 3 times grid resolution, which depends on sounding coverage and the beam footprint ). That means sediment sorting in bedforms with a wavelength of ≫150 cm can be detected in this study. Sediment sorting in megaripples with a wavelength from 10 to 40 m (Passchier and Kleinhans, 2005;Koop et al., 2019) or sand waves with a wavelength of a few hundred meters (Passchier and Kleinhans, 2005;Damen et al., 2018) are detectable, whereas ripples with a typical wavelength of tens of centimeters (Blondeaux, 2012) contribute to the seabed roughness (Scenario 1) or degrade the quality of the BS data (Scenario 2). For example, Scenario 3 is valid for the megaripples and sand waves in Profile 1 measured in October 2017 and March 2018, respectively, and the megaripples in Profile 2 measured in March 2018 with wavelengths between 20 and 100 m (see green reference scale in Fig. 12, representing the 150 cm reference). However, the megaripples, visualized by Profile 2 for the October 2017 measurement, are probably too short (lower end with ∼5 m) and fall in Scenario 2. This would explain the more chaotic distribution of the ACs along the profile. Using MBES BS bedform is the wavelength of the bedform, is the instantaneously ensonified signal footprint, is the beam footprint and slope presents the spatial resolution of the slope correction.
to identify sediment sorting in bedforms requires the determination of the three scenarios for the specific MBES and the study area and the comparison to the length-scale of the bedforms.

Sedimentary patterns in a nourished area
Cross-and longshore sediment sorting is often observed along beaches (Blondeaux, 2012) which agrees with the observed sedimentary pattern during the pre-nourishment survey in our study area. In our study, recovery of the nourishment area occurred within a few days and the pre-nourishment state was fully reached after roughly 12 months (11 to 16 months). Considering the mud and sand fractions, the composition of the nourished sediment ( 50 205 to 217 μm) is very similar to the original sediment ( 50 220 to 280 μm) which might yield, in combination with the continuous hydrodynamics, to the original sedimentary patterns. Based on both the degradation-aggradation pattern and BS observations, we hypothesize that a washing-out of the fine sediments of the nourished sand and a redeposition northwest from the nourished area towards the Borndiep channel took place. According to Elias (2018) this channel is governed by an ebbdominated flow, and the outflow of Borndiep would therefore have an seaward-directed sediment transport, supporting our observations. In addition, shell fragments and gravel being less easy to erode remain on the nourished area, resulting in a similar sedimentary pattern compared to the pre-nourished state.
For sorting along bedforms, field observations in the Belgian, Danish and Dutch North Sea have shown varying sorting patterns over tidal sand waves with wavelengths ranging from 100 to 750 m in water depths of 3 to 30 m (Van Lancker and Jacobs, 2000; Anthony and Leth, 2002;Passchier and Kleinhans, 2005). In these studies some locations showed coarser sediments towards the crests while other locations indicated a reverse pattern of coarser sediment in the troughs. Theoretical approaches by Van Oyen and Blondeaux (2009) and Roos et al. (2007) indicated that sediment sorting processes are controlled by a balance between reduced mobility effects, supporting transport of fine grains, and hiding effects, yielding to higher transport rates of coarser grains. During moderate tidal currents finer sediments accumulate at the crest while strong currents lead to a coarsening of the crest (Van Oyen and Blondeaux, 2009), which would explain the contrasting field observations. Further modeling results also indicated that grain-size sorting requires a wide range in grain sizes but the location of the accumulation (crest or trough) of the different grain-size fractions highly depends on the fluid and sedimentary conditions (Blondeaux, 2012). Although the inlet (Elias, 2018). Therefore, the theoretical requirements are given that the sedimentary and hydrodynamic conditions in our study area result in contrasting sediment sorting patterns at different times and locations as revealed from the ASC.

Conclusion
In this study, we applied a time series of seven high-resolution and full coverage multibeam echosounder (MBES) datasets of both bathymetric and backscatter (BS) measurements to monitor underwater nourishments in the Ameland inlet (The Netherlands). On the BS data, the nourished sediment could be distinguished from the natural sediment in the study area. Degradation and aggradation, as determined from the bathymetric times series, and changing sediment classification maps obtained from relative BS strength values, revealed that the original seabed sediments with a high shell content were covered by the sandier nourishment material. Within approximately 12 months, the pre-nourishment shell-rich sediment state was recreated by washing out finer sands and mud, which were deposited north-west of the nourishment, in the direction of the ebb-dominated tidal current. These findings were supported by seismic data, where a discontinuous highamplitude reflector identified the shell-rich pre-nourishment seabed sediments in the subsurface, corresponding to the pre-nourishment bathymetry and corroborated by sediment cores.
Rapid generation of bedforms was observed in the nourished material in relatively shallow water (10 to 14 m), whereby megaripples with a height of 2.5 m and wavelength of about 40 m developed into 3 m high sand waves with wavelength of up to 120 m within 6 months. Grain-size sorting patterns over bedforms were exhibited by the systematic variation of acoustic sediment classes, showing coarse sediments (high acoustic class (AC)) on the stoss sides, sandy sediments near the crests and sandy mud on the lee sides and in troughs of the bedforms. The combination of the artificial sediment supply with a relatively large range in grain sizes and the varying tidal currents over bedforms might cause these phenomena.
Employing an unsupervised Bayesian classification of relative BS strengths, five ACs resulted in all datasets of the time series. This allowed assigning sediment type to the ACs for the ground truth dataset at the time of one MBES survey and propagating the assigned sediment classes to the other surveys in the time series without a system calibration. When comparing the time series of BS measurements, evaluation of the external sources of variation in BS measurements is necessary. In this study, we demonstrated that in our surveys, water column properties and survey azimuth in relation to organized seabed roughness (i.e., environmental sources) as well as survey settings and biofouling (i.e., system sources), with a variation in BS of less than 1.4 dB, do not interfere with changes in BS over time. However, for some scenarios external sources such as sediment suspension, survey azimuth, varying transducer sensitivity of the sonar heads and changing equipment may affect the BS level significantly (up to 3.5 dB) and, consequently, may hamper the use of MBES BS data for environmental monitoring. Keeping these factors constant during the monitoring campaign is therefore important.
This investigation demonstrates that the combination of bathymetric and BS time series measurements from MBES is highly valuable in monitoring the evolution of underwater nourishments. This approach helps to explain morphodynamic and sedimentary processes, such as transport, deposition and sorting, and thereby contributes to understanding of how natural systems respond to anthropogenic disturbances. Sediment composition is a main factor in benthic ecology and thus highly relevant in ecological impact assessments, which are increasingly important in coastal maintenance and its surroundings. In addition, MBES BS can be used to support the interpretation of seismic reflectors and, inversely, sub-bottom profiling can be used to establish the subsurface structure of the nourished seabed. Furthermore, the sediment classification maps can be used as input for hydrodynamic and morphodynamic modeling studies, providing information on the spatial distribution of the sediment types, to improve sediment transport calculations.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

CRediT authorship contribution statement
Timo C. Gaida