Biogeomorphology in the marine landscape: Modelling the feedbacks between patches of the polychaete worm Lanice conchilega and tidal sand waves

Tidal sand waves are dynamic bedforms found in coastal shelf seas. Moreover, these areas are inhabited by numerous benthic species, of which the spatial distribution is linked to the morphological structure of sand waves. In particular, the tube‐building worm Lanice conchilegais of interest as this organism forms small mounds on the seabed, which provide shelter to other organisms. We investigate how the interactions between small‐scale mounds (height ∼dm) and large‐scale sand waves (height ∼m) shape the bed of the marine environment. To this end, we present a two‐way coupled process‐based model of sand waves and tube‐building worm patches in Delft3D. The population density evolves according to a general law of logistic growth, with the bed shear stress controlling the carrying capacity. Worm patches are randomly seeded and the tubes are mimicked by small cylinders that influence flow and turbulence, thereby altering sediment dynamics. Model results relate the patches with the highest worm densities to the sand wave troughs, which qualitatively agrees with field observations. Furthermore, the L. conchilegatubes trigger the formation of sandy mounds on the seabed. Because of the population density distribution, the mounds in the troughs can be several centimetres higher than on the crests. Regarding sand wave morphology, the combination of patches and mounds are found to shorten the time‐to‐equilibrium. Also, if the initial bed comprised small sinusoidal sand waves, the equilibrium wave height decreased with a few decimetres compared to the situation without worm patches. As the timescale of mound formation (years) is shorter than that of sand wave evolution (decades), the mounds induce (and accelerate) sand wave growth on a similar spatial scale to the mounds. Initially, this leads to shorter sand waves than they would be in an abiotic environment. However, near equilibrium the wavelengths tend towards their abiotic counterparts again. © 2020 The Authors. Earth Surface Processes and Landforms published by John Wiley & Sons Ltd


Introduction
The underwater landscape of sandy coastal shelf seas, such as the North Sea, is characterised by fascinatingly rhythmic bed features of various shapes and sizes (McCave, 1971;Stride, 1982). Tidal sand waves (Figure 1) are highly relevant to study from an engineering perspective, as their occurrence, size and dynamic behaviour interfere with offshore engineering applications (Németh et al. 2003;Roetert et al. 2017). Sand waves have heights up to 10 m, wavelengths of the order of hundreds of metres and migrate several metres per year (Van Dijk and Kleinhans, 2005;Damen et al. 2018). Moreover, after dredging they regenerate in several years' time (Knaapen and Hulscher, 2002). These offshore areas are also the home of numerous benthic organisms (collectively referred to as benthos) living in and on top of the seabed (Witbaard et al. 2013). The habitat distribution of these benthic organisms is shown to be related to temporal (e.g. seasonal) and spatial variations in the physical environment (Barros et al. 2004;Baptist et al. 2006;Van Dijk et al. 2012;Markert et al. 2015; Van der Wal et al. 2017;Damveld et al. 2018; Van der Reijden et al. 2019). At the same time, benthos may modify their physical habitat (Eckman et al. 1981;Jones et al. 1994;Rabaut et al. 2009), through which they shape the landscape around them on a wide scale.
These links between organisms and their physical environment are generally referred to as biogeomorphology (Viles, 1988;Murray et al. 2008;Reinhardt et al. 2010;Corenblit et al. 2011). Biogeomorphological processes shape a wide range of landscapes, for instance salt marshes (Marani et al. 2007), mangroves (Friess et al. 2012) and coastal dunes (Baas, 2002). Also, in subtidal areas these processes can leave distinct marks in the landscape. Field observations show that the highly abundant tube-building polychaete worm Lanice conchilega ( Figure 2) is able to form sandy mounds on the seabed, both in intertidal as in subtidal areas Rabaut et al. 2009). In turn, these small-scale features interact with large-scale morphology, such as sand waves, but to an extent that is yet unknown.
A wide range of anthropogenic activities, such as the construction of wind farms, sand extraction, shipping and fisheries are increasingly pressuring the coastal environment (Halpern et al. 2008). To ensure a sustainable use of the resources in these areas, knowledge of the impacts of these activities is essential. Sand extraction, for instance, impacts both the biological and physical landscape, whereas only little is known about the effects of such an intervention (De Jong et al. 2015a). One way to study these effects is by means of process-based biogeomorphological models, i.e. models that are based on a mathematical representation of the fundamental laws of physics. However, current state-of-the-art numerical models (Van Gerwen et al. 2018;Campmans et al. 2018)used for finite-amplitude sand wave modellingdo not account for biological processes. Conversely, idealised modelling studies have studied the feedbacks between sand waves and benthos Damveld et al. 2019). Although these model studies showed that including biophysical interactions can significantly change the initial formation properties of sand waves, they did not provide information on finite-amplitude properties of sand waves due to the applied modelling techniques. As a consequence, these models are less suitable for studying the spatiotemporal effects of offshore engineering operations on the biological and morphological landscape.
Hence there is a need for biogeomorphological models that provide a mechanistic understanding of the processes governing the nonlinear feedbacks between biology and  finite-amplitude sand waves. In this respect,  proposed a method to describe mound formation by L. conchilega in a process-based manner. Combined with a readily available nonlinear sand wave model (Van Gerwen et al. 2018), this allows for the exploration of the feedbacks between tube-building worm patches and sand waves. Therefore, in this work we present a two-way coupled biogeomorphological model in which both (i) the influence of patches of L. conchilegaon the morphological evolution of sand waves and (ii) the nonlinear effects of the physical environment on the population distribution of this organism can be studied. Using this model, we aim to determine how the feedbacks between patches of tube-building worms and sand waves affect the morphodynamics of the marine landscape on varying spatial and temporal scales. This paper is organised as follows. The next section provides further background information on the marine underwater landscape and biogeomorphological models. The model formulation is then given, which is followed by the results. Finally, the discussion and conclusions are described.

