A Mega-Nourishment (Sand Motor) Affects Landscape Diversity of Subtidal Benthic Fauna

The Sand Motor is a very large (20 million m3) nourishment constructed along the coast in The Netherlands. The huge volume of sand is redistributed along the coast by natural forces stemming from tidal currents and waves. For environmental evaluation of this large construction, the benthic subtidal fauna has been sampled prior to the construction of the Sand Motor, and at 1, 2, 4, and 6 years after construction. Although some significant differences between years were detected, overall the total density, total biomass and average number of species per sample were surprisingly constant over this time period. However, large differences were found in the species accumulation curves over samples, and in the rank-biomass and rank-abundance plots. These were related to two important trends in the communities. First, the invasive mollusk Ensis leei, the biomass dominant in the years before construction of the Sand Motor, dwindled in importance in later years. Recruitment of the species failed, but it is unclear whether, and how, this is related to the construction of the Sand Motor. Second, the correlation structure between depth, grain size, bottom shear stress due to waves and currents, which is very tight along a linear coast, was disrupted by the Sand Motor. The community composition was shown to depend strongly on these physical factors. The nature of the dependencies did not change, but the range of different combinations of factors after construction of the Sand Motor was widely larger than before. Although samples had similar number of species per sample before and after construction, the average difference between samples after construction was much larger than before. The Sand Motor is a very large construction, leading to loss of a substantial area (order 100 ha) of submarine area, which recovers at a long time scale. Total disturbance of benthos by burial, expressed as area∗(time before full recovery) was shown to be similar for the Sand Motor and for other coastal nourishment schemes when expressed per unit volume of sediment applied. However, in contrast to beach and shoreface nourishments, the Sand Motor led to a habitat diversification in the coastal zone.


INTRODUCTION
Approximately one quarter of the world's beaches is experiencing significant erosion (Luijendijk et al., 2018). As the coastal zone is squeezed between the rising sea level and urban development (Schlacher et al., 2007), future climate conditions and coastal urban growth will exacerbate the coastal erosion problem and increase the need for coastal defense (Temmerman et al., 2013;Tessler et al., 2015). Hard coastal defenses, such as groins and breakwaters, require high capital investment and have several problems, such as inflexibility and scouring (van Rijn, 2011) that can lead to failures or require very high maintenance costs. Their use is increasingly restricted to areas with sand shortage and availability of rocky substrate, while application of "soft" solutions (sand nourishments) is gaining importance worldwide where sand is available and reliable dunes are present (Basco, 1999;Hamm et al., 2002;Hanson et al., 2002). Nourishments became a standard coastal defense method only a few decades ago and different forms and techniques are still being developed. The two most commonly used forms are beach and shoreface nourishments. In beach nourishments, a mixture of sand and water is pumped through pipelines onto the dry beach and redistributed with heavy machinery. Shoreface nourishments are applied in deeper water (typically around 5 m below mean sea level). Per unit of sand volume applied they are cheaper, as the sand can be transported by ship to the dumping zone, but their effectivity (% of the sand applied reaching the beach and dunes) is lower (van Rijn, 2011). In addition, these measures need to be repeated frequently and applied along the entire length of the eroding shoreline. In the Netherlands nourishments have been common practice since 1990. In the first years beach nourishments were used, but nowadays mainly shoreface nourishments are applied. A single shoreface nourishment generally has a volume close to 2 million m 3 of sand. Return times are between 2 and 10 years, depending on local conditions. The coastal zone, including dunes, beaches, and shallow shoreface subtidal zones, is rich in ecological values and ecosystem services (Barbier et al., 2011). These systems are known to harbor unique ecosystems (McLachlan and Brown, 2006) constituted of unique plant and benthic communities, as well as fish, birds and reptiles that depend on the coastal landscape and the benthic resources. Repeated disturbance of these ecosystems is a potential threat, warranting the search for alternative ways to defend sandy coasts.
A sand nourishment can have several devastating effects on the coastal ecosystem (Speybroeck et al., 2006). The macrozoobenthic community is most directly affected, but higher trophic levels may also be influenced as a result (Peterson et al., 2006). Damage to the ecosystem is caused by the construction activities and the physical alteration of the environment. We specify these effects below.
Construction activities may cause visual and audible disturbance to shorebirds and other terrestrial organisms (Peterson et al., 2006) and cause temporal murkiness of the water which could hinder visual activities of animals and hamper the feeding of filter feeding organisms (Essink, 1999). Yet, the most detrimental effect of a nourishment construction is the death of the benthic fauna by burial. The mortality of benthic organisms is usually extremely high, depending on the thickness of the newly added sediment layer, the mobility of the species and their resistance to an anoxic environment (Löffler and Coosen, 1995;Essink, 1999). It is generally supposed that deposition of layers in excess of 50 cm kill most animals in the sediment. Recovery of the benthic community from such disturbance depends on the ability of animals to recolonize the sediment. It is therefore dependent on season of application (Manning et al., 2014), differential mobility of the species as adults and/or pelagic larvae, development time of the species and extent of the nourishment (Essink, 1999). The latter is particularly important for species without pelagic life stages. In general, it is observed that the benthic community of shallow coasts is capable of recovery from such disturbance because the environment is intrinsically dynamic. As an example, during severe storms the morphology of an exposed shore can change significantly overnight. Short-term studies of post-nourishment recovery show that the nourished area is recolonized the next year by mostly juvenile benthic organisms, but full recovery of the community might take longer, as the first colonizers are mainly opportunistic species and some of the dominant species in late successional stages may require several years to reach sexual maturity (Gorzelany and Nelson, 1987;Adriaanse and Coosen, 1991;Peterson et al., 2000;Menn et al., 2003;Peterson et al., 2006;Jones et al., 2008;Fanini et al., 2009;Leewis et al., 2012).
Nourishments may also lead to physical alteration of the environment. If the nourishment changes morphological characteristics, e.g., grain size of the sediment or slope of the coastal profile, recovery may not be possible before these variables return to their original values. An increase in grain size is technically preferred by engineers as it enhances the stability of nourishments (van Rijn, 2011), but results in habitat changes that hamper recovery over many years (Peterson et al., 2000). A high mud content can temporarily change the murkiness, but more permanently change the community structure as species have differential responses to mud content of the sediment (Cozzoli et al., 2013). Furthermore, an alteration of the grain size may also alter the beach slope and morphology and thereby the complete ecosystem (Degraer et al., 2003;Rodil and Lastra, 2004;McLachlan and Brown, 2006). A consequence of these physicalbiological interactions is that apart from the biological time scales determining the pace of recolonization, also morphodynamic time scales play an important role. The present study investigates the ecological impact of a novel nourishment design, the Sand Motor. The Sand Motor is a mega-nourishment pilot that was constructed in 2011 in the southern part of the Dutch coast (Figure 1). With a total volume of about 19 million m 3 , it contains 10 times the volume of a regular shoreface nourishment and takes the shape of a peninsula with a lagoon and a constructed lake. It originally stretched over 2 km of coastline. Waves and currents distribute the sand over the neighboring coast, a process that is expected to extend over 10 km and to last for approximately 20 years (Stive et al., 2013). The pilot is referred to as a "building with nature" solution (de Vriend et al., 2015), where the natural forces redistribute sand alongshore. The Sand Motor will thus protect a much larger stretch of the coastline than the original 2 km that has been nourished. The design was inspired by a natural phenomenon observed at the Dutch barrier Islands. In this system large ebb tidal sandflats slowly migrate toward the adjacent barrier islands, then merge with the coastline and form hook-shaped peninsulas that feed the adjacent shore (Israël and Dunsbergen, 1999). Until 2018, the Sand Motor morphology has developed as predicted (Stive et al., 2013;Luijendijk et al., 2017).
One of the explicit aims of the Sand Motor was to reduce the disturbance frequency of the benthic fauna, compared with regular repeated nourishments (Stive et al., 2013). The expected life-time of the Sand Motor is several decades, thus reducing the need for renewed disturbance approximately every five years. In addition, it was hypothesized that the construction of the Sand Motor would create new gradients and niches, resulting in a more diverse benthic community. The hook-shaped design of the mega-nourishment creates a distinct feature along an otherwise straight and homogeneous coastline that is likely to increase the environmental heterogeneity and hence biodiversity. The benthic community composition, strongly dominated by vertical zonation patterns (Dahl, 1952;Raffaelli et al., 1991;McLachlan and Jaramillo, 1995;Degraer et al., 2003;Janssen and Mulder, 2005), is likely to change with changing wave exposure, flow velocity and sediment characteristics caused by the Sand Motor.
Due to its size, expected long lifetime and scale-dependent alterations to the physical environment, the Sand Motor is a unique "unaffordable experiment" (Peterson et al., 2006). Regular shoreface nourishments are known to shift the area of wave breaking seaward and thereby modify longitudinal and crossshore sediment transport rates (van Rijn, 2011), but the effect is lasting for only a few years and little interference with tidal currents is usually observed because the nourishments are narrow and oriented in the direction of the tidal currents. In contrast, the Sand Motor has already been shown to significantly interfere with tidal currents, contracting the flow lines near the peninsula (Radermacher et al., 2015) and causing recirculating eddies at both sides in the residual flow. This influences nearby sediment sorting (Huisman et al., 2016) and diversifies the correlation structure between depth, waves, tidal current and sediment composition. This diversification could lead to a higher biodiversity if species are able to respond to new combinations of environmental conditions. However, more extreme conditions could also lead to a decrease in biodiversity, as it is known that. only a few, well adapted species can survive in highdynamic environments (Degraer et al., 2003). As a very big experiment, the Sand Motor shares the disadvantages common to such measures (Peterson et al., 2006): there is a lack of replication and moreover control areas are difficult to identify as the extent of the changes in the environment is difficult to predict beforehand. However, in contrast to many studies on smaller-scale disturbance, this unique experiment was studied intensively from a multidisciplinary point of view. Where many studies of shoreface benthos are restricted to correlations with sediment grain size and topology, we have access to the results of a detailed hydrodynamic model (Luijendijk et al., 2017) that summarizes (in its careful calibration) many of the field studies of the last years.
In this paper we address the question whether nourishment effects on benthos are dependent on the scale of the nourishment. If effects are scale-dependent, one can expect that different processes and effects become important as one scales up the operation. In order to address that question, we first describe the spatial and temporal patterns observed in the subtidal macrozoobenthic community surrounding the Sand Motor during the first six years of its existence. We relate these patterns to observed and modeled changes in the environment. We address the question if, and how, biodiversity has changed due to the Sand Motor and whether the scale of the operation has influenced the effects. The analysis of speciesenvironment relations specifically tests the following hypotheses: (1) hydrodynamic forces are the main drivers for community composition in the shallow subtidal coastal fringe and (2) these forces have been reshuffled spatially in their correlation structure, explaining the shifts in community composition and regional-scale biodiversity in the macrobenthos. We use the information from these analyses to contribute to the evaluation of advantages and disadvantages of mega-nourishments as a coastal defense strategy.

Study Area
The Sand Motor is situated along the Holland coast, the slightly curved, 120 km long stretch of sandy shore in the center of the Dutch coastline (Figure 1). The area is wave-dominated with a semi-diurnal tide and a mean tidal range of about 1.7 m. Waves are wind driven and mainly come from the southwest and northnorthwest, with an average significant wave height of 1 m in summer and 1.7 m in winter (Wijnberg, 2002). The subtidal profile is characterized by 1 or 2 sand bars. Details of the design, original morphology and early development are given by de Schipper et al. (2016). Sediment grainsize along the Holland coast is variable, both in space and time. Beach and dune sediments are fine sands (100-200 µm), swash and surf zone (0-8m) have moderate sand (200-400 µm) while finer sands (100-300 µm) occur in the 8-10 m depth zone (Huisman et al., 2016). The sand used for construction of the Sand Motor had an average median grain size of 278 µm.

Macrozoobenthos Samples
The macrozoobenthic community composition, density and biomass were monitored in a wide area surrounding the Sand Motor. Before construction of the Sand Motor in 2010 and after construction in 2012, 2013, 2015, and 2017, the benthic community was sampled by ship during fall using a Van Veen grab sampler with a surface area of 0.1 m 2 . Each year 120 samples were taken, evenly distributed over 12 fixed transects perpendicular to the coastline of 2010, stretching out over 11 km of coastline surrounding the Sand Motor. Starting from 2013, an additional transect was sampled 1 km further southwest and the fourth transect from the south was skipped from the program. The "lagoon" inside the Sand Motor was sampled only in 2013 and 2015 and is excluded from this analysis. Along each transect a single sample was taken at 1. 5, 3, 4, 5, 7, 8, 9, 10, 11, and 11.5 meter below NAP, i.e., Dutch Ordnance Level which is approximately Mean Sea Level (RWS). Positions were predefined based on recent bathymetry measurements. A sediment sample was taken from each grab sample before sieving and stored in a freezer. The remainder of the sample was sieved over a 1 mm mesh size sieve and stored in 4% formaldehyde to be analyzed in the lab. Macrobenthic organisms were identified up to species level where possible and the ash free dry weight (AFDW) was measured for each species per sample. A field campaign in 2011, directly after the construction of the Sand Motor unfortunately failed. Only half of the 2010 samples were analyzed; the samples were selected evenly over the area to ensure comparability with other years. In all years it turned out especially hard to take a sample in the surf zone (1.5-3 m under NAP), resulting in an under sampled depth zone.

Identification and Data Pretreatment
All individuals were identified to the lowest possible taxonomic level. Damaged or very small individuals were sometimes identified to genus, family or higher level, while intact larger specimens could be identified to species level. In order to avoid bias in the diversity measures (e.g., number of species in a sample) higher taxonomic levels were discarded from the counts of number of species whenever a lower taxon belonging to the high-level taxon was found in the sample, but it was counted as a species if this was not the case. All taxon names were checked against the World Register of Marine Species (WoRMs 1 ) and synonymized with the accepted taxon name where needed.
In every year, the sorting of the samples was contracted out to at least two different laboratories, but as this was regulated by commercial tender mechanisms, not all laboratories participated in the sorting and identification during each of the sampling years. Within a year, samples were allotted equally over the sampling depths and transects to the participating laboratories. We checked for differences in identification between laboratories in the following way. For every species we considered all samples of all years and performed a chi-square test calculating the probability that presence of the species was laboratorydependent. This test was, however, approximate at best since some species showed strong increasing or decreasing trends in presence over years, and laboratories were not equally spread over years. Therefore, when a significant laboratory-presence correlation was detected, we used expert judgement to appraise the probability that the correlation was caused by differences in identification practice. For this judgement we took into account the difficulty in identifying the species, as well as the presence of a trend in the species. When doubts on the identification could not be ruled out, species were lumped into a higher taxonomic level (usually genus level) to avoid laboratory bias.

Environmental Variables
In addition to monitoring the macrozoobenthic community, several environmental parameters were measured or calculated for each sampling location.
The sediment samples were analyzed using a laser diffraction Malvern Mastersizer to determine the median grain size and the grain size distribution. In every benthic sample, a subsample (35 ml) was taken to a depth of 5 cm for grain size analysis. Samples were freeze-dried, homogenized and sieved over a 1 mm sieve before analysis by laser diffraction. No other pretreatment of the samples, such as removal of shell remains or organic material, took place.
The bathymetric evolution of the Sand Motor was monitored on a monthly basis via a purpose-built Jetski mounted with RTK-GPS and an echo sounder [vertical accuracy ∼10 cm;van Son et al. (2010)]. This information was used by Luijendijk et al. (2017) to construct the model bathymetries for simulation of hydrodynamics in September. Here we adopted bathymetry from these model bathymetries.
Hydrodynamic parameters (flow velocities and bed shear stress) were derived from numerical modeling using Delft3D (Lesser et al., 2004). The model construction and analysis is detailed in Luijendijk et al. (2017). In short, a 2D tidal and wave model were setup (referred to as ZM model) covering a 10 km area around the ZM. The resolution of the computational grid cells was uniformly set at 50 m. The tidal boundary conditions of the ZM model were derived from a validated large-scale model for the Holland coast. Both the tidal and wave model were calibrated and validated with in-situ measurements, described in Luijendijk et al. (2017). The computed tidal propagation was in very good agreement with the tide gauges at Hoek van Holland and Scheveningen. The computed tidal and wave driven currents were generally in good agreement with two ADCP measurements in the vicinity of the ZM, with correlation coefficient (R) values of 0.88 and 0.96. The model tends to slightly underestimate the peak tidal velocities with deviations up to 10%. The wave propagation was computed with the wave model SWAN (Booij et al., 1999). The computed nearshore wave heights, periods and direction were validated with the measurements from a nearshore wave buoy. The significant wave height showed good reproduction of the observed wave heights (R = 0.96).
For each sampling period a separate one-month simulation was conducted to obtain the tidal and wave-driven current velocities on the measured bathymetry corresponding to the sampling period. For each simulation the observed hydrodynamic forcing conditions (tide, waves, and surge) in September 2011 were applied. The wave fields were computed every 20 min during the 1-month simulation and fed back to the hydrodynamic model to resolve the wave-driven currents. From these coupled current-wave simulations the bed shear stresses were extracted. All model output variables used are given in the Supplementary Material File "modeloutput.csv."

Statistical Analysis
All analyses have been performed in R3.5 (R Core Team, 2019). Basic analysis of the benthos data have been performed using the package vegan (Oksanen et al., 2019) and packages raster (Hijmans, 2019) and gstat (Pebesma, 2004) for spatial analysis.
A Principal Component Analysis (PCA) of benthic data was performed on the log(x+1)-transformed density data. Only the first axis of the PCA was used in subsequent analyses. The second PCA axis was strongly quadratically related to the first one, suggesting that there was only one strong gradient present in the data. PCA analyses of the data of the separate years have been compared to the results of the overall PCA and showed generally strong similarity. Consequently, all conclusions have been based on the overall PCA.
In order to confirm the robustness of the PCA analysis, results of this analysis have been compared to two other ordination methods: correspondence analysis and nMDS. Correspondence analysis has been based on log(x+1) transformed values, whereas nMDS has used Wisconsin double standardization and squareroot transformation, with Bray-Curtis dissimilarity. All three methods were applied to the numerical abundance, as well as to the biomass data. Comparison of the site scores of all six analysis shows that for the first axis (Supplementary Information Figure 1), strong but non-linear correlations exist between all methods. Within the methods, biomass-based and abundancebased analyses were also strongly correlated. For the second axis (Supplementary Information Figure 2) such correlations were not found. This confirmed the robustness of the results for the first axis and lends support to the decision not to focus on the second axis.
The scores of the sites on the first PCA axis have been analyzed by multiple linear regression on the environmental factors depth, bottom shear stress and sediment grain size (median grain diameter). In order to linearize the relationship, the scores and environmental factors have been transformed appropriately. Two model approaches were used. In the first, the model was fitted on the 2010 (pre-Sand Motor) samples and used to predict the PCA scores in the subsequent (post-Sand Motor) years. In the second approach, all data were used to estimate the linear model. The best fitting, most parsimonious model was selected based on Akaike Information Criterion (AIC). The terms kept in this model can be considered significant, as the selected model was substantially better than a model containing only an intercept. The fitted model was subsequently used to spatially predict the scores over the entire area, using kriging with linear co-factors. A single variogram model for the different years was used. It contained an anisotropy term reflecting greater spatial correlation in the alongshore direction than cross-shore. The strength of the anisotropy was estimated based on a directional variogram.
In comparing different groups of samples, some variables have been log-transformed to homogenize variances and make the data more normally distributed. A significance level of 0.05 was used throughout the text.

Macrozoobenthic Community
Out of a total of 519 samples 196 unique taxa were identified, with a total count of 67791 individuals. The average log-transformed total macrobenthic density per sample was significantly different between years (1-way ANOVA, F 4,513 = 4.31, p < 0.05). Tukey HSD post-hoc test showed that it was significantly (at p = 0.05) higher in 2012 than in the other years, which were not different amongst each other and also not different from 2010, before the construction of the Sand Motor. The average (log-transformed) total biomass per sample (AFDW g m −2 ) also showed significant differences between years (1-way ANOVA F 4,506 = 7.71, p < 0.05). It was slightly (and significantly) higher in 2012 than in 2010, peaked in 2013 when the average individual weight was maximal, and then decreased again in 2015 and 2017 to values not significantly different from 2010. The average number of taxa per sample showed significant differences between years (1-way ANOVA F 4,513 = 5.70, p < 0.05). It was slightly lower in 2010 than in the years after construction, although it was not significantly different from the number in 2013. Apart from this slight difference, no obvious time trend in the mean or median number of species per sample was found. However, the ranges show that after 2010 more samples with exceptionally high number of taxa (>30) were found than in 2010, and this may have been an effect of the Sand Motor. These results are shown as box-and-whisker plots in Figure 2.
The most abundant species of the macrobenthic community was the American razor clam Ensis leei. The species made up 17% of the total macrobenthic abundance and 71% of the total biomass (Table 1) In 2017 this species became more dominant in abundance, but also the mollusk species Donax vittatus and Spisula subtruncata became more numerous. In terms of biomass, Ensis leei remained the dominant species in 2015, despite a factor 6 decrease compared to 2013. In 2017 it was replaced as a biomass dominant by Spisula subtruncata, with also Donax vittatus reaching a considerable average biomass.
The species accumulation curves expressing the expected number of species to be found as a function of the number of samples, show a marked difference between 2010 and the other sampling years (Figure 3). In 2010, the expected number of species for any number of samples is significantly lower than in the other years, which differ very little between them. This is not an artifact caused by the lower number of samples in 2010. Random subsampling of the data sets of the other years to the same number of samples as 2010 (62), showed that the average number of species in 62 samples was 84, with a standard deviation of 3. This is much larger than the observed number of 61 species in the 62 samples of 2010. However, the random subsampling suggested that the small differences between the curves of the post-2010 years are not significant, as some relative swapping of the years was observed visually. The difference between 2010 and the following years points to a substantially higher diversity in the years after the construction of the Sand Motor. This contrasts to the impression given by the distribution of number of species per sample (Figure 2c) showing only a very small (or in 2013 even non-existent) increase in number of species per sample after construction of the Sand Motor. The explanation for the difference in species accumulation curves is primarily in the between-sample, rather than in the withinsample diversity.
The Rank Biomass plot ( Figure 4B) shows that the dominance structure in terms of biomass was drastically changed by the near-disappearance of Ensis in the area. It clearly offsets the years 2010-2013 from the years 2015-2017. In the latter years, total biomass was divided over many more species without a clear dominant. The Rank Abundance plot ( Figure 4A), however, offsets 2010 from the other years as a year with a much more expressed dominance of few species. This alters in 2012 and 2013, not due to changes in the dominance of Ensis, and only to a very limited degree by changes in the average number of species per sample as discussed above, but most probably by increased differences between samples. After 2013, the near-disappearance of Ensis probably affected also the rank abundance structure, but more dominance emerged again in 2017 when other species became more dominant.
A PCA indicated one dominant Eigenvector/Axis explaining 38% of the variance. The median value of the scores along this axis ( Figure 5) increased slightly between 2010 and 2012, then steadily decreased over the years to values below that in 2010, but none of these changes were significant. However, a considerable change in the range of the scores is seen after the construction of the Sand Motor: there are consistently more higher values in the years after 2010. This peaks in 2012 and decreases afterward. The lower values of the range are approximately similar over the years.
The first PCA axis (PCA1) correlates positively both with the logarithm of total abundance (R = 0.84, p < 0.001) and with the number of species in the sample (R = 0.93, p < 0.001). Both aspects are summarized in the sum of log-transformed abundances, which correlates extremely well with the first PCA axis (R = 0.95, p < 0.001).
On one side of the spectrum (negative scores on the first axis) we found mobile species that are typically found in shallow and (very) dynamic coastal waters ( Table 2). The most frequently found species were the polychaetes Scolelepis (Scolelepis) squamata and Nephtys cirrosa, and the mysid Gastrosaccus spinifer. On the positive side of the first axis more species with approximately equal loading on the axis were found. Many mollusks, as well as tubeworms dominated this side of the spectrum, and most are known to be sessile suspension feeders. The spatio-temporal distribution of some species with strong negative and positive scores on the first axis is shown in Figure 6.
The species with negative scores were constrained to the surf zone and very shallow subtidal in 2010. After the construction of the Sand Motor, their distribution moves more offshore in many locations, increases in the surf zone, and concentrates in the area Average values across all samples for abundance, resp. biomass of the species are given between brackets. in front of the Sand Motor. The reverse is true for many species with strong positive scores. They move to relatively shallower locations north and south of the Sand Motor, are less abundant in the zone in front of the Sand Motor, and generally become more abundant after the construction of the Sand Motor, although this increase is spatially variable and certainly not true for Ensis. On top of these systematic changes that relate to the Sand Motor and the changed conditions in the environment, there is considerable year-to-year variability in the abundance of most species. This variability is not synchronous for all species, and remains largely unexplained in this analysis.

Relation With Environmental Variables
We obtained spatial maps for the different sampling years of the following environmental variables: depth, grain size distribution, maximum current velocity, mean current velocity, mean bottom shear stress derived from tidal currents and waves, maximum bottom shear stress derived from currents and waves. Of these variables, the first two were obtained by direct observation, the others were derived from numerical modeling. There was some correlation between the scores of the sampling stations on the first PCA axis with each of these environmental variables. The strongest correlation was obtained with depth (R = 0.75, p < 0.001), median grain size (R = 0.23, p < 0.001) and mean bottom shear stress (R = 0.56, p < 0.001). These variables are shown in Figure 7. Note that the axis scale for bottom shear stress has been truncated in order to show the changes over time in the largest part of the domain. Considerably higher values than 1 N.m −2 occur in the surf zone, without any apparent change over the years. The depth maps show the morphological development of the Sand Motor over time. Details of this development are discussed in de Schipper et al. (2016) and will not be repeated here. Overall, sand is eroded from the tip of the Sand Motor and transported to the sides where the shoreface accretes. The erosion at the tip is restricted to the upper part of the profile only. Below a depth of 5-10 m little erosion has taken place. The initially quite steep profile at the tip is therefore flattened to a certain extent. The surf zone around the Sand Motor has become broader and morphologically more active, with temporary sand bars developing that move relatively fast. Apart from these developments, the bathymetric map also shows traces of the foreshore nourishments north and south of the Sand Motor that were applied as the very first step in the Sand Motor construction. They were planned as a precautionary measure in case the Sand Motor would have caused unexpected coastal erosion in the immediate surroundings.  The maps of bottom shear stress and especially of median grain size show that the environmental influence of the Sand Motor goes well beyond the constructed peninsula itself, even when its immediate offshore zone is included. Increased bottom shear stress is found to a considerable distance offshore of the tip of the Sand Motor. Closer inshore, and to the north and south of the peninsula, extensive zones with decreased bottom shear stress have developed. In these zones, there is a rapid transition from these lowered values of bottom shear stress into a widened surf zone with increased wave influence and 2 | Species scores in the PCA. The table shows the 10 species with the most negative scores and the 10 species with the most positive scores, and documents their frequency (% of samples in which the species occurs) and mean abundance (ind.m −2 ).

Taxon
Axis 1  related bottom shear stress. The largest affected far-field area is demonstrated by the median grain size. There is a consistent coarsening of the sediment offshore of the peninsula, where also bottom shear stress has been increased due to the presence of the Sand Motor. North and South of the flanks of the Sand Motor, where bottom shear stress was decreased, a fining of the sediment can be seen. However, in addition to this pattern, extensive very coarse patches have developed in the years after construction. These patches had median grain sizes larger than 600 µm and have largely remained in place during subsequent years. As they developed outside of the direct Sand Motor area, they were not discovered in the research program devoted to grain size developments (Huisman et al., 2016). In fact, the areas several km south and north of the peninsula were originally included into the benthic monitoring program as "reference areas." For this reason, the benthic monitoring had a wider extent than the other monitoring programs. It is difficult to understand the development of these coarse patches. They are consistent in spatial location over several years and are not the result of single deviant samples. They may relate to changes in the tidal currents caused by the Sand Motor, or alternatively may relate to the shoreface nourishments that have taken place just onshore of the location of these patches. Note, however, that the pattern is not caused by the sand used for these nourishments, as the nourishments formed part of the Sand Motor construction project and used exactly the same sand sources as the Sand Motor itself.

Modeling PCA1 as a Function of Environmental Variables
We have tested a large number of different linear models relating PCA1 to different combinations of environmental variables. The most dominant relation is between PCA1 and depth, but the relation is not linear. In order to linearize it, we have transformed the PCA axis 1 scores logarithmically, after adding the maximum negative score. The linear model was fitted on the observations of 2010 and then validated using the other sampling years.
The most parsimonious linear model explaining the scores on the first PCA axis uses median grain size and the depth * bottom shear stress interaction: Details of the model are given in Supplementary Table 1, and compared to the model for a biomass-based analysis in Supplementary Table 2. As was the case for the scores on the first axis of the PCA themselves, the linear modeling of the biomass-based scores was very similar to the results of the abundance-based analysis, and is therefore not further detailed. The linear model was favorably validated on the observations of the subsequent years. The unexplained variance in these years was 28% of the total variance, to be compared to 22% unexplained variance on the calibration year. There was a very small bias in the validation (mean absolute deviation = 0.12 on a range of 2.38). As a check, we fitted the same model on all observations over the full time period. It had an R 2 = 0.72 and virtually the same coefficients as the model fitted on the 2010 observations only. We conclude that the nature of the relationship between the PCA scores and environmental variables was not qualitatively changed between 2010 and the later period.
By adding year as a factor to the model, we tested for any significant year effects. Only the year 2012 was found to differ significantly from the other years, implying that the deviations from the linear model in this year were significantly larger than the deviations in the other years, which were themselves not significantly different from the deviations in 2010. The year 2012 also had the highest abundances of all years (section 3.1). The spatial pattern of the observed and predicted PCA scores was different between the year 2010 and the other years (Figure 8).
For the year 2010 a rather uniform depth-dependent pattern of community composition was predicted. The slight lowering of the predicted PCA scores along the northern and southern borders of the sampled domain are a consequence of the slightly higher bottom shear stresses predicted here. There is also a minor effect of grain size in the middle of the domain. For the following years, first a two-lobed structure is predicted for 2012, with low values in front of the tip of the Sand Motor, and uniformly high values to the north and south of it. From 2013 onward, and mostly influenced by the grain size pattern, this develops into a three-lobed structure, with high predicted values in a patch offshore and north of the Sand Motor, and two predicted highs to the south, interspersed with the area of large median grain size FIGURE 6 | Spatial and temporal variation of the abundance of some species with strong negative scores on the first axis of the PCA (left-hand column) and with strong positive scores (righthand column). The radius of the dots is proportional to the logarithm of the abundance in the individual samples; absence is indicated by the smallest points.
where predicted PCA scores are lower. This is also the pattern shown by the observed scores. There is very little spatial pattern to be discovered in the maps of deviations between model and observations. Only the high values in 2012 stand out as a definite feature (results not shown).
A further insight into the structure of the data was gained by plotting (transformed) PCA1 scores as function of depth for the different sampling years, while indicating median grain size and bottom shear stress by the size and color of the symbols (Figure 9).
In 2010, we see the regular pattern of a straight coastline. Bottom shear stress is high in the shallow surf zone, and relatively low and homogeneous in deeper zones. Some samples have (moderately) higher median grain size than others, at various depths. These samples always show, relative to their depth, the lowest scores on the PCA axis and are thus expected to be poorer in species and abundance.
Construction of the Sand Motor has disrupted this pattern. Whereas bottom shear stress remains high in the surf zone, much like in 2010, we now observe both higher and lower values of bottom shear stress in deeper points outside this surf zone. In general, the benthic community reacts to this bottom shear stress, with higher PCA scores where bottom shear stress is lower and vice versa. This is an illustration of the (highly significant) depth: shear stress interaction term in the linear model of the PCA scores. We also observe lower minima of bottom shear stress in 2012 than in 2010, over several depths. In accordance, also the range of PCA scores is wider in 2012 than in 2010. With respect to grain size, samples with high and very high median grain size occur at a diversity of depths in 2012. Median grain size is significantly related to local bottom shear stress in most years, but not in 2012. When applying the linear model of eq. (1) by year, 2012 is the only year where median grain size has no significant influence on the axis 1 scores of the PCA. In subsequent sampling With respect to median grain size, we still observe high values over the entire depth range, but gradually they acquire a tighter relationship with the PCA scores of the community.

DISCUSSION
In this study, we have used an extensive monitoring database to show that the application of a mega-nourishment, the Sand Motor, has resulted in substantial changes in the benthic community. Although only minor changes in the mean abundance, biomass and number of taxa per sample have been recorded, the total biodiversity of the Sand Motor has increased considerably. We have shown that the major environmental gradient in the area relates to depth, bottom shear stress and composition of the sediment. By fitting a model on the observations of 2010 (before construction of the Sand Motor) and applying it to the data of later years, we show that the fundamental relations between community composition and physical conditions has not changed. No new and previously unimportant physical variables dominate the community composition.
In addressing the question whether the effects of coastal nourishments are dependent on the scale of nourishment, and thus whether a mega-nourishment causes qualitatively different effects compared to smaller nourishments, we have to consider the two types of effects of nourishments on benthic fauna: direct disturbance by burial and indirect effects on biodiversity through modification of the environment. These are specified below.

Direct Disturbance by Burial
The Sand Motor is a very large-scale nourishment operation, in comparison with more traditional nourishments. In most countries, nourishments are restricted to the beach and involve limited (∼10 5 m 3 ) volumes of sediment, although larger volumes are sometimes used for construction of berms or when very long stretches of beach are nourished (Hanson et al., 2002). This practice was largely replaced in The Netherlands by shoreface nourishments at water depths of 5-8 m, involving larger (10 6 m 3 ) volumes of sediment. Compared to these nourishments, the Sand Motor (2 10 7 m 3 ) differs another order of magnitude in scale. Here we discuss whether the ecological disturbance caused by different types of nourishment scale linearly with the volume of sediment applied. Burial disturbance is generally one of the most prominent effects of nourishment activities. Benthic animals do not survive sudden burial with 50 cm or more of sand, causing the area to be void of animal life shortly after nourishment (Löffler and Coosen, 1995;Essink, 1999). Recovery by recolonization takes time and constitutes a temporary loss of biodiversity and ecosystem services. In the case of the Sand Motor, the situation is different from classical nourishments. An area of approximately 120 ha was covered with (on average) a 10 m layer of sand. The final height of this area exceeded the tidal window, effectively turning it into terrestrial habitat. This terrestrial environment on top of the Sand Motor remained hostile to any macrofaunal or plant life for several years, although some primary dune formation eventually took place. From the point of view of benthos, all benthic life originally present over this large area was lost, and recovery in this case does not depend on larval dispersion or other recolonization strategies, but on the morphological lifetime of the supratidal sand mass. In comparison to that long time scale, recolonization of the areas that are re-integrated into the marine realm following erosion of the top layer is fast, as can be seen from the data presented in this paper. Sampling sites in our study were selected based on depth, and thus moved with the submerged area in the part of the Sand Motor where erosion took place (mostly at the tip of the construction). In the years after the construction of the Sand Motor, we did not observe anomalously low abundance or number of species, even in the shallow points that were only recently recolonized. Therefore, in order to estimate burial disturbance, we used the estimated lifetime of the Sand Motor as a whole.
In a comparison of different nourishment types, the area buried does not scale linearly with the volume of sediment applied. Beach nourishments have a typical thickness of 1 m, shoreface nourishments are typically 2 m thick, but the Sand Motor had an average thickness of order 10 m (Stive et al., 2013). Per unit of sediment volume applied, this causes 10 times less area buried than a beach nourishment, and 5 times less than a shoreface nourishment.
However, the time needed to recover also differs greatly between the affected habitats in the different nourishment types. On beaches, typically highly energetic habitats with a characteristic mobile and opportunistic, short-living fauna, recovery is usually obtained within a year, although recovery can be slow or incomplete when inappropriate (e.g., very muddy or very coarse) sediment is used for the nourishment (Speybroeck et al., 2006). Shoreface nourishments in deeper water with somewhat less dynamic habitats compared to beaches, take a longer time for recovery, typically 1-2 years but more time is needed for longer-lived animals that can occur here (Baptist et al., 2008). In the Sand Motor, recovery of the entire buried area is only possible after the complete disappearance of the bulky volume of nourished sediment. This will happen on a time scale of decades, currently estimated at 20-40 years (Stive et al., 2013), resulting in an average recovery time of appr. 10 year. Therefore, the disturbance by burial, when expressed in units of area * (time to recovery), is more or less linearly related to the volume of sediment applied. Summarizing the order-ofmagnitude estimates for area affected and time to recovery given above, beach nourishments have a disturbance of 10 ha * years for 10 5 m 3 ; shoreface nourishments have a disturbance of 100 ha * years for 10 6 m 3 ; the Sand Motor has a disturbance of 1000 ha * years for 10 7 m 3 of sediment. Scaling up the nourishment in terms of volume of sediment applied, does not greatly influence the burial disturbance to the benthic community per unit of sediment volume applied. In that sense, the Sand Motor is not improving on smaller shoreface nourishments and not even on beach nourishments.
In evaluating disturbance, one should also take into consideration the efficiency of the nourishment operations. Beach nourishments bring the sediment directly where it is aimed at (beaches and dunes), while shoreface nourishments only contribute 25 % of the volume to beach and dune volumes (Van Der Spek and Lodder, 2015). In that sense, a beach nourishment has a four times higher efficiency than a shoreface nourishment. The efficiency of the Sand Motor to bring sediment to beach and dunes is still poorly known, as we are not yet far enough in its life-time and much of the sand is still laying where it has been applied. This efficiency, once it has been properly estimated, may determine whether more human profit (in terms of beach protection) has been gained per unit of ecological disturbance than in regular shoreface nourishments. This reasoning could be extended to include more societal costs and benefits of different coastal defense schemes, but this consideration is outside the scope of the present paper.

Effects on Biodiversity Through Modification of Environment
The species accumulation curves (Figure 3) show a marked difference in biodiversity between the year before the Sand Motor was constructed (2010) and the years after the construction. Many species were only found in the latter years, and the expected cumulative number of species in a set of samples is much higher post-construction than before. Two causes, that are not mutually exclusive, can contribute to this difference. Either the species diversity within samples is different, with the samples in 2010 on average poorer in species and more dominated by a single species, or the diversity among samples is different, with samples in 2010 being more similar to one another than in later years. Our analysis has clearly pointed to the latter cause as the dominant trend in the data.
The average number of species per sample has shown that most post-2010 sampling years had only a slight (but usually significant) increase in the average number of species per sample (Figure 2c). However, this analysis also showed stronger spatial differences in species richness than in 2010, as there were more positive outliers in the post-2010 years, with 30-50 species per sample. The very rich samples were located offshore and did not occur right in front of the Sand Motor, but in side areas indirectly affected by it. They were characterized by a lower than average median grain size, and a low bottom shear stress. The species composition of the samples with more than 20 species, showed that they all belonged to the Abra alba community type. Several species were present in over 85 % of the samples with more than 20 species: Abra alba, Fabulina fabula, Lanice conchilega, Magelona johnstoni, Scoloplos armiger, Spio, Spiophanes bombyx, Spisula subtruncata.
It is also unlikely that the Sand Motor has attracted species that were not occurring in the area before. In total 82 taxa (excluding higher taxa) were found in later years but not in 2010. Most of these taxa were rare and their absence could have been caused by chance alone. We calculated the probability of not finding these species in the 2010 samples, using the zero term of the Poisson distribution as a measure of probability, and applying a Bonferroni correction to the critical probability level of 0.05. Only two taxa, Caprellidae and Donax vittatus, passed the probability threshold: they were absent in 2010 but became relatively frequent in the following years. Overall, this provides little evidence for a substantial qualitative change of the regional fauna. The relatively high number of species not found in 2010 is mostly a consequence of the fact that a small group of samples (62 samples in 2010) is compared to a much larger group (457 samples in total in later years).
The effect of the Sand Motor on biodiversity has also not been effectuated through a qualitative change of the relations between species and their environment. The first axis of the PCA represents a very recognizable, often-described gradient between high-energy shore bound communities characterized by low diversity and low density, with species such as Scolelepis (Scolelepis) squamata, Nephtys cirrosa, and mobile crustaceans at one extreme, and a much richer community composed of mostly sedentary suspension feeders with a high dominance of mollusks at the other extreme. Most benthic assemblages occur in between these extremes (Van Hoey et al., 2004). The position of samples along this continuum is very well described based on physical characteristics of the environment: grain size and the depth-bottom shear stress interaction term. We showed that the explanatory model for this relation, established based on the 2010 samples only, was very well capable of predicting the sample scores in later years. Only the year 2012, shortly after the construction of the Sand Motor, was exceptional in that it showed higher abundances of most species, especially species with a high score on the first PCA axis, than expected based on the relations in 2010. The different position of the year 2012 may have been a transition effect caused by the changes brought about by the Sand Motor. However, the overall picture is one of a relatively tight, and fundamentally unchanged, relation between environmental conditions and composition of the fauna. Such relation can be expected based on other studies that have measured or estimated physical hydrodynamic stress. Kröncke et al. (2018) describe a community gradients across a hydrodynamic stress gradient off Norderney, Germany. They found similar communities as described in our study, and also documented a tight relationship with hydrodynamic stress along the cross-shore gradient. Other studies have used proxies for hydrodynamic factors, such as morphological units (Holzhauer et al., 2019), or depth and grain size of the sediment (Janssen and Mulder, 2005), but also these studies describe essentially the same communities and the same gradient. At larger scales including deeper coastal waters, the range of communities found widens, but the link with sediment composition remains strong and the communities found near the beach are also part of the continuum of communities along the stress gradient (Breine et al., 2018).
Summarizing, the increase of biodiversity after the Sand Motor construction, as shown by the species accumulation curves, was not caused by a much higher number of taxa per sample (although a small but significant increase was observed between 2010 and 2012, 2015, and 2017), nor by a qualitative jump in the regional species pool or a fundamental change in the relation between species and their environment. In contrast to these invariant factors, we observe a larger range in scores of the first axis of the PCA after construction of the Sand Motor, and we find combinations of the important environmental variables depth, grain size and bottom shear stress that did not occur before (Figure 9). The diversification of habitats provides more opportunities for species to have high local abundances in some places. More moderately abundant species find places where they can have a higher abundance than found in the 2010 samples. The Sand Motor has transformed a linear coast with strong and simple correlations between distance from shore, depth, grain size, and bottom shear stress, into a patchwork of subareas where these relations are not necessarily the same as along a linear coast. It has been established globally that habitat heterogeneity is one of the strongest factors determining patterns in biodiversity (Stein et al., 2014). The present case illustrates that general principle.

Scale Effects of Mega-Nourishments
The environmental changes brought about by a nourishment are most probably scale-dependent. In contrast to beach nourishments and regular shoreface nourishments, it has been shown that the Sand Motor has a size that is sufficient to significantly change the tidal currents. Field and model data presented by Radermacher et al. (2015) and Radermacher et al. (2017) show a large tidal eddy on either side of the Sand Motor. Grain size data presented by Huisman et al. (2016) confirm the image, and it is also reproduced by the model (Luijendijk et al., 2017). The large size of the Sand Motor and its shape force the tidal flow to contract around the peninsula, creating areas with high bed shear stress and coarse sediment near the tip of the peninsula and low bed shear stress with fine sediment directly on either side of the peninsula. However, none of these studies described the patches of very coarse sediment on either side of the Sand Motor, at a considerable distance north and south. We do not know very well what has caused these patches of sediment with median grain size well above 500 µm. Possibly, the phenomenon is related to the divergence of the tidal flow caused by the eddies, that may locally result in net offshore transport of fine sediment fractions. This requires, however, closer study to be confirmed. It is also worthwhile to investigate whether any relation exists with the foreshore nourishments applied to the areas north and south of the Sand Motor, as these nourishments are spatially closer to the very coarse patches.
Most changes in the benthic community were observed in the relatively deeper water, beyond the outer breaker bar. This suggests that the wave-induced dynamics in the shallow water overrule any changes caused by the Sand Motor. Organisms living in the shallow water are already selected and adapted to a very dynamic and stressful environment. This observation matches with that (Van Egmond et al., 2018) of who found no clear spatial patterns in the intertidal macrozoobenthic community along the coastline surrounding the Sand Motor, with the exception of the wave-sheltered lagoon. The fact that the Sand Motor mostly influences benthic ecology in the deeper water, and over areas considerably larger than its own surface, is the most important scale-dependent aspect of its environmental effect. Modification of the coastal environment by the indirect effects of a mega-nourishment on tidal flows introduces new features that are not observed in linear coastal systems. The spatial range over which these effects extend is much larger than the size of the mega-nourishment. However, as far as the macrobenthic community is concerned, the scale-dependent changes to the coastal environment lead to diversification of the habitat and the fauna, rather than to loss of natural values.

The Changes in Occurrence of Ensis leei
We observed a strong decrease in dominance of Ensis leei in the Sand Motor area between 2013 and 2015. The species was dominating the community, both in abundance and biomass, prior to this change, but in 2015 and 2017 had lost this status. It is difficult to ascertain whether this change is related to the construction of the Sand Motor. The population of Ensis leei is monitored yearly as part of the regular monitoring of shellfish along the Dutch coast. Although this monitoring revealed yearto-year changes in the density and biomass of Ensis over the past 25 years, there is no sign of a dramatic decline between 2013 and 2015 in the Dutch coastal area in general (Perdon et al., 2019). The distribution maps of Ensis leei show that the area around the Sand Motor does not belong to the core distributional area of the species. Much higher density and biomass is found south of Rotterdam and along the Wadden Islands, as well as along the Holland coast north of the Sand Motor area. Between 1995 and 2010, the population has grown exponentially in the Dutch coastal zone. Since 2010, it fluctuates strongly (differences of a factor 3-4 between consecutive years are found) but without an apparent trend. It is surprising, therefore, that it has decreased so strongly around the Sand Motor. The timing of the strong decrease, appr. 3 years after the construction of the Sand Motor, does not suggest a direct link between the Sand Motor and this population change. It cannot be excluded, however, that the Sand Motor somehow interferes with the larval recruitment of the species. The dispersal and differential survival of larvae is an important part of the population processes of bivalves (Huxham and Richards, 2003). In 2010 the mean individual weight of the species was 32 mg. This increased to 95 and 79 mg in 2012 and 2013, respectively, and after the drop in density to 370 mg in 2015 and 140 mg in 2017. The steep increase in mean individual weight together with the steep drop in density suggests that recruitment may have been much worse than survival and growth, probably already from 2012 or 2013 onward. It remains to be seen whether this change will continue in the future. Any link with the construction of the Sand Motor remains speculative but considering the lifetime of the species of several years, it can also not be excluded. Tulp et al. (2020) describe a response of Ensis leei to experimental shrimp fisheries. Large specimens of the species invaded an experimentally fished area and remained there. This could be a response to competition for space with other species that are more vulnerable to disturbance by fisheries, but the exact mechanism is not known. It suggests, however, that the species is not particularly sensitive to environmental disturbance and can, in contrast, profit from it. Witbaard et al. (2015), in their study of the reproduction and recruitment of the species close to the shore of Egmond along the Holland coast, found no response of the species to a shoreface nourishment close by the sampling location. In a follow-up study, Witbaard et al. (2017) found a relation between the condition index of Ensis and sediment grain size distribution, in particular the silt content of the sediment. Using an experimental approach, they demonstrated that the causal mechanism is not a response of Ensis to an existing silt content, but rather an active role of the species in the silt incorporation into the sediment. This suggests an active adaptability of the species and does not provide an explanation for its strong decrease near the Sand Motor.
A remarkable aspect of the community dynamics after the sudden decrease in abundance of Ensis leei is that the average total biomass of the macrobenthic community has shown very little change, despite the disappearance of the biomass dominant. After the relative peak in total biomass in 2012 and 2013, caused by the large individuals of Ensis leei, there was a relative dip in 2015 followed by a substantial increase in the biomass of the mollusks Spisula subtruncata and Donax vittatus. Together with Ensis leei, these three species make up approximately the same average biomass as Ensis leei alone in 2010. The sequence suggests that, even in this very dynamic and stressful coastal habitat, resource limitation may play a role in structuring the community. The recruitment failure of Ensis leei may have created an opportunity for other species to take a more dominant role in the community. In their study where the Ensis density increased in fished areas, Tulp et al. (2020) found no change in total density between the situations with high and low Ensis density. Unfortunately, no biomass data are available in this study.

Species-Environment Relations
The main purpose of the species-environment statistical model in our analysis was not to provide detailed predictive habitat suitability models at species level. For such purposes, nonparametric approaches may be better suited (Reiss et al., 2015). Here, the emphasis was on testing whether hydrodynamic forces are the main drivers for community composition in the shallow subtidal coastal fringe and, more specifically, how the Sand Motor has reshuffled the causal factors (hydrodynamics, grain size distributions) spatially, so as to increase the spatial habitat diversity. Hydrodynamic variables have previously been shown to be important determinants of macrobenthic community structure (Ysebaert and Herman, 2002;Ysebaert et al., 2003;Cozzoli et al., 2017), but the use in habitat suitability modeling is restricted because of restricted availability of models at the correct resolution (Reiss et al., 2015). With a thoroughly validated hydrodynamic model, we did not suffer from this disadvantage. Nevertheless, we were unable to fully describe the variability in benthic community structure using only hydrodynamic variables such as current velocity magnitude or bottom shear stress (separately for currents and waves or summed). In addition to these variables, it was clear that depth played a major modulating role. Similar values for bottom shear stress at larger depth had different influence than in the shallow surf zone. This corresponds partly with a shift from current-dominated (deeper) to wave-dominated shear stress, but the effect was not fully described by using these separate values. Within the surf zone, communities are adapted to high physical stress, almost independently from exactly how high this stress is, whereas in deeper waters much more subtle reactions to (relatively smaller) variations in bottom shear stress are observed. The general pattern (Figure 9, right-hand side) shows the PCA scores to increase with decreasing stress within each of the depth zones, but the absolute responses to change with depth. The mechanistic way in which depth plays a role remains unresolved, despite the availability of many data and the experimental manipulation of the correlation structure.
The relation between depth, bottom shear stress and median grain size is strongly disrupted by the construction of the Sand Motor. It takes 5-7 years for the expected relations (grain size increases with bottom shear stress; community responds in its species composition to grain size) to re-establish. Especially in the first years after the construction of the Sand Motor, community response to grain size as such was weak. At the same time, no relation was found between bottom shear stress and grain size composition. In later years the relation between grain size and community composition (PCA scores) re-established. The mutual correlation between bottom shear stress and grain size composition remained weak (correlation coefficient r < 0.3, significance levels around p = 0.05 for the different years), as it also was in 2010. With respect to predictive capability, we showed that a model fitted to the 2010 situation was able to capture the general evolution of the system after the strong influence of the Sand Motor. Some transient phenomena were not entirely captured as discussed before, but the general pattern was predictable. However, it is also clear that for similar values of the first PCA axis, actual species composition may differ. Such was the case with the replacement of Ensis biomass and abundance by other species.

Design Considerations
The Sand Motor was designed based on a similarity with naturally advected shoals in the Wadden Sea, but in comparison with these phenomena also had some differences. The height of the nourishment, well above the high-water line, has created relatively large areas that remained uninhabited by macrofauna and vegetation for very long periods. In part this causes the high estimate of "burial disturbance" as used in our analysis. This effect could have been avoided if the nourishment had been made lower, leading to a faster re-integration with the subtidal and intertidal marine area. The construction of a lagoon and a (freshwater) lake within the Sand Motor has not been proven very successful from an ecological point of view. The lagoon had too limited water exchange with the sea and has turned anoxic very soon, which was the reason why benthic sampling in the lagoon was stopped after 2012. Improvements of the design to better connect the main body of the mega-nourishment with the marine ecosystem seem possible and promising, although any changes should be evaluated thoroughly with respect to morphodynamic development before they can be implemented.

CONCLUSION
In this paper, we have analyzed the effects of a mega-nourishment on gains and losses in terms of natural (biodiversity) values of the macrobenthic community. We have shown that, in terms of direct ecological losses by burial of habitats and organisms, different nourishment strategies have comparable effects that scale linearly with the volume of sediment applied in the nourishment. In that sense, mega-nourishments present no obvious advantages compared to other forms of nourishment, and final judgement of the ecological implications will have to wait for the quantification of the balance of losses (burial of habitats and organisms) to gains (amount of sediment effectively contributing to coastal defense). It is too early to evaluate the gains of the Sand Motor, as the effectiveness to coastal defense will have to be judged over the life-time of the nourishment.
Mega-nourishments have scale-dependent effects that lead to qualitative differences in the correlation structure between depth, bottom shear stress and grain size. In the linear coast the cross-shore depth profile is similar at all points along the shore, and with depth also grain size and bottom shear stress vary predictably across the profile. A mega-nourishment changes this simple linear habitat structure. Different cross-shore profiles develop between the tip and the flanks of the Sand Motor. Bottom shear stress due to currents is more variable than is the case along a linear shore. A much larger range of sand grain size is found than in the original situation. In addition, the four factors depth, bottom shear stress due to currents, bottom shear stress due to waves and sediment grain size distribution, now occur in combinations that were unseen before. These large-scale effects have led to a diversification of habitats and a correlated diversification of macrobenthic communities. Similar enrichment of the habitat is not expected to occur in smaller-scale nourishments and can be considered an advantage of meganourishments.
The Sand Motor was the first attempt at the construction of a mega-nourishment. If the idea is to be implemented again in the future, the design could be improved compared to this attempt. Such improvements will obviously take into account many more aspects than macrobenthic biodiversity, but experiences gained in this study may contribute to critical evaluations of merits and disadvantages of alternative designs.

DATA AVAILABILITY STATEMENT
All data used for this study are publicly available through Rijkswaterstaat-Informatiehuis Marien. The code and auxiliary data files used in the data analysis are available from 4TU repository (https://data.4TU.nl).

AUTHOR CONTRIBUTIONS
JW was responsible for sampling and sample treatment. AL provided hydrodynamic model calculations. PH, JJM, TY, and JW performed statistical analysis. PH and JJM wrote the text. All authors corrected the text and contributed to the discussion.