Characterizing bedforms in shallow seas as an integrative predictor of seafloor stability and the occurrence of macrozoobenthic species

In soft‐bottom marine ecosystems, bedform variation is induced by wind‐ and tidal‐driven hydrodynamics. The resulting megaripples, sand waves and sandbanks form a spatially and temporally heterogeneous seafloor landscape. The strong physical forces imposed by the migration of these bedforms are important determinants for the occurrence of different macrozoobenthic species. Quantifying the effect of these forces can help in differentiating natural‐ and anthropogenically induced physical stressors. However, large‐scale mapping of seabed morphology at high resolution using multibeam echosounder is challenging, costly and time‐consuming, especially in shallow seas, prohibiting wide swaths. Instead, their bathymetry is typically studied using single‐beam transects that are interpolated to bathymetric grids with a relatively coarse resolution (20 m). However, this leaves out information on smaller scale (<20 m) bedforms that can be ecologically relevant. In the Dutch Wadden sea, a shallow tidal system, we characterized bedform variation at high resolution using single‐beam data for the first time. We calculated a 2‐D Terrain Ruggedness Index (TRI) at sub‐meter resolution along the single‐beam transects and interpolated the results to a full 3‐D grid. We then validated the result by relating TRI to independently modeled hydrodynamic parameters and to the distribution of macrozoobenthic species. We found that TRI successfully integrates the variation of tidal‐driven bed shear stress and wave‐driven orbital velocity. In addition, we found TRI to be a good predictor of the occurrence of macrozoobenthic species. The inferred small‐scale bedforms provide valuable information for separating the relative importance of natural dynamics versus anthropogenic disturbances such as dredging and bottom trawling activities. We discuss that by repurposing already available single‐beam data in this way, bedforms can be characterized at high resolution without the need for additional equipment or mapping campaigns, yielding novel input to decision‐making on marine management and conservation.