Background
Biota and bedformsthe bottom of the sea Coastal areas are inhabited by a rich community of benthic species Callaway et al. 2002). Knowledge about these species is of great importance for both ecology and engineering. Benthic communities are, for instance, an indispensable link in the marine food cycle (Petersen and Exo, 1999). Moreover, benthic organisms are able to influence sedimentary processes on a large spatial and temporal scale (e.g. Widdows and Brinsley, 2002, and references therein). Some may even reshape the abiotic environment in such a way that they create and maintain their own habitat and that of others: the so-called ecosystem engineers (Jones et al. 1994;Zühlke, 2001;Meadows et al. 2012).
The species of interest in this paper (the polychaete L. conchilega; see Figure 2) is such an ecosystem engineer. It lives in dense patches in sandy environments. The patch densities show great variation, but are often found to be in hundreds up to thousands of individuals per metre squared (Rabaut et al. 2007;Van Hoey et al. 2006;Degraer et al. 2008;Van Hoey et al. 2008;De Jong et al. 2015b). The tubes of this worm protrude a few centimetres into the water column, affecting the thin boundary layer above the bed (Friedrichs et al. 2000). As such, sedimentation and erosion patterns around the patches change with respect to an abiotic environment, leading to the formation of sandy mounds of several tens of centimetres in height (Rabaut et al. 2009;. Rabaut et al. (2009) defined the biogenic structures created by L. conchilega as reefs, and suggested that the communities could persist for years. Other studies corroborate on this and show that the lifespan of the worms is indeed several years (Beukema et al. 1978;Ropert and Dauvin, 2000). Conversely, since the tube-building worms are sensitive to low water temperatures, high mortality rates within the communities are observed during severe winters (Strasser and Pieloth, 2001;Alves et al. 2017). These long-lasting sandy mounds also attract other species through which they alter the local benthic communities, making them very valuable ecosystem services (Zühlke, 2001;Rabaut et al. 2010).
One of the most important indicators of a healthy functioning ecosystem is its biodiversity (Naeem et al. 2012). The biogenic structures formed by L. conchilega patches provide shelter for other species through which they stimulate the local biodiversity. The protection of these structures is therefore crucial in sustainable marine management (Braeckman et al. 2014). Bottom trawling is shown to have a destructive effect on these reefs . Moreover, dredging plumes potentially lead to higher mortality rates (Mestdagh et al. 2018). Hence information on the spatial distribution of these communities is essential in designating areas protected from anthropogenic habitat modification, such as dredging and demersal fishing (Sheehan et al. 2015).
It is well known that certain abiotic variables (e.g. temperature, sediment characteristics and chlorophyll-α content) drive the composition of benthic communities Künitzer et al. 1992;Reiss et al. 2010). Also, for tidal sand waves a clear relation has been established between their morphological structure and the occurrence of benthic organisms. Baptist et al. (2006) investigated seasonality in spatial variations and found evidence that endobenthic organisms occur in higher densities in sand wave troughs. Furthermore, Damveld et al. (2018) reported significantly higher population densities of both epi-and endobenthic organisms in sand wave troughs, compared to the crests.

Biogeomorphological modelling
Organisms are often categorised in functional groups by their effect on morphological processes (i.e. stabilisers and destabilisers Widdows and Brinsley, 2002). The inclusion of these biological processes has made significant improvements in geomorphological modelling (Corenblit et al. 2011). Originally, biogeomorphological research was only done in a one-way fashion. For instance,  investigated the influence of biological processes on the formation of sand waves in an idealised process-based model. Using the complex process-based model Delft3D (Lesser et al. 2004),  modelled the formation and degradation of sandy mounds by L. conchilega. On the other hand, Willems et al. (2008) and Cozzoli et al. (2014) used statistical models to predict habitat composition based on environmental characteristics. However, instead of such one-way interactions, including two-way coupled feedbacks between biology and geomorphology could even further improve our understanding of the shaping of marine landscapes (Corenblit et al. 2011).
Recently, Damveld et al. (2019) modelled the two-way coupling between organisms and tidal sand waves. In particular, they coupled a morphodynamic model to a general logistic growth model in which they investigated the organisms' spatial and temporal distribution over sand waves. Their mathematical solution method was based on a linear stability analysis (LSA), a common technique used for sand wave modelling (e.g. Hulscher, 1996;Besio et al. 2006;Campmans et al. 2017). A major drawback of LSA is that it only provides information on the initial formation stage of sand waves (preferred wavelength, growth/migration rate) and not on the height and shape. The nature of this technique is also not ideal for reproducing the patchy spatial distribution of L. conchilega that is observed in the field. To obtain these kinds of results, a nonlinear process-based model such as Delft3D can be used. Van Gerwen et al. (2018) already demonstrated that it is possible to model the equilibrium height of sand waves with this modelling software. Moreover,  showed that the grid-based structure of Delft3D is well suited for the representation of tube-building worm patches. Hence combining these methodologies enables us to model the two-way coupling between tidal sand waves and the biogenic activities of L. conchilega. 2574 DAMVELD ET AL.

Model formulation
In this study, the evolution of sand waves and tube-building worms is modelled using the numerical process-based model Delft3D (Lesser et al. 2004). The model concerning the evolution of sand waves is comprehensively discussed in Borsje et al. (2013 and van Gerwen et al. (2018). In the following subsections a summary of the Delft3D sand wave model is given, followed by a description of the evolution of benthic organisms and the two-way coupling of the biogeomorphological system. Finally, the model set-up is presented.

Hydrodynamics
The hydrodynamics are described by the shallow-water equations a continuity equation, supplemented with appropriate boundary conditions. Turbulence is represented by a spatiotemporally varying vertical eddy viscosity. The model equations are solved by applying sigma layering in the vertical. Following a 2DV approach we consider variations in the horizontal x and vertical z direction only, thus ignoring typical 3D processes such as the Coriolis effect, wind forcing and the interaction with sand banks. In terms of σ-coordinates, the model equations read where u is the horizontal flow velocity in the x-direction, ω is the vertical velocity in the moving σ-coordinate system, h is the total water depth, ζ is the free surface elevation, ρ w is the water density, P u is the hydrostatic pressure gradient, F u is the horizontal momentum exchange due to turbulent fluctuations, and ν T is the vertical eddy viscosity. The vertical eddy viscosity ν T is calculated using the kÀϵ turbulence model (e.g. Burchard et al. 2008).
The boundary conditions at the bed σ ¼ À1 ð Þand at the free respectively. Here, τ b is the bed shear stress, g is gravitational acceleration, u b is the horizontal velocity in the first layer above the bed and C is the Chézy coefficient, related to the roughness of the bed.

Morphodynamics
Both bedload and suspended load sediment transport are included in the model. Bedload sediment transport S b [kgs À1 m À1 ] is calculated according to in which α s is a slope correction parameter (based on Bagnold, 1966), calculated by with α bs being a user-defined tuning parameter, which is set to 3 following van Gerwen et al. (2018), and Φ the internal angle of friction of sand, set to its default value (30°). Furthermore, ρ s is the specific density of the sediment, d the sediment grain size, u * the effective bed shear velocity, D * the dimensionless particle diameter and T the dimensionless bed shear stress. The suspended sediment concentration in the water column is calculated by solving an advection-diffusion equation: in which c is the sediment concentration, w is the vertical velocity (now in a Cartesian coordinate system, see below) and w s is the sediment settling velocity. Furthermore, ϵ s,x and ϵ s,z are the sediment diffusivity coefficients in the x and z direction, respectively. To distinguish bedload from suspended load, a reference height a=0.01h is defined, which marks the transition between the two modes of transport (Van Rijn, 2007). Below this reference height the concentration is assumed to equal a reference concentration. Above the reference height, transport of suspended sediment S s [kgs À1 m À1 ] is calculated by For further details on initial and boundary conditions (e.g. pickup and deposition), and the relation between the vertical velocities w (Equation (7)) and ω (Equation (1)), we refer the reader to Deltares (2012).
Finally, the evolution of the bed elevation z b is governed by the sediment continuity equation (Exner equation), which reads in which ϵ p =0.4 is the pore volume fraction in the bed. This equation states that a local convergence (divergence) of sediment transport leads to a rising (falling) bed profile.

Biodynamics Evolution equation
After Damveld et al. (2019), we employ the following general equation for logistic growth to describe the spatial and temporal evolution of benthic organisms, while ignoring biological dispersion: Here, φ [ind m À2 ] is the tube-building worm population density (in contrast to Damveld et al., 2019, where φ denoted the biomass density). Furthermore, α g is the logistic growth parameter and φ eq the carrying capacity, which is a function of the bed shear stress τ b (see below). 2575 The solution to Equation (10) can be found analytically and reads where φ init denotes the initial patch density (at t=0). To ensure a smooth introduction of patches in the morphodynamic model, we have set this value to φ init =10 indm À2 .
Carrying capacity The carrying capacity φ eq describes the maximum population density an area can sustain without degradation, given the environmental attributes of that area. A suitable way of describing this aspect is by letting the carrying capacity be a function of the bed shear stress (Damveld et al. 2019), as bed shear stress has been shown to be a good predictor for the distribution of macrobenthos (Herman et al. 2001;De Jong et al. 2015b). Therefore, the carrying capacity is written as where φ eq,flat and τ b,flat are the flat bed carrying capacity and bed shear stress, respectively. Furthermore, · h i tide denotes tide-averaging and κ eq is a parameter controlling the effect of shear stress variations.
We have set the reference value for the flat bed carrying capacity to φ eq,flat =200 indm À2 , which is of the order of the densities observed in subtidal areas De Jong et al. 2015a). Furthermore, the logistic growth rate has been estimated to be α g =0.01 m 2 ind À1 d À1 , such that a newly seeded patch, without being subjected to mortality, reaches its carrying capacity density in around 2 years. The value of κ eq is somewhat arbitrary and has been set to 5 ms 2 kg À1 . A positive value of this parameter indicates that a higher bed shear stress leads to a lower carrying capacity.
Delft3D rigid vegetation module  combined flume experiments and field observations to validate a model tool in which the interaction among the hydrodynamics, sediment dynamics and tube-building worms was quantified. This model tool is the so-called 'rigid 3D vegetation model' of Delft3D, which explicitly accounts for the influence of cylindrical structures on flow and turbulence (Temmerman et al. 2005;Dijkstra and Uittenbogaard, 2010). The influence of the tubes upon the flow is represented by an extra source term of friction force F bio in the momentum equation (Equation (1)), and is given by Here, C d is the cylindrical resistance coefficient and d tube is the tube diameter. Furthermore, as F bio is a function of height, also the tube length protruding from the sediment L tube needs to be defined in the model. In the following simulations, the resistant coefficient is set to 1.0 (default value), the tube protruding length is set to 2 cm and the diameter to 0.5 cm, corresponding to the values used by . The spatial and temporal variation in tube-building worm density φ is based on the variation in bed shear stress, as described above.
The influence of the tubes on turbulence is reflected in an extra source term of total kinetic energy T bio and turbulent energy dissipation T bio τ À1 in the kÀϵ turbulence closure model, where the term T bio ¼ F bio u ð Þrepresents the power spent by the fluid working against the tube drag. Furthermore, τ is a timescale, which reads where τ free is the dissipation timescale of free turbulence and τ tube is the dissipation timescale of eddies in between the tubes. These eddies have a typical size that is limited by the smallest distance between the tubes: Here, C l is a user-defined tuning coefficient, which is set to its default value of 0.8, and A tube is the horizontal cross-sectional tube area. For further information about this module, we refer the reader to the Delft3D user manual (Deltares, 2012).
Initial seeding and mortality In order to be actually able to model the evolution of L. conchilega, we need to specify the initial seeding behaviour, as well as the lifespan of the organisms. Following field observations along the Belgian coast, recruitment occurs in spring, summer and fall (Van Hoey et al. 2006). During these seasons we let the patches evolve according to the logistic growth model in Equation (10). Conversely, the winter season is characterised by a high mortality rate (Alves et al. 2017). Therefore, for a single patch we assume a yearly mortality rate of M 1y =50% with a maximum lifespan of 5 years (M 5y =100%). This dynamic growth process is illustrated in Figure 3, where the full life cycle of a single patch is presented, given a constant carrying capacity.
Field observations by Rabaut et al. (2009) showed that in areas where tube-building worms are abundantly present patches can cover up to 10% of the total area. Given the maximum lifespan of a patch (set to 5 years), we assume that each year in spring 2%(=S y ) of the model domain is randomly seeded with new patches. As a result, the average domain area covered with patches equals 10%. Finally, we set the horizontal length of the seeded patches to L patch =2 m, which again is in agreement with observed patch dimensions in the field (Rabaut et al. 2009).

Model set-up
The numerical model set-up generally uses the same settings as Borsje et al. (2013) and van Gerwen et al. (2018). The horizontal model domain has a length of around 50 km, with a variable horizontal distribution. The central part of the domain is 1 km long with a uniform grid spacing of 2 m, which gradually increases to 1500 m at the lateral boundaries. Vertically, the grid comprises 60 σ-layers, with a resolution of 0.05% of the local water depth near the bed, gradually decreasing towards the free surface boundary. At the lateral boundaries, so-called Riemann boundary conditions are imposed. This type of boundary allows outgoing waves to cross the lateral boundary without reflecting back into the domain (Verboom and Slob, 1984). At these boundaries, the system is forced by the S 2 tidal component with an angular frequency of 1.45×10 À4 s À1 and a depth-averaged velocity of U=0.65 ms À1 . Following van Gerwen et al. (2018), the Chézy roughness coefficient has been set to C=75 m 1/2 s À1 . The hydrodynamic time step of the model is 12 s, and the the computational spin-up time is set to one tidal cycle during which no bed-level changes are allowed. 2576 DAMVELD ET AL.
The initial bed perturbations are situated in the fine-spaced central part of the domain at an average depth of H 0 =25 m, which is representative of a North Sea sand wave location. Three bed configurations are used in this paper: a flat bed, a sinusoidal sand wave pattern with an amplitude of A 0 =0.25 m and a random bed pattern. For the sinusoidal bed, the initial wavelength has been set to λ 0 =204 m. According to van Gerwen et al. (2018), this corresponds to the fastest growing mode (FGM, i.e. the preferred wavelength) for the model parameters used in this work. Because of the stochastic nature of the model we have run each (biotic) case 24 times, where for each simulation we have varied the initial bathymetry and patch seeding in a randomised manner. Exceptions are the abiotic runs, which are run for one (sinusoidal initial bed, so fully deterministic) and 12 times (random initial bed).
A morphological acceleration factor (MORFAC) is used to speed up the geomorphological changes. Van Gerwen et al. (2018) found that a MORFAC of 2000 is still sufficient to obtain accurate results. However, since we are interested in the seasonal growth of benthic organisms, the MORFAC is limited by the biological timescale. Therefore, we use a MORFAC of 182.5, such that one tidal cycle equals exactly one season (13 weeks). Using this MORFAC, one model run of 100 years takes approximately 7 days on a high-performance computing facility. Finally, we rerun the model for each single season using the restart files generated from the preceding season. This allows us to calculate the biological evolution offline, i.e. separately from Delft3D.
For an overview of the above parameters and assumptions, we refer the reader to Table 1.

Flat bed
To understand the influence of patches of L. conchilega on its direct environment, we first present the results of a flat bed . Growth occurs during the first three seasons of the year; the last season is characterised by mortality (indicated by the diamonds). The dashed lines illustrate the logistic growth profile without mortality, starting every winter season. Parameters according to À2 † α g is estimated based on the assumption that a patch evolves to its carrying capacity in roughly 2 years. ‡ These parameters reflect qualitative estimations based on field observations. 2577 with a single patch of tube-building worms. The patch is located in the middle of the domain with a fixed density of φ=100 indm À2 . Figure 4a shows the initial response of the tide-averaged absolute value of the bed shear stress, as well as the response after 5 years. Initially, the bed shear stress shows a small but sharp decrease at the location of the patch. In the lateral direction the bed shear stress reaches its flat bed reference value after several tens of metres. After Other parameters are as presented in Table 1 5 years, the combined presence of the patch and the mound leads to variations in bed shear stress on an even larger spatial scale. These findings are in line with . The mound height after 5 years is plotted in Figure 4 for four different fixed patch densities. It can be seen that for an increasing patch density the mound height increases. For a densely populated patch (φ=500 indm À2 ), the mound height can be up to 50 cm. The trough-to-trough distance is around 70 m for all four patch densities. Furthermore, in the case of simulations with multiple consecutive patches in close proximity of one another (not shown here), one single mound forms eventually where this distance is even longer.

Sinusoidal initial bed
Next, for the initially sinusoidal bed topography the twoway coupled model was run for 100 years (400 seasons). Figure 5 shows three snapshots (of a single run) for the self-organisation of both sand waves and population density. To better visualise the density distribution, the maximum occurring density in the different decades is plotted. During the first decade (Figure 5a), the patches are randomly distributed over the domain and the maximum density is comparable between the patches. Although no distinctive sand wave pattern is visible yet, the local response to the patches clearly stands out, as several mounds are formed. However, this is not yet reflected in the coefficient of determination (R 2 ), which denotes the relation between bed level and patch density. Next, Figure 5b shows that, while the morphological pattern evolves over time, the distribution of the organisms is altered. Higher densities occur in the troughs, rather than on the crests, which is now also reflected in the coefficient of determination. This is even more pronounced after 100 years (Figure 5c), when almost all patches are situated in the troughs (R 2 =0.85). Furthermore, here the (maximum) population densities are higher compared to densities on the initial bed.
The temporal evolution of the bed profile is visualised in Figures 6a, 6b and 6c, where the biotic (green, including biology) as well as the abiotic (brown, excluding biology) evolution is shown. The light-green shade in the figure denotes the standard deviation of all simulations (grey lines), whereas the dark-green line denotes the average. Initially, for the biotic case the crests are slightly higher and the troughs are lower than for the abiotic case. This is the result of the presence of the mounds, as their dimensions dominate those of the sand waves. Furthermore, when evolving towards the equilibrium profile, it is visible (Figures 6b and 6c) that the average biotic sand wave height is a few decimetres lower than in the abiotic case. Figure 6d shows the spacing regularity 1 (i.e. wavelength distribution) of the sinusoidal bed case. The effect of the mounds clearly stands out, as the wavelength distribution is large during the first 200 seasons. Moreover, the average wavelength of the biotic runs is initially smaller than that of the abiotic run. Later, towards the equilibrium stage, the standard deviation of the biotic simulations decreases and its average wavelength becomes similar to the abiotic wavelength, which is based on the FGM, as explained earlier.
For larger values of the flat bed carrying capacity (φ eq,flat = 400 indm À2 ), it can be seen that the same observations hold as for the previous case, albeit more pronounced (Figure 7). In general, the standard deviation of the crest and trough development is larger than in the previous case. Also, the average equilibrium sand wave height is again lower (∼50 cm) when taking into account the effects of the tube-building worm patches. Finally, the difference in average wavelength between the biotic runs and the abiotic run is larger than the previous case, while this occurs over an even longer time period here, as well.

Random initial bed
This case considers the evolution of an initially random bed profile over 150 years (600 seasons). Here, the morphodynamic and population evolution are presented in Figure 8. Again, roughly each 50 years a snapshot is presented for a single run. Note that because of the different initial bathymetry, this case has a generally slower time-toequilibrium compared to the sinusoidal bed case. Similar to the sinusoidal bed case, initially the randomly distributed patches lead to small mounds, whereas patch densities are even throughout the domain. After 50 years ( Figure 8b) the mounds become higher and they evolve in morphological patterns on a spatial scale larger than that of the patches themselves. These morphological features eventually evolve into distinct sand wave patterns (Figures 8c and 8d). Here also the effect on the population distribution is clearly visible, with 1 The spacing/wavelength is defined as the crest-to-crest distance. The crests that are considered here are only those that are higher than 50% of the average wave height of all peaks in the domain at that point in time. The wave height is defined as the vertical difference between the crest and the average of the two adjacent troughs.  higher densities in the troughs compared to the crests. Similar to the sinusoidal bed case, the coefficient of determination (R 2 ) clearly shows the increasingly strong relation of bed level and patch density.
To further analyse the bed development of this case, we present a timestack of the evolution of the random bed (Figure 9). For the biotic case (Figure 9a), the results show that a large number of mounds develop randomly throughout the domain. Remarkably, it can be seen that the mounds trigger the growth of sand waves with a smaller wavelength than that in the abiotic case. The total number of sand waves in the domain around t=400 seasons (emphasised by the dashed white line) is seven for the biotic case, and five for the abiotic case. Moreover, the biotic sand waves feature a faster growth than the abiotic ones. Over time, when the sand waves evolve towards their equilibrium in the biotic run, several crests merge and, hence, the number of sand waves decreases again (see also Figures 8c and 9d). Eventually, the equilibrium wavelengths are in agreement with those of the abiotic run (Figure 9b). Finally, the eventual wavelengths in both simulations are in good agreement with the FGM as used in the sinusoidal case ( Figure 5).
Since the above results only represent a single run, we present in Figures 10 and 11 the crest and trough development, and the spacing regularity of all simulations in the random bed case, respectively. Figure 10 now clearly shows that the growth rates of the biotic runs are larger than the abiotic runs. Also, initially the spread of the biotic crest and trough development is higher, as expressed through the standard deviations. However, towards equilibrium, both the crest and trough levels, and the standard deviations of both cases, become comparable. The average bed levels at the equilibrium stage fall in the same range as observed in the sinusoidal case.
Also, the spacing regularity shows clear differences between the runs with and without tube-building worm patches. Figure 11 shows that the average wavelength of the biotic runs is significantly smaller than that of the biotic runs, in particular in the initial growth stage. Albeit less distinct, also during the secondary growth stage this difference is still present. While evolving towards equilibrium, this effect vanishes and both average wavelengths become similar again. Finally, due to the self-organisational properties of the random bed, towards equilibrium sand waves of different wave heights and wavelengths can be observed, as already shown by Campmans et al. (2018). This is reflected in the increased standard deviation for the abiotic as well as the biotic runs (Figures 10 and 11), compared to the sinusoidal bed case, where the wavelength distribution decreased significantly over time.
Notably, the average wavelengths tend to that of the FGM (∼204 m) only in the final stages over the morphological development, suggesting that the system has not yet reached a full equilibrium after 600 seasons. Moreover, one of the properties of this dynamic system is the merging of crests towards the end of the simulations (as shown in Figure 9). The result is a sudden increase in bed levels and wavelengths, which is visible in some of the single runs in Figures 10 and 11. The spike in crest level is of a temporary nature, as the crest level quickly stabilises again in these situations.

Discussion
The spatial distribution of L. conchilega In this study we developed a process-based model, able to describe the two-way coupling between the morphology of finite-amplitude sand waves and tube-building worm patches of L. conchilega. The model results related the sand wave troughs to the areas most densely populated with tube-building worms. These results qualitatively agree with field observations on general community distribution patterns, as Baptist et al. (2006) and Damveld et al. (2018) showed that benthos prefer to settle in sand wave troughs rather than on the crests. However, since these studies focused on community assemblages in general, it is questionable to what extent these observations may be translated to the distribution pattern of a single species. Yet, for bed forms of larger scale, such as shoreface-connected ridges (Markert et al. 2015) and sand banks (Van der Reijden et al. 2019), it has been shown that tube-building worms are significantly more abundant in the troughs than on the crests, although it should be noted that the study by van der Reijden et al. (2019) focused on the tube-building worm Sabellaria spinulosa. Furthermore, preliminary results from a recent field campaign suggest that also for sand waves L. conchilega is more likely to be found in the trough region than on the crests (Cheng et al. 2020). Still, insight into the distribution of benthic organisms in subtidal areas, particularly in relation to marine bed forms, is scarce. Therefore, we encourage future monitoring campaigns to focus on the population distribution of benthic species on the spatial scale of sand waves.
The effects on sand wave morphology The fastest growing mode (FGM) is assumed to represent the wavelength as observed in the field (Dodd et al. 2003). To bypass the initial growth stageand shorten calculation time we used the FGM as calculated by van Gerwen et al. (2018) for the initial sinusoidal bathymetry. However, as pointed out by Damveld et al. (2019), the inclusion of biological processes may affect properties of the FGM. Using a linear stability analysis, they demonstrated that benthic activity can result in a change in the preferred wavelength over time. Also, the growth rate of the bedforms can be significantly affected by the presence of biota. Indeed, our results led to similar findings.
We showed that the wavelengths of the emerging bedforms were significantly smaller than that of the supposed FGM (up to several tens of metres), which was particularly true for the random initial bed topography. Moreover, the growth rate of these smaller bedforms was higher than that of the FGM, resulting in a much shorter time-to-equilibrium. Here, we illustrate that it is in fact the spatiotemporal interaction between small-and large-scale morphology that leads to these results. Revisiting Figure 4b, we saw that mound evolution occurs on a much shorter timescale (years) than sand wave evolution (decades). However, their wavelengths act on an almost similar spatial scale. Hence it is the formation of biogenic mounds that initially triggers the growth of sand waves with smaller wavelengths than that of the FGM. Yet, in the biotic cases the abiotic FGM did eventually prevail, implying that these processes mainly affect the formation stage and to a lesser extent the equilibrium stage. The driving mechanism in sand wave formation consists of tide-averaged circulation cells in the water column (Hulscher, 1996). They display near-bed flow velocities directed towards the crests, supporting a tide-averaged sediment transport in that same direction. Conversely, the biogenic mounds are formed due to a local decrease in bed shear stress, which leads to sedimentation in the area around the patch (Figure 4), as already demonstrated by . Moreover, here we have shown that the majority of the benthic population eventually settles in the sand wave troughs. As a consequence, the crest-directed sediment transport is decreased, and so sand wave heights are lowered. However, this effect was only small, with equilibrium sand wave height differences of the order of a few decimetres. Moreover, it was only visible for the sinusoidal initial bed profile, as the standard deviation of the bed level in the random case was larger than the actual wave height decrease. Nonetheless, we demonstrated that including biological processes potentially leads to a decrease in sand wave height.

Assumptions and limitations
In this study we let the bed shear stress influence the biological carrying capacity. This resulted in significantly higher population densities in the sand wave troughs, compared to the crests. Although these results generally agree with field observations on benthic community patterns, other indicators for benthic life are available and might perform better (Reiss et al. 2010). In particular, Willems et al. (2008) showed that sediment type and size are good predictors for the presence of L. conchilega. Also organic matter concentration in the water column might be important, as L. conchilega is a suspension feeder (Buhr, 1976). Besides, in contrast to bed shear stress, these other predictors are easier to measure in the field. However, the majority of these indicators are not incorporated into biogeomorphological models, making bed shear stress the most suitable as of yet. Further development of such models should thus focus on the inclusion of other predictors for the occurrence of benthic organisms.
The combination of different model parts resulted in a rather comprehensive parameter space, which inherently leads to uncertainties in the model results. We were able to provide evidence for most parameters based on field observations and other modelling studies. However, information regarding the value of the parameters is particularly scarce for those concerning the biological evolution (Equation (10)). In fact, according to Strasser and Pieloth (2001), population densities reported in field surveys do not necessarily have to correspond to the actual carrying capacities of these areas. Hence this illustrates the difficulties in translating field observations to model parameters. In order to gain more insight into the robustness of the results, we ran some additional simulations where we varied the logistic growth rate (not shown here). As already reported by Damveld et al. (2019), varying the logistic growth rate significantly influences the initial growth rate of bedforms. Indeed, our simulations revealed that a higher logistic growth rate results in faster mound formation and, as such, a shorter time-to-equilibrium of sand waves as well. Moreover, this slightly increased the dominance of the smaller wavelengths during the initial growth stage, similar to the effect of increasing the flat bed carrying capacity (see Figure 7). Although our simulations showed that some results may change due to varying input parameters, they all showed qualitatively the same outcome. In particular, the time-to-equilibrium turns out to be sensitive, but in all simulations the biotic bed development was significantly faster than in abiotic ones.
Additional uncertainties originate from the chosen numerical set-up, as shown by van Gerwen (2016) and Krabbendam et al. (2019). In particular, the horizontal grid spacing, number of vertical layers and boundary conditions may lead to sand waves of different shapes and sizes. In this work we identified another model aspect that affects the results. In contrast to sand waves, biogenic mounds are three-dimensional patterns. In other words, their length and width are of the same order. Consequently, flow is not only diverted over the mounds and patches, but also around them. Flow velocities are thus likely to be overpredicted due to the adopted 2DV approach, and with that also the growth rate of the mounds. However, the validation study by  showed that the flow deflection around the patch is 10 times smaller than over the patch and, hence, focusing on two-dimensional flow only is justified. Nevertheless, the uncertainties originating from these numerical choices and assumptions indicate that the results presented here should indeed be interpreted in a qualitative fashion.

Practical implications
Previous biogeomorphological sand wave models could only be applied to the initial stage of formation Damveld et al. 2019). The model presented in this paper extends these studies by focusing also on the finite-amplitude stage. Consequently, the model can serve as a tool to assess the biological and morphological response to engineering interventions in a sand wave environment. As the model results indicated a clear difference in wavelength distribution and growth rate between biotic and abiotic situations, including biophysical processes would be highly relevant for, for instance, predicting the regeneration of sand waves after dredging activities. Moreover, sand extraction policies aim at a quick recovery of the biological environment (Rijks et al. 2014). The model can help to gain more insight into species composition in sand wave areas by using the tube-worm distribution as an indicator, and thereby the formulation of sustainable engineering strategies. It should be noted that, although the model is able to represent the general patchy distribution of the mounds, it is not suitable to provide their exact location owing to its statistical components. It should rather be interpreted as an overall population density. Hence we would like to stress that biogeomorphological modelling and field monitoring should always be complementary. Still, future planning of offshore engineering projects and marine management directives could seriously benefit from the predictive capacities of two-way coupled biogeomorphological models such as the one presented in this work.

Conclusions
In this work we studied a two-way coupled model of patches of the tube-building polychaete worm Lanice conchilega and finite-amplitude sand waves, using the process-based numerical model Delft3D. We related spatial variations in sand wave morphology to the population distribution of the organisms. It turns out that sand wave troughs are characterised by a higher population density, compared to the crests. At the same time, the presence of the tube-building worm patches alters sedimentary processes in such a way that small-scale sandy mounds are formed on the seabed. In turn, the model results have shown that, compared to a situation without biota, the combined presence of the worm patches and sandy mounds lead to sand waves with both smaller wavelengths and faster growth rates during the initial growth stages. However, near the equilibrium stage the wavelengths tend towards their abiotic counterparts again. Finally, the equilibrium wave height of sand waves turns out to be a few decimetres lower (∼5%) in a biotic situation. However, this is only visible if the initial bed topography consists of sinusoidal sand waves, since the equilibrium wave height differences in a randomly perturbed bed situation are larger than the actual decrease due to benthic activity.