Introduction
Coastal morphodynamics studies the complex interplay between hydrological properties, seafloor dynamics and biological and anthropogenic processes in the coastal zone (Wang et al., 2012). Coastal ecosystems are dynamic and diverse systems of large ecological and societal importance (Barbier et al., 2011;Costanza et al., 1997) and are among the densest populated areas worldwide, resulting in extensive exposure to anthropogenic stressors (Harley et al., 2006;He & Silliman, 2019). In soft-bottom ecosystems, benthic marine macrofauna are important determinants of marine biodiversity and a range of ecosystem services (Snelgrove, 1999). In these systems, the relative importance of natural dynamics versus negative additional impacts of bottom trawling fisheries, dredging for shipping and sand supplements for coastal protection on biodiversity and ecosystem functioning have been the subject of much debate (Pitcher et al., 2022;Rippen et al., 2021;van Dalfsen et al., 2000). Anthropogenic disturbances on top of natural dynamics may significantly alter benthic communities. Mapping the key natural environmental determinants of macrobenthic communities is therefore important in understanding the status and restoration potential of these ecosystems, contributing to sustainable marine spatial planning. This requires a clear separation of anthropogenic and natural determinants of biodiversity and ecosystem functioning.
Species distributions are strongly linked to underlying abiotic conditions. The potential occurrence of a species is dependent on their response traits in relation to the possible physiological conditions they can withstand (Kraft et al., 2015). On a large spatial scale, it is generally accepted that with increasing depth, differences in light, temperature and salinity are important physiological constraints for the occurrence of benthic species (Snelgrove, 1999). However, on a more local scale, variation in hydrodynamic conditions and sediment dynamics are an additional important determinant of species distributions (de Jong et al., 2015;Rosenberg, 1995;van der Wal, Lambert, et al., 2017a). In shallow near-shore zones of softbottom systems, differences in hydrodynamic conditions lead to variation in seabed morphology through a combination of near-bed flow processes, wave-driven orbital velocities and sediment sorting and transport (Barman et al., 2016;Dodd et al., 2003;Morales et al., 2015).
Different bedforms thus reflect variation in natural hydrodynamic conditions and in the nature and intensity of sediment transport (Barnard et al., 2011;Lefebvre et al., 2016). Bedforms in turn impose strong selective forces on their associated fauna as with fast moving bedforms, species have to deal with strong variation in burial and exposure (van der Wal, Ysebaert, et al., 2017b). Migration rates of regular bedforms are dependent on their wavelengths (Ernstsen et al., 2005;Wever, 2004) which are highly correlated with wave height (Flemming, 1988). Different sized bedforms thus migrate at different speeds, resulting in different temporal-and spatial scales of ecological impact (Fig. 1). Small bedform patterns, such as sand ripples can move at rates of several meters per hour (Knaapen, 2005). However, the small wave heights of several centimeters make their burial effect negligible leading to low ecological impact. In contrast, large-scale bedform patterns, such as sandbanks and sand waves can reach heights of several to tens of meters (Knaapen, 2005). Nevertheless, their migration rates are only several meters per year, timescales that macrobenthic communities can more easily adapt to. Intermediate bedforms, such as megaripples with wave heights of generally half a meter, are known to move with several meters per day up to a meter per hour (Brakenhoff, 2021;Knaapen, 2005). For macrozoobenthic species, this would mean burial of several decimeters of sediment for hours up to days at a time. Considering that epifaunal suspension feeders and large deep-burrowing siphonate suspension feeders are generally not able to adapt to more than 1 cm of burial due to a restricted mobility (Kranz, 1974) consequences for associated macrozoobenthic fauna. Highly dynamic environments are then generally dominated by mobile and short living, resistant and/or resilient species, while low dynamic environments are often characterized by sessile, long living species that are more vulnerable to disturbances and are generally low-resilient (Holzhauer et al., 2022). By mapping small-scale bedform variation, we can identify areas that are potentially highly diverse and vulnerable systems based on underlying natural dynamics and prioritize these areas for conservation. For offshore soft-bottom systems, the relation between bedform variation and macrobenthic communities is relatively well studied (e.g. Baptist et al., 2006;Damveld et al., 2018;Markert et al., 2015). Some studies indeed show that bedform variation may explain important variation in macrobenthic species distribution in near-shore systems as well (Holzhauer et al., 2019(Holzhauer et al., , 2022van der Wal, Ysebaert, et al., 2017b), but seascape-wide mapping studies are still lacking, especially for the more shallow parts of soft-bottom coastal systems. While the seafloor of deeper parts of coastal shelves and seas is often mapped using multibeam echosounders, this technique is more challenging to apply in more shallow estuaries and barrier island inlet lagoons, as the obtained swath width is limited by the depth, making the mapping of large areas expensive (Landero Figueroa et al., 2021). Single-beam echosounding is a cheaper method for mapping bathymetry but has a smaller coverage and is usually interpolated to larger areas at coarser resolution. For example, single-beam data have been used to create bathymetric maps of the entire Dutch subtidal Wadden Sea every 6 years at 20 m grid cell resolution. However, small bedforms of high ecological relevance, such as megaripples, are not represented during gridding and interpolation of these extensive mapping surveys (Brakenhoff, 2021), complicating the assessment of the relationship between natural bedform occurrence and macrozoobenthic communities. Hence, we aim to (1) determine small-scale bedform variability by calculating a Terrain Ruggedness Index (TRI) based on single-beam survey data and interpolating this to the whole landscape, (2) check if this TRI reflects temporal dynamics by relating it to historic bathymetry data and (3) relate high-resolution bedform data to macrozoobenthic species. For this, we analyzed and combined a single-beam and benthic species survey dataset from the Western Dutch Wadden Sea.

Study area
The Wadden Sea is one of the largest intertidal ecosystems in the world and consists of 33 tidal inlet systems along the Dutch, German and Danish coast (Lotze et al., 2005;Reise, 2005). Barrier islands separate the North Sea from the intertidal flats that are divided up into tidal basins (Elias, 2006). For this study, we focused on the Marsdiep tidal basin in the Western part of the Dutch Wadden Sea (Fig. 2). This is the largest tidal basin (690 km 2 ) of the Dutch Wadden sea (Oost et al., 2018) and is the most western tidal basin situated between the island of Texel and the provinces Noord-Holland and Friesland (Fig. 2).

Terrain ruggedness index calculation
The Terrain Ruggedness Index (TRI) is a measure for terrain variability that can be calculated from bathymetric data (Wilson et al., 2007). TRI is a measure of the variability of depth around a central pixel and is calculated as the average absolute difference of the central pixel with each of its neighbors for a particular window size (Wilson et al., 2007). It is a focal statistics measure, meaning that windows applied to neighboring pixels are overlapping while the original resolution of the data is retained. Here, TRI was derived from the raw single-beam echosounder data also used for the interpolation of the bathymetric map for 2015 by Rijkswaterstaat (Fig. 2) (de Kruif, 2001). To accomplish this, the single-beam data was first rasterized to gridded line transects at 1 m resolution, as done before for the bathymetric interpolation (van Prooijen et al., 2020). For each raster cell, the TRI was calculated using a 55 meter window size with the "tri" function in the spatialEco package (Evans et al., 2021) in R-Studio, Rversion 4.1.0. (R Core Team, 2021). Megaripples, a mesoscale bedform in current-dominated conditions (Hulscher, 1996) -that is larger scale than sand ripples but smaller than sand waves -have typical wavelengths between 0.6-20 m (Idier et al., 2004), but wave lengths of up to 28 meters have been measured for the Ameland tidal inlet in the Wadden Sea (Brakenhoff, 2021). The singlebeam transects likely do not run from crest to trough of each megaripple, as the transects are perpendicular to the direction of the main gulley. Still, the megaripples are not uniform in shape and crests do not line up perfectly parallel. The transects therefore do cross crests and troughs of different ripples but with likely larger wavelengths than 28 meter. By taking a window size of 55 m around a focal pixel, we ensure that at least one entire megaripple can potentially be included in the calculation. To discern the effects of the transects not running perpendicular to the megaripple wave direction, single-beam transects were simulated on available multibeam data for the Marsdiep inlet. For this simulation, the transects were created to run perpendicular to the megaripple formations and parallel to the main gulley (Appendix A). The TRI calculation results in a raster containing 100-200 m spaced transects of equal configuration to the original single-beam echosounder dataset (Fig. 2c). These transects now contained a value of TRI (in m) instead of depth to mean sea level for each pixel, reflecting the mean of the absolute differences of all depths in that window to the mean depth of the window, both measured in meters. For example, a TRI of 5 implies for a focal pixel that all pixels within a 27 m radius ((55-1)/2) around it are on average 5 m higher or lower than that focal pixel. The obtained value was multiplied by 100 and natural log transformed to normalize the statistical distribution of the residuals when used in a regression. Because the result is still an index, this log (TRI*100) will be abbreviated as TRI from here onwards. This TRI value was interpolated to a seascape-covering grid at 20 m resolution using the same interpolation application (DIGIPOL) used for the Wadden Sea bathymetric maps. A 20 meter grid was chosen as this is also the standard grid size that is used for the regular bathymetric interpolation. Since the transects are 100-200 meters apart, a 20 m resolution gives the best mathematical approximation of the unknown bathymetry. DIGIPOL is an iterative linear interpolation method based on triangulation within the MARIA software package (Rijkswaterstaat, RIKZ, 2021). As high TRI values can indicate both small-scale variation in bedforms on large-scale flat seafloor surfaces, as well as steep slopes on edges of gullies (what we are not interested in here), we masked the map at locations with a slope in bathymetry higher than 65 degrees. Slope was derived from the original 20 m resolution bathymetry map using QGIS version 3.18.1 (QGIS Development Team, 2021). Slope and TRI were then compared in predictive power for the occurrences of species to highlight the added value of TRI as an indicator for natural dynamics. Masked areas are not taken into account for the analyses.

Abiotic data
Bathymetry Bathymetric data has been collected since 1985 by the Ministry of Infrastructure and Water Management of the Netherlands (Dutch: Rijkswaterstaat) (de Kruif, 2001). For subtidal areas, this is done by combining single-beam echosounder (NAVISOUND 200) data collected at a resolution of approximately 30 cm along transects spaced approximately 200 meters apart (Fig. 2c). Lidar data of 1 m resolution is used for the intertidal flats. These combined data were interpolated to a 20 m resolution grid (van Prooijen et al., 2020). The echosounding was typically done in relatively calm weather conditions, which may underestimate bedform heterogeneity during stormy conditions. Hence, the obtained TRI values are probably a conservative estimate, and observed patterns might be stronger under more severe weather conditions. On the other hand, bedform patterns may be strongly modified (higher or lower amplitude, different phase and direction) by stormy conditions and then the obtained TRI could be an overestimation. Nevertheless, macrofauna still have to cope with these severe sediment dynamics.

Bathymetric instability
We used the 20 m resolution bathymetric maps derived with the previously described procedure as measured in the periods 1985-1989, 1990-1994, 1995-2000, 2001-2006, 2007-2013 and 2014-2020 to calculate the absolute mean annual changes in depth as an indicator for seabed stability. For the different years of observation, the absolute bathymetric difference was calculated for each period (t) with its consecutive period (t + 1). This difference was standardized to mean absolute depth change per pixel per year for each period by dividing it by the period length in years. Finally, the mean depth change per pixel per year over all periods was calculated (See Appendix B for further details).

Orbital velocity, flow rate and bed shear stress
We related variation in bedform and bathymetric stability to several modeled hydrodynamic variables. Orbital velocity was predicted using the SWAN-Kuststrook wave model using March 2020 wave measurements as a representative for average annual wave conditions (Van Weerdenburg & Vroom, 2021). The root mean square of the amplitude of the predicted orbital movement was interpolated to a 100 9 100 meter grid as representative value for orbital velocity (

Sediment properties
A large-scale sampling campaign was conducted in 2019 throughout the subtidal Dutch Wadden Sea. In this campaign, sampling points were arranged in a 1 9 1 km grid with additional random points (Franken et al., n.d.). Points were sampled for sediment properties and benthic macrofauna. Deeper parts were sampled using a 20 x 30 cm boxcore (0.06 m 2 ). More shallow areas were visited using a Zodiac and four benthos cores (∅ 10.4 cm) were taken and combined at these locations (0.035 m 2 ). At each sampling location, a sediment sample was taken of the upper 4 cm and immediately frozen at À20°C. In the lab, the samples were subsequently freeze-dried for up to 96 h and then homogenized with a mortar and pestle. Prior to grain-size analysis between 0.5 and 5 grams, depending on the estimated grain size, of homogenized sample was weighed over a 2 mm sieve, in 13 mL PP Autosampler tubes. RO water was added and the sample was shaken vigorously on a vortex mixer for 30 s. Median particle size and the percentage silt (fraction <63 lm) of sediments were determined using a Coulter LS 13 320 particle size analyzer and Autosampler. The optical module "Gray" was used for the calculations (Compton et al., 2013). Median grain size (D50) and silt percentage per sample were used for further analysis.

Biotic data
We tested the predictability of TRI for the occurrence of macrobenthic species to see if TRI is a useful indicator of seabed stability and natural dynamics. The species data were also obtained from the sampling campaign in 2019. Macrobenthic samples were sieved over a 1 mm mesh in the field. Larger shellfish were kept separate and stored frozen. All other living organisms were stored on a pHbuffered 6% formaldehyde solution and stained with Rose Bengal (CAS Number: 4159-77-7) for sorting. In the lab, all individuals were identified to the lowest taxonomic level possible, in most cases to species level (Compton et al., 2013). We used presence-absence data of these species for further analysis.

TRI vs. abiotics
The relation between Terrain Ruggedness Index (TRI), as calculated in this study, and various bathymetric and hydrodynamic variables was explored by calculating the Pearson correlation coefficient between all factors using the "cor" function in base R. Prior to this, variables were checked for linear relationships. Bed shear stress, orbital velocity and bathymetry were ln(x + 1) transformed, and slope was square-root transformed to linearize their relationship with TRI. P-values were determined for these correlations using the "cor_pmat" function in the "ggcorrplot" package (Kassambara, 2019). Significant correlations were then plotted using the same package. Subsequently, these correlations were explored through multiple linear regression as well as additive models on the untransformed variables using the "mgcv" package (Wood, 2011). Predictor variables were first tested for collinearity by determining VIF scores (Variable Inflation Factors). Variables with VIF scores higher than 5 were sequentially removed until no collinearity could be detected anymore (Alin, 2010). Flow rate and bathymetry had high multicollinearity with each other as well as with orbital velocity and bed shear stress. Flow rate and bathymetry had highest VIF scores and were therefore removed. In the final model, orbital velocity, median grain size (D50), percentage of silt and bed shear stress were used as predictor variables in a linear model. Additive models did not significantly enhance model predictions (Appendix C) and linear models were chosen for ease of interpretation. The model variables were tested for significance and model assumptions were checked. The best model was selected through stepwise backward selection.
Subsequently, the newly obtained TRI map was related to the bathymetric instability through a linear model. For this, the mean annual absolute bathymetric difference calculated over the last 40 years was ln(x + 1) transformed and used as a response variable in a linear model, checked for assumptions and significance.

TRI vs. species occurrences
To test if the TRI was related to macrozoobenthic species occurrence, we extracted the 20 m-interpolated value from the TRI raster at the geographic location of each benthos core sample. We used the presence-absence data of all identified taxonomic groups that were found in more than 10 samples. We then explored how their probability of occurrence depended on the TRI using generalized linear models with a binomial distribution, testing for a quadratic effect of TRI (reflecting an ecological optimum in probability of occurrence). For each species, the model was tested for overdispersion, in which case a quasibinomial distribution was used instead. From the set of predictor variables, the best model was chosen using stepwise backward selection (so retaining the linear term when the quadratic term was significant). Then, a weighted average of TRI (ecological optimum (Ter Braak & Prentice, 2004)) was calculated for all taxonomic groups that showed a significant relationship between TRI and the probability of occurrence. For this, the TRI value was averaged for all taxonomic groups over all sampling locations where this group was present. Finally, six species that showed low, intermediate and high weighted averages of TRI were further explored. This was done by plotting a regression line based on the significant binomial model determined above.
As mentioned before, high TRI values can indicate both small-scale variation in bedforms on large-scale flat seafloor surfaces, as well as steep slopes on edges of gullies. To show the added relevance of TRI compared with using slope, the same generalized linear models were used with slope as a predictor instead.

TRI calculation
By interpolating the Terrain Ruggedness Index (TRI) values obtained from the single-beam transect data, we were able to create a map of the TRI values for the whole Marsdiep tidal basin (Fig. 3). Even after masking areas where TRI values were high due to the effect of steep slopes of gullies, we observed that TRI values were highest in the deeper gullies ( Fig. 3; Appendix D). These patterns were also observed in the original single-beam transect, high-resolution data, where single-beam transects crossing deeper gullies showed larger small-scale variation in the deeper parts than in the shallower areas (Fig. 3). The 20meter interpolated bathymetry map did not capture this small-scale bathymetric variation, which shows that the calculated TRI values presented in this paper were able to distinguish between small-scale bedforms that are typically not resolved when processing raw single-beam data to bathymetric maps (Fig. 4). Simulating transects that run perpendicular to the bedform migration does not result in major differences of TRI compared with the original transects either (Appendix A).

TRI in relation to hydrodynamic and bathymetric factors
We tested how the derived TRI values relate to other independently modeled and derived abiotic data. The correlogram of the Pearson correlation coefficient showed strong correlations between TRI and the hydrodynamic and bathymetric factors. There was a strong negative correlation between TRI and bathymetry (À0.72) as well as with orbital velocity (À0.75). On the other hand, TRI was strongly positively correlated with bed shear stress (0.76), current velocity (0.75) and slope (0.88). Median grain size was only weakly positively correlated with TRI (0.42), whereas median grain size was strongly positively correlated with bed shear stress (0.63) and current velocity (0.62). Hydrodynamics and bathymetry were likewise strongly correlated, indicating the strong codependence between all variables (Fig. 5a).

TRI as an indicator for the occurrence of species
We found significant relationships between our TRI and the occurrence of 31 out of 50 tested taxa ( Fig. 6; Table D2, logistic regressions). In the majority of cases, TRI was a better indicator in predicting the occurrence of taxa compared with slope (Appendix E), highlighting the added value of TRI in modeling species occurrence and habitat mapping. High TRI values were associated with mobile scavengerpredator species like Nephtys cirrosa and Eumida sanguinea (Fauchald, 1977;Hayward & Ryland, 2017) as well as mobile deposit-feeding species, such as Bathyporeia spp. and Microphthalmus sp. (Fauchald, 1977;Nicolaisen & Kanneworff, 1969) (Fig. 6). Bathyporeia spp. and N. cirrosa, with high weighted averages of 4.38 and 5.13, respectively, showed an increase in occurrence with increasing TRI  Table D2). Less mobile deposit-feeders, such as Scoloplos armiger and highly mobile suspension feeders like Ensis leei (Kranz, 1974;Shumway et al., 1985), showed high occurrence rates at intermediate TRI (weighted average of 3.94 for both) (Fig. 6) Figure 7; Table D2). Sedentary deposit-feeders like Aphelochaeta marioni (Bolam, 2011), suspension feeders like Mya arenaria (Kranz, 1974) and Cerastoderma edule (Andr e & Rosenberg, 1991) are very sensitive to burial (Bolam, 2011), hence their association with low TRI (Fig. 6). C. edule and Peringia ulvae are both associated with the lowest weighted average of TRI of 3.06 and 3.13, respectively and show decreasing occurrence with increasing TRI (binomial regression, B = À1.00, 95% CI = [À1.36, À0.60] and    Table D2).

Discussion
We have outlined in this paper a novel method to determine fine-scale seabed morphology using existing singlebeam echosounding data which can be used to predict macrozoobenthic species occurrence. We find that Terrain Ruggedness Index (TRI) accurately captures small-scale bathymetric variation. Variation in TRI related to both spatial and temporal dynamics of the seafloor and was found to be highly correlated with bed shear stress and orbital velocity. In addition, TRI significantly predicted the occurrence of 31 out of 50 macrozoobenthic taxa that were found in a large-scale sampling campaign. We

Nephtys cirrosa
Nemertea  indeed observed a trophic shift from sedentary depositfeeders and suspension feeders in low dynamic (indicated by low TRI) areas to mobile scavenger-predator and mobile deposit-feeding macrozoobenthic species in highly dynamic (indicated by high TRI) areas. Dynamic areas with high burial rates, such as megaripple areas characterized by high TRI, are expected to be dominated by species with short lifespans and fast reproduction (Holzhauer et al., 2022), mobile predator-scavengers or deposit-feeders (Robert et al., 2021;van der Wal, Ysebaert, et al., 2017b). Detailed occurrence, dynamics and biotic relationships of multiscale bedforms have been studied extensively by multibeam echosounding for a smaller part of the Dutch Wadden Sea, the ebb-tidal delta of the Ameland inlet (e.g. Brakenhoff, 2021;Holzhauer et al., 2019Holzhauer et al., , 2022van Prooijen et al., 2020). Holzhauer et al. (2022) highlight strong responses of the benthic community to morphometric seabed properties. However, high-resolution mapping of shallow soft-bottom systems, such as the Dutch Wadden Sea, to characterize these morphometric seabed properties was not done due to the complications of using multibeam sonar in shallow seas. We show that the TRI we have calculated using existing single-beam data is a good indicator of morphometric seabed properties such as megaripples. Comparing the calculated TRI values to single-beam transects shows clear patterns between TRI and bathymetric variations. Small bathymetric variation, indicating smaller megaripples, found in the shallower gullies under influence of intermediate hydrodynamic forcing, is associated with intermediate TRI values ranging from 4 to 5. Shallow sand flats, with lower hydrodynamic forcing through bed shear stress, show less bathymetric variation, and hence are associated with low TRI values (<4). Contrastingly, deeper gullies that are known for their rough hydrodynamic conditions show the largest variations in bathymetry and similarly have highest TRI values (>5). Despite the fact that the single-beam transects do not line up crossing megaripple formations, we are still able to discern these bedform formations and the strength of their formation. Prior knowledge of the general direction of megaripples (e.g. by adding high-resolution single-beam transects in the direction of flow) would enable more accurate mapping of their effects as well as give information of their specific morphometrics (i.e. wave length, migration speed). Still, by taking a larger window size when calculating TRI, we ensure the inclusion of at least one crest and trough to be able to identify the presence of megaripple formation and gauge the strength of their migration. Areas with high TRI and thus largest bathymetric variation are then hypothesized to have the largest sediment displacement. Species in these areas have to be adapted to constant burial by large amounts of sediment. They need to be highly mobile and resilient enough to withstand constant disturbance, as confirmed by our results.
Using historical bathymetric data, we found that spatial heterogeneity of the seafloor as measured by TRI is a good indicator for its larger scale temporal bathymetric instability with large potential implication for macrozoobenthic species distributions. These dynamics are largely driven by the interplay between sediment and hydrodynamic conditions, such as bed shear stress and orbital velocity. The intensity of hydrodynamic forcing affects rates of sedimentation and erosion, so natural dynamics of seafloor habitats, with largest ecological impacts are expected at intermediate forcing (Fig. 1). Here, we found high TRI to be indicative of high dynamic areas as high TRI was strongly related to hydrodynamic forcing (i.e., bed shear stress and orbital velocity). As this relationship was log-linear, the untransformed fitted models show an exponential increase of TRI (and thus seabed dynamics) with higher bed shear stress. However, it can be argued that with even higher hydrodynamic forcing the system would shift from softbottom to hard substrate with only rocks being able to remain. Consequently, bedform dynamics would get lost and thus only able to be sustained at intermediate hydrodynamic forcing. High TRI was also related to high temporal bathymetric instability. High temporal instability indicates large and fast movement of sediment, resulting in large annual differences in depth through the constant migration of bedforms. The strong relationship between the spatial heterogeneity captured by TRI and the temporal variability in bathymetry emphasizes the quality of TRI as an indicator of ecologically relevant bedforms, such as megaripples. High TRI values thus reflect the massive sediment transport and constant burial that species have to deal with. This imposes a strong selection for species that are highly mobile either through the water-and/or sediment column (Hinchey et al., 2006;van der Wal, Ysebaert, et al., 2017b;Z€ uhlke & Reise, 1994).
The physical forces imposed by hydrodynamics and related migrating bedforms, as captured by the TRI, is expected to be reflected in the functional trait makeup of a macrozoobenthic community. We find sedentary depositfeeders and suspension feeders like Aphelochaeta marioni (Bolam, 2011), Mya arenaria (Kranz, 1974) and Cerastoderma edule (Andr e & Rosenberg, 1991) in low TRI areas, whereas high TRI areas are characterized by mobile deposit-feeders and scavenger--predators like Bathyporeia spp. (Fauchald, 1977;Nicolaisen & Kanneworff, 1969), Nephtys cirrosa and Eumida sanguinea (Fauchald, 1977;Hayward & Ryland, 2017). Other studies on macrozoobenthic communities in response to natural dynamics find similar results. van der Wal, Lambert, et al. (2017a) found a trophic shift from a deposit-feeding community (e.g., Heteromastus filiformis, Scoloplos armiger and A. marioni) in low hydrodynamic conditions to a predator and omnivore-dominated community with high mobility (e.g., Nephtys spp.) in the Westerschelde estuary. Highly dynamic areas in the ebb-tidal delta of the Ameland inlet were characterized by B. pelagica, B. elegans and N. cirrosa (Holzhauer et al., 2022). Contrastingly, species like Urothoe poseidonis and Nephtys hombergii were indicative of low dynamic areas (Holzhauer et al., 2022). Likewise, van der Wal, Ysebaert, et al. (2017b) found high energetic sites with large bathymetric variation to be distinguished by Bathyporeia spp. and Nephtys spp., whereas H. filiformis and S. armiger were mainly found in low-to-intermediate energetic sites. Pygospio elegans and Hediste diversicolor were especially indicative of low-to-intermediate energetic sites (van der Wal, Ysebaert, et al., 2017b). The high degree of similarities between bedform-species patterns observed in previous studies and this study using TRI, further indicates the potential to use this proxy for large-scale bedform characterization for habitat modeling.
Seafloor morphodynamics not only predicts species occurrence, but is also subject to dominant macrobenthic species that influence sediment dynamics. Interactions between the bioturbating activities, or the presence of ecosystem engineers might alter the formation of bedform patterns (e.g. Borsje et al., 2009;Damveld et al., 2020;Malarkey et al., 2015), and thereby provide protection from hydrodynamics and facilitate more sensitive species. Van der Reijden et al. (2021) found patches of Sabellaria spinulosa, an epibenthic reef-forming polychaete, in the troughs of megaripple waves in the North Sea. They hypothesized the S. spinulosa reefs to be capable of reducing megaripple migration to rates these reefs can cope with (van der Reijden et al., 2021). Modeling studies show the two-way coupling of Lanice conchilega reefs and sand wave morphology (Damveld et al., 2019(Damveld et al., , 2020. Hydrodynamic attenuation of other reef-forming species might also influence the formation of bedform structures. Blue mussel (Mytilus edulis) beds are known to greatly influence near-bed flow properties (Folkard & Gascoigne, 2009;van Leeuwen et al., 2010). This may alter bedform formation and migration and facilitate more sensitive species in otherwise high dynamic areas. Discrepancies in TRI and hydrodynamics might be due to the presence of biogenic reefs. Unfortunately, our sampling method was too coarse to properly assess the spatial extent of biogenic reefs to make this comparison. Additionally, boxcoring is not a suitable method to sample hard substrates like bivalve beds as the core will not penetrate the reef. This may have caused an underrepresentation of species in bivalve beds. Nevertheless, low TRI values in areas despite strong hydrodynamic conditions might then hint at the occurrence of sediment stabilizing species and warrant further inspection as these areas might be of high natural value.
Additional anthropogenic disturbances within highly dynamic areas are usually considered negligible (Bergman et al., 2015;Rijnsdorp et al., 2018), but this has been poorly tested until now for shallow coastal waters. Hydrodynamic conditions can shape macrozoobenthic communities along gradients of resistance to physical disturbances (Hershkovitz & Gasith, 2013;van der Wal, Lambert, et al., 2017a;Warwick & Uncles, 1980). High physical forcing in dynamic areas then predispose macrozoobenthic communities to be resistant against additional anthropogenic perturbations (Rijnsdorp et al., 2018). Alternatively, a long history of continued anthropogenic disturbance may have altered benthic communities to be predisposed to withstand disturbances, regardless of their natural dynamics (Bergman et al., 2015). In addition, highly energetic hydrodynamic conditions in soft-bottom systems can give rise to sand wave structures (Passchier & Kleinhans, 2005). These bedforms in turn alter hydrodynamic conditions, creating microclimates with calmer conditions that allow more sensitive species to persist (Hoziaux et al., 2008;Mestdagh et al., 2020;van der Reijden et al., 2019). These communities might then not be able to persist under high anthropogenic pressures. Still, bedform structures such as megaripples are often superimposed on larger bedforms such as sand waves (Koop & Amiri-Simkooei, 2019;Lindenbergh et al., 2006), adding to the spatial complexity of species responses to hydrodynamic conditions as well as their resistance to additional anthropogenic perturbations. Characterizing all natural physical forces allows for the appropriate incorporation of the natural dynamics background when studying the impact of anthropogenic disturbances on ecosystem functioning and biodiversity in shallow coastal ecosystems.

Conclusion
High-resolution marine habitat mapping can not only help in understanding species distributions but also in prioritizing conservation areas by identifying areas with the potential to develop or sustain communities of sensitive species. Here, we have shown a method to repurpose existing single-beam data to map seafloor dynamics on a large spatial scale. By repurposing the data of an existing program there is no need for additional mapping programs that are usually expensive and often logistically challenging. The strong correlation with hydrodynamics highlights TRI as an indicator for natural dynamics and can help in discerning the effects of natural-and anthropogenic dynamics in relation to questions about the impact of anthropogenic activities such as bottom trawling. The obtained information can then be used for conservation and marine spatial planning purposes. Integrating measures of natural dynamics in marine spatial planning and the decisions for placing Marine Protected Areas will aid in a sustainable division of seafloor habitats between nature conservation and anthropogenic use. the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

Data Accessibility Statement
Data supporting the findings of this study are available from the authors and will be made archived and publicly available in the University of Groningen Research Data Repository: https://doi.org/10.34894/PZJSSY.