Benthic species as mud patrol - modelled effects of bioturbators and biofilms on large-scale estuarine mud and morphology.

Sediment-stabilizing and -destabilizing organisms, i.e. microphytobenthos (biofilms) and macrozoobenthos (bioturbators), affect the erodibility of muddy sediments, potentially altering large-scale estuarine morphology. Using a novel eco-morphodynamic model of an idealized estuary, we investigate eco-engineering effects of microphytobenthos and two macrozoobenthic bioturbators. Local mud erodibility is based on species pattern predicted through hydrodynamics, soil mud content, competition and grazing. Mud resuspension and export is enhanced under bioturbation and prevented under biostabilization through respective exposure and protection of the supra- and intertidal. Bioturbation decreases mud thickness and bed elevations, which increases net mud fluxes. Microphytobenthos reduces erosion, leading to a local mud increase of intertidal sediments. In multi-species scenarios, an effective mud-prone bioturbator strongly alters morphology, exceeding that of a more abundant sand-prone moderate species, showing that morphological change depends on species traits as opposed to abundance. Altering their habitat, the effective mud-prone bioturbator facilitates expansion of the sand-prone moderate bioturbator. Grazing and species competition favor species distributions of dominant bioturbators. Consequently, eco-engineering affects habitat conditions while species interactions determine species dominance. Our results show that eco-engineering species determine the mud content of the estuary, which suggests large effects on the morphology of estuaries with aggravating habitat degradation.

INTRODUCTION

As transition zones between the river and the sea, estuaries provide valuable ecosystems services, such as natural protection of coastlines, fishing and nursery grounds, and constitute an important habitat for a wide range of organisms (Barbier et al., 2011). Estuarine morphology evolves from a combination of riverine, wave and tidal energy that controls erosion and deposition pattern of sand and mud (Dalrymple & Choi, 2007; Dam et al., 2016; Van der Wegen & Roelvink, 2012; Van Ledden et al., 2004). At the same time, a variety of species thriving in muddy sediments modify their habitat by so‐called ecosystem engineering (eco‐engineering), a mechanism where organisms directly or indirectly change the physical conditions of their habitat (Crooks, 2002; Jones et al., 1994). As a result, eco‐engineering activity leads to changes in the stability of intertidal sediments and can at the same time affect biodiversity through increasing habitat heterogeneity and resulting niche development (Crooks, 2002). These eco‐engineering effects suggest that the distribution of mud flats and the morphology of the estuary may change due to the in situ modifications of physical sediment properties. The active reworking of suspended sediments may also have large‐scale effects on the mud content and the sediment balance of the estuary. These would, in turn, provide feedback on the eco‐engineering species. This general hypothesis will be tested in this paper by numerical modelling.

Species abundance depends on both abiotic conditions and biotic interactions. Species occurrence in multi‐species environments has been previously related to abiotic factors, such as temperature, sediment properties, hydrodynamic stresses and salinity, allowing for species predictions based on statistical relationships found through field measurements (Cozzoli et al., 2017; Herman et al., 2001; Singer et al., 2016; Thrush et al., 2003; Ysebaert et al., 2002; Ysebaert et al., 2003). In general, Ysebaert et al. (2002) and Fujii (2007) reported a reduction of bioturbator abundance from the mouth towards the upstream parts of the Western Scheldt and Humber estuaries. Similarly, along the inundation gradient abiotic and biotic parameters lead to distinct species zonation (Ysebaert et al., 2002). However, the spatial pattern and density of species is also a function of species interactions, such as competition for space and resources, and predation. Foraging pattern as well as the species' individual life‐cycles control reproduction, migration and feeding activity (Wilson & Parker, 1996). Hence, disturbances of the ecosystem that cause habitat alteration or degradation will change the community structure, evoking direct effects on the morphology of the estuary (Heck et al., 2008). We therefore require a more holistic understanding of the interacting biotic and abiotic processes that determine species distribution and abundance, especially in view of increasing pressures by climatic changes and human impacts, species invasion or biodiversity loss.

Mud, a mixture of silt and clay (m), can have strong effects on the morphology of estuaries through stabilization of shorelines and tidal bars (Braat et al., 2017; De Jorge & Van Beusekom, 1995; Mitchener & Torfs, 1996; van de Lageweg et al., 2018). The presence of mud in the bed shapes the morphology of tidal systems through the counteracting effects of sediment refinement and cohesiveness: with increasing mud content, the bed becomes first more erosive because of sediment fining; however, when the bed comprises large mud fractions, the bed becomes cohesive and erodibility is strongly reduced (Le Hir et al., 2007; Van Ledden et al., 2004). As mud settling is largely governed by the hydrodynamics, the largest mud fractions can be found in sheltered intertidal areas leading to higher bed elevations and infilling of small tidal channels (Braat et al., 2017; Brückner et al., 2020; Kleinhans et al., 2018; van de Lageweg et al., 2018). This local stabilization of the estuarine morphology leads to less dynamic channel networks with stable bars and steeper bank slopes that can be colonized by marine organisms and vegetation (Braat et al., 2017).

However, organisms that live within or on top of intertidal mud alter the stability of the mud cover. Microphytobenthos growth enhances seasonal mud stability of intertidal mudflats through secretion of extracellular polymeric substances (EPS). These EPS form biofilms which protect the sediment matrix of the top bed from erosion during spring and summer and can be washed away during storms and floods in autumn (Chen et al., 2017; Daggers et al., 2018; Le Hir et al., 2007; Paterson, 1994; Van der Wal et al., 2008; Vos et al., 1988; Widdows & Brinsley, 2002; Yallop et al., 1994). Their stabilizing effect strongly reduces mud resuspension on intertidal flats and can lead to a local increase in mud layer thickness (Widdows & Brinsley, 2002). As a result, the presence of microphytobenthos alters the resuspension threshold of mud, the erosion by waves and currents and hence the stability of mudflats.

Several studies on intertidal flats reported a shift between stabilization of sediments in spring through encroachment of microphytobenthos and destabilization by grazing macrobenthos during summer and autumn, leading to seasonally altered mud thickness and suspended sediment load on intertidal flats (Herman et al., 2001; Widdows et al., 2004; van de Koppel et al., 2001). Deposit‐feeding macrozoobenthic species can reduce the stabilizing cover of the microphytobenthos through grazing (Widdows et al., 2000). In addition, bioturbation induced by the movement of macrozoobenthos can directly destruct the cohesive mud cover through the creation of burrows and tracks (De Deckere et al., 2001; Montserrat et al., 2008). This effect was measured to be greatest on muddy sediments where bioturbation reduces sediment cohesion (Li et al., 2017). Consequently, macrozoobenthic bioturbation not only decreases the erosion threshold of the sediment but also increases resuspension rates. Systems dominated by biostabilizing or bioturbating species are therefore affected by contrasting resuspension rates of mud that lead to varying mud coverages and protection of intertidal sediments.

Macrobenthic bioturbation is thought to be proportional to the species‐specific metabolic rate, defined by the rate of biological processing of energy and material, that scales with animal biomass (Cozzoli et al., 2017; Cozzoli et al., 2019). Several empirical studies quantify species effects on erodibility, either in flume experiments where the effect of species abundance on resuspension was measured or in field studies where morphological change and macrobenthos presence were correlated (e.g. Cozzoli et al., 2019; Widdows et al., 2004). Especially, the mud fraction in the top layer of the bed was found to be important for species behaviour. On the one hand, fine sediment is rich in nutrients compared to sand that can be grazed by deposit‐feeding organisms, whereas cohesion can hamper locomotion of species (Li et al., 2017). Furthermore, pollutants such as heavy metals and microplastics mainly accumulate in muddy sediments and reduce species productivity (Kröncke et al., 2013). As a result, the bioturbation efficiency of macrozoobenthic organisms, in terms of volume of sediment that is processed per unit of time, depends on their habitat quality and conditions.

While increasing numbers of fieldwork‐based studies and physical experiments confirmed the effects of various biodestabilizing species on sediment erodibility (Cozzoli et al., 2019; Le Hir et al., 2007), only few studies have tried to disentangle local and large‐scale morphological effects of bioturbation. Existing numerical studies model eco‐engineering effects of key species, such as stabilization by biofilms [e.g](Le Hir et al., 2007; van de Koppel et al., 2001), bioturbation [e.g.](Knaapen et al., 2003; van Prooijen et al., 2011), changes in bed roughness (Borsje et al., 2009; Coco et al., 2006) or sediment mixing in the upper layer of the sediment bed (Paarlberg et al., 2005). Many such models are one‐dimensional or cover small‐ to medium‐scale areas, simulating the effect of a single species or few key species. In recent numerical modelling studies large‐scale effects of marine vegetation led to a stabilization of the estuary shores, which limited lateral erosion of the estuary (Brückner et al., 2020; Kleinhans et al., 2018; Lokhorst et al., 2018). Possibly microphytobenthos induces a similar confinement of the estuary, whereas destabilization by macrobenthic organisms might promote erosion of the estuarine flanks facilitating lateral expansion. We hypothesize that the dominance of one or more species controls mud coverage and morphology, where biostabilizers promote confinement and bioturbators expansion of the estuary.

To understand how large‐scale estuarine morphology evolves under the influence of common bioturbators and biostabilizers, we developed an eco‐morphodynamic model that couples a species model with a hydro‐morphodynamic model in Delft3D. The species model computes species occurrence based on environmental parameters and feeds biomass‐dependent effects on mud erodibility into the hydro‐morphodynamic computations based on empirical relations from the literature. We investigate microphytobenthos and two generic macrobenthic species: a mud‐prone effective bioturbator inspired by Corophium volutator and a sand‐prone moderate bioturbator based on Arenicola marina that are both abundant organisms in estuaries in northwestern Europe.

First, we validate the macrobenthic species predictions on a calibrated hydro‐morphodynamic model of the Western Scheldt Estuary as a reference system by comparison to field data. Then, the eco‐morphodynamic model is applied to an idealized model domain to analyze evolving mud distribution and morphology under bioturbators, biostabilizers, and their combination. To quantify the importance of eco‐engineering effects for species abundances, we compare our results to a control run excluding eco‐engineering effects. Our results will inform future management and conservation strategies by improving the understanding of feedbacks between biodiversity and morphology.

METHODOLOGY

In order to quantify large‐scale morphological effects of microphytobenthos and macrobenthos in estuaries, we applied a novel eco‐morphodynamic model consisting of a dynamic species model coupled to the hydro‐morphodynamic model of Delft3D. The species computations and the hydro‐morphodynamic model were coupled each morphological month.

The species model is literature‐based and comprises two parts: a macrobenthos module and a microphytobenthos module. The former determines species abundance with biomass‐dependent effects on resuspension threshold and rate. The latter computes presence‐absence of microphytobenthos based on habitat conditions and increases the resuspension threshold. The model includes species interactions, such as competition between species and grazing of macrobenthos on the microphytobenthos.

Below we first describe the hydro‐morphodynamic domain before we outline the equations that determine temporal and spatial benthic growth, species interactions and eco‐engineering effects.

The Western Scheldt domain

To compare species predictions of the macrobenthic species with field data, we used a calibrated and optimized hydro‐morphodynamic model of the lower Western Scheldt estuary (Nevla‐model; for more information see Vroom et al. (2015), Schrijvershof & Vroom (2016), and Brückner et al. (2020)). The model was adapted from the generic vegetation simulation in Brückner et al. (2020) (bathymetry of 2008) with enhanced mud deposition on the bars. We computed the species distributions on the initial hydro‐morphodynamic domain that provides values for the environmental parameters that define species growth. This one‐time simulation was used to validate trends in species predictions before the idealized model domain was applied to test our hypotheses.

The idealized hydro‐morphodynamic domain

We run the hydro‐morphodynamic computations on a depth‐averaged (2DH) configuration of Delft3D. Delft3D is a throughout applied and validated software package solving the shallow water equations (Lesser et al., 2004). Sediment transport is computed using the Engelund‐Hansen equation for total transport of sand. The total load predictor of Engelund‐Hansen is chosen since it results in the most realistic bar and channel patterns in long‐term morphodynamic models (Baar et al., 2019). The transverse bed slope effect was enhanced by using an of 0.2 to increase downslope sediment transport and counteract unrealistic grid size‐dependent incision (Baar et al., 2018; Baar et al., 2019). The Partheniades‐Krone formulation (Partheniades, 1965) computes the sediment flux of mud (kg ) as with M (kg ) the erosion rate, () the critical bed shear stress, maximum bed shear stress (N ) induced by the hydrodynamics and S erosional step function that is when the bed erodes and for static or sedimentary beds.

Delft3D solves morphology through an advection‐diffusion scheme with sand‐mud mixture that includes a critical volumetric mud fraction in the bed set to . For mud fractions above this threshold, erosion of both mud and sand is predicted by the Partheniades‐Krone formula (Van Ledden et al., 2004).

The idealized estuary was inspired by a trumpet‐shaped estuary (Dalrymple & Choi, 2007) and the model developed in Braat et al. (2017). Tides and geometry are based on the Western Scheldt Estuary, a meso‐ to macrotidal, tide‐dominated estuary. The initial estuarine bathymetry evolved from a constant river discharge and M2‐tide (lunar semi‐diurnal), equilibrium sand transport and mud supply at the river boundary. Through this set‐up, we assume a well‐mixed estuary and neglect salinity effects (Olabarrieta et al., 2018). Waves were excluded from the model for reasons of computational cost as well as to mimic the tide‐dominated Western Scheldt Estuary. The initial condition for the eco‐morphodynamic simulations was defined after the development of a realistic bathymetry with sandy bars and muddy shores (Braat et al., 2017; van de Lageweg et al., 2018). First, an initial sandy domain was obtained by computing 2000 years of morphodynamics. After the development of a realistic sandy morphology, riverine mud was included for a simulation time of 280 morphodynamic years, leading to muddy shores and tidal bars (Figure 1).

FIGURE 1

Initial bathymetry of the domain with river boundary (right) and tidal boundary (left) in yellow‐blue colours. The mud deposits are displayed as contoured mud fractions along the intertidal floodplains and on tidal bars. The grey line denotes the sub‐intertidal boundary [Colour figure can be viewed at wileyonlinelibrary.com]

The ecological and morphological time‐scales were set equal. To upscale the hydrodynamics to morphological time‐scales, a morphological acceleration factor of 60 was applied (Lesser et al., 2004; Ranasinghe et al., 2011). This means that one tidal period represents one morphological month and 12 tidal cycles represent one morphological year. After each tidal cycle, the ecological and the hydro‐morphodynamic models were coupled to update the maps of species abundance and mud erodibility used in the hydro‐morphodynamic computations. Total simulation time was 50 morphological years, hence, a hydrodynamic simulation time of 600 tidal cycles. The specific model parameters are defined in Table 1.

TABLE 1

Model parameters as defined in the Delft3D model and the ecological module

Parameter	value	unit
Numerical settings
Simulation time ecological year	6	days
Numerical time‐step	0.2	min
No. ecological time‐steps	12	per ecological year
Grid cell size	50 by 80 ‐ 125 by 230	[m by m]
Boundary conditions
Tidal amplitude	2	m
Principle tidal period	12	h
River discharge	100	m3/s
Mud input at river	100	mg/l
Sand
Median sand diameter	0.3	mm
Dry bed density	1600	kg/m3
Mud
Settling velocity of mud	2.5e−4	m/s
Crit. bed shear stress for erosion	0.2	N/m2
Erosion parameter	1e−4	kg/m2/s
Dry bed density	1600	kg/m3
Morphology settings
Transverse bed slope parameter Ash	0.2	[‐]
Active layer thickness	5	cm
Max. storage layer thickness	5	cm
Morphological acceleration factor	60	‐

The macrobenthos computations

We parameterized two generic intertidal bioturbators that are based on the two macrobenthic species Arenicola marina (AM) and Corophium volutator (CV), both abundant species in estuaries of NW‐Europe, such as the Western Scheldt Estuary (Cozzoli et al., 2013; Ysebaert et al., 2003). Moreover, the two species differ in their distinct habitat preferences and bioturbation efficiency, specifically sediment properties and inundation gradient (Beukema & Flach, 1995; Ysebaert et al., 2002).

A. marina (lugworm) is a deep‐burrowing polychaete that lives in intertidal sediments with low silt content (Beukema & Flach, 1995). Even though they are large in size, their bioturbation efficiency is limited by the typically low abundances (< 100 ind./: individuals per unit area). It forms J‐shaped burrows with feeding pits and pseudo‐faeces that increase sediment exposure to near‐bed flow (Beukema & De Vlas, 1979; Wendelboe et al., 2013; Volkenborn et al., 2009). C. volutator (mud shrimp) is an intermediate‐burrowing Amphipode that builts U‐shaped burrows and actively irrigates the sediment (De Backer et al., 2011). Additionally, C. volutator acts as a deposit‐feeder that grazes microphytobenthic diatoms (De Backer et al., 2010). With its preference for high mud content in the bed, its habitat is typically located at the higher intertidal and under low hydrodynamic energy conditions. With high densities up to 20,000 ind./ its bioturbation efficiency is very high (De Backer et al., 2011). In habitats where both bioturbators occur, A. marina outcompetes C. volutator, leading to a distinct boundary between A. marina and C. volutator occurrence (Beukema & Flach, 1995; Herman et al., 2001).

The macrobenthos module, detailed below, comprises two parts: (i) species distribution computations based on environmental parameters and (ii) relations between species abundance and bioturbation effects used for calculating the spatially varying erodibility of the mud. Each coupling, a new abundance of the macrobenthic species is computed for each grid cell based on the output from Delft3D. The resulting eco‐engineering effects are fed back into the hydro‐morphodynamic model as a modified spatially varying, species‐ and biomass‐dependent critical bed shear stress and erosion parameter of the mud. The updated erodibility of the mud was fed back into the Delft3D model.

Species distribution model based on environmental parameters

To allow for realistic abundances of various macrobenthic species, the model determines species occurrence based on the relations reported in Ysebaert et al. (2002) from the Western Scheldt Estuary: The relative macrobenthos biomass is computed based on maximum flow velocity , mud content and inundation period as which determines habitat suitability based on the limiting environmental parameter (Figure 2). As a result, is the smallest biomass that can survive based on the three environmental parameters and is relative to the maximum biomass that can occur if the conditions were ideal (stars in Figure 2).

FIGURE 2

The species occurrence as a function of three environmental parameters for the moderate sand‐prone bioturbator AM (left) and the efficient mud‐prone bioturbator CV (middle). The star represents the maximum biomass occurrence for each parameter. The grey area shows the habitat area that the species thrives in. Right: MPB occurrence depends on inundation period and mud fraction in top layer as bimodal distribution (presence‐absence) [Colour figure can be viewed at wileyonlinelibrary.com]

Even though sediment type, hydrodynamic forcing and salinity are the main determinants of habitat suitability (Cozzoli et al., 2017; Ysebaert et al., 2002; Ysebaert et al., 2003), we neglect salinity as we assume the generic species to dominate along the entire estuarine gradient. Moreover, we focus on intertidal bioturbation, which is a fair assumption as mud settling is limited within the channels of deep estuaries due to high flow velocities and sediment instability (Ysebaert et al., 2003; Heip & Herman, 1995). Seasonality is neglected as we focus on the maximum bioturbation effect of the species persistent throughout the year.

Relations between species abundance and bioturbation effects

We identified several mechanisms induced by macrobenthos that affect the erodibility of the sediment, such as bioturbation through feeding and reproductive movement, sediment sorting and faecal pellet formation that are expelled at the sediment surface, and grazing of diatoms (Le Hir et al., 2007). These different mechanisms are species‐dependent behaviours that can be approximated by an alteration of the erosion rate M (kg ) and the critical bed shear stress () of the sediment (Cozzoli et al., 2019; Wrede et al., 2018). As a result, the sediment flux of cohesive sediment computed by the Partheniades‐Krone formulation (eq. 1) is affected by the presence of the organisms. A change in erosion parameter (M) represents the altered erosion flux from the top layer of the bed induced by species motility, sediment sorting or pellet formation. A reduction of the critical bed shear stress () leads to enhanced resuspension by sediment fining, bioturbation or pelletization.

Bioturbation effects depend on species density and size (De Backer et al., 2011). To derive a linear relation between species fraction and increasing erodibility, we used the relations reported in Cozzoli et al. (2019), who measured species‐, size‐ and biomass‐dependent resuspension in controlled flume experiments. First, we determined the critical bed shear stress under the influence of bioturbation for a given species at a low resuspension rate of 25 mg and set this as the critical bed shear stress (Widdows et al., 1998). To derive the erosion rate of the sediment under bioturbation, we then fit a linear function between this value and the maximum resuspension at an intermediate bed shear stress set at 0.5 Pa in the experiments. The slope of this function is the erosion rate M, describing the amount of resuspended sediment per time unit:

This derivation of the erosion rate assumes that the resuspension rate is constant along the measured time in the experiment. However, empirical values of critical shear stress of the sediments used by Cozzoli et al. (2019) differ from those in the model. To scale the values from the sediment used in the experiments to our model sediments, we computed the factor between the derived values for and M and those of the uncolonized sediment. The values of the model sediment (0.2 Pa for and 0.0001 for M) were multiplied with the factor representing the maximum bioturbation effects at a bioturbator fraction of 1 (for more information see supplementary material). For a macrozoobenthos fraction below unity, and as used in the model were calculated as

The microphytobenthos computations

The microphytobenthos (MPB) is modelled as presence‐absence and recalculated every coupling interval based on inundation period and mud fraction in the top bed (Figure 2). Based on Widdows and Brinsley (2002) we prescribed MPB growth in numerical cells that were inundated between 35‐80% and contained mud fractions above 0.3 (Daggers et al., 2018; van de Koppel et al., 2001), indirectly selecting sheltered areas with limited erosional and depositional processes (Herman et al., 2001).

Linked to the secretion of extra‐polymeric substances (EPS), MPB stabilizes the sediment and reduces local erosion (van de Koppel et al., 2001). We account for this effect by an alteration of the critical bed shear stress of the mud fraction. Hereby, we assume that MPB live on top of the sediments and increase the critical bed shear stress for erosion by a factor of four during summer (Le Hir et al., 2007) and averaged over one year.

Consequently, the updated abundances emerging from the hydro‐morphodynamic computations alter local mud erodibility and allow for a direct response of the morphology to species dynamics. Moreover, our modelling species are based on contrasting species both in terms of their habitat and their eco‐engineering effect. All values for critical bed shear stress and erosion parameter can be found in Table 2.

TABLE 2

Physical parameters of the abiotic reference run REF, the two bioturbators AM and CV and microphytobenthos and MPB

Parameters	REF	AM	CV	MPB
τcr, e [N/m2]	0.2	0.11	0.02	0.35
M [kg m−2s−1]	0.0001	0.0005	0.0004	0.0001

Notes. and M represent the value for a maximum biomass under ideal environmental conditions. For lower biomasses this value is linearly adjusted between that maximum and the abiotic value (REF).

Competition between species and grazing pressure

We introduced competition between the macrozoobenthic species, by assuming that A. marina is dominant over C. volutator following Beukema and Flach (1995) and Herman et al. (2001). This assumption results in the disappearance of CV in model cells that are suitable for both AM and CV. This presence‐absence relation is independent of the biomass of each species. Competition is computed at each coupling, meaning that the abundances based on the environmental parameters are instantaneously updated.

Moreover, we considered species‐specific interactions between the macrozoobenthic species and MPB. CV grazes on MPB, which leads to a linear reduction in MPB cover depending on the biomass of CV. The new of MPB is defined as the mean value between , which is the value when only MPB is present, and the bioturbator value of the present biomass as

Moreover, MPB is dominant over AM, setting the critical bed shear stress to . However, we assume that AM bioturbates the lower sediment layers underneath the biofilm and set the erosion parameter of AM.

Model scenarios

The eco‐morphodynamic model was applied to two different set‐ups. First, we simulated the calibrated hydro‐morphodynamic model domain of the lower Western Scheldt Estuary to compare model predictions with field data. Second, we simulated an idealized domain to investigate and quantify the interactions between eco‐engineering effects and estuarine morphology. A variety of scenarios were executed, including reference, single species and multiple species scenarios (Table 3).

TABLE 3

Model scenarios for species in isolation and combination, their eco‐engineering effects and species interactions

Model scenario	Species	Eco‐engineering effects	Interactions
E‐REF	‐	‐	‐
E‐AM	AM	yes	‐
E‐CV	CV	yes	‐
E‐AMCV	AM & CV	yes	competition
E‐AMCVMPB	AM, CV & MPB	yes	competition & grazing
E‐MPB	MPB	yes	‐
N‐AM	AM	no	‐
N‐CV	CV	no	‐
N‐AMCV	AM & CV	no	competition
N‐AMCVMPB	AM & CV & MPB	no	competition & grazing
N‐MPB	MPB	no	‐
N2‐AM	AM	no	‐
N2‐CV	CV	no	‐
N2‐AMCV	AM & CV	no	no
N2‐AMCVMPB	AM & CV & MPB	no	no & no
N2‐MPB	MPB	no	‐

Notes. The dynamic eco‐morphodynamic model (E) represents the scenarios with eco‐engineering, the control model (N) without eco‐engineering effects allows for quantification of the eco‐engineering effects for all species and their combination. To investigate the role of species interactions a second control model excludes both eco‐engineering effects and species interactions (N2).

In order to quantify the importance of eco‐engineering effects for emerging species abundance, we compare the results from the dynamic eco‐morphodynamic model with control runs (N‐scenarios). The controls consist of models with species distributions calculated as a function of the physical conditions, but without the feedback between eco‐engineering effects and the hydro‐morphodynamics.

RESULTS

Validation of species predictions

To test the predicted species abundances by our model, we first compare the results obtained from the Western Scheldt model domain (bathymetry 2008) with data from Ysebaert et al. (2002) that includes sampling data from 1978‐1992 (Figure 3). The model predictions represent the trends observed in the data well. As shown from the data, A. marina occurs along tidal bars towards the main channels, whereas C. volutator is locally constrained to the bars and flanks of the estuary. The model predicts AM on the small and dynamic tidal bars and lower flanks. CV is restricted to the higher elevations in the centre of bars and along the shores. The considerable data sampling provided extensive information on species occurrence over more than a decade but was limited to the sampling locations, which is why the model predicts comparably larger abundances. However, species predictions deliver satisfying results to characterise general species occurrence of the two generic species. As a next step, we apply the model to the idealized domain to investigate species effects on morphology.

FIGURE 3

Occurrence of A. marina (AM) and C. volutator (CV) as presence‐absence in the Western Scheldt Estuary. (a) is data partly provided by Ysebaert et al. (2002); (b) is the results from the Nevla‐model domain for AM and CV (more information on the model in Braat et al. (2019) and Brückner et al. (2020) showing similar trends in species predictions [Colour figure can be viewed at wileyonlinelibrary.com]

Spatial patterns of species abundance and eco‐engineering effects

In the idealized estuary model, each scenario resulted in varying distributions of bioturbators and biostabilizers as a result of the species‐specific habitat requirements and eco‐engineering effects (Figure 4a). While AM was widely abundant in the centre and mouth reaches of the estuary (red colours), CV was restricted to the shores and the bars (purple colours). The combination AMCV led to an expansion of AM towards the upper shores and restricted CV abundance to the centre shorelines of the estuary (red and purple colours Figure 4a). In the scenario combining the three species (AMCVMPB), CV and MPB co‐existed on the bars and along higher shores (orange colours). In the scenarios with only MPB, the biofilms occurred on bars and along the main channel in the centre and mouth reaches (black colours).

FIGURE 4

(a) shows bathymetry at the end of the simulation period with presence‐absence for the bioturbators (red and purple), MPB (black), and cells containing both CV and MPB in AMCVMPB (orange). The grey line (MLW) and black line (MHW) indicate the limits between subtidal‐intertidal‐supratidal domain. Bioturbator occurrence is species‐dependent: AM is highly abundant while CV inhabits the highest intertidal cells along the shores. MPB mostly establish on the bars in the mouth and along the estuary shorelines. In (b) and (c), mud fraction difference in the top layer and bed level difference to a reference run without biota, respectively, show erosion for all bioturbating species (blue colours) both at the fringes and on the bars. The bed level changes (c) lead to dynamic variations between accretion and erosion with parts of the mud being transported to the mouth of the estuary. Arrows define the reaches mouth, centre and upstream of the estuary used in the subsequent figures [Colour figure can be viewed at wileyonlinelibrary.com]

Differences in mud fractions arose between the scenarios with eco‐engineers and the reference without species. Scenarios that included bioturbators contained reduced mud fractions, especially along the estuary fringes in the centre and on the bars in the mouth (Figure 4b). MPB locally enhanced mud fractions along the main channel and on the bars outside the estuary. Compared to AM, scenarios that included CV led to a larger mud reduction, including the erosion of mud fractions upstream. This suggests that CV determined the resuspension of mud compared to AM. However, the stabilizing properties of MPB in the combined AMCVMPB‐scenario lowered erosion along the shorelines in the centre and partly even stabilized the mud on the bars in the mouth.

50 years of morphological evolution led to a significant reduction of the bed elevations along estuarine shores. This effect was most pronounced in scenarios that included CV, implying that CV exerts a major control over the morphology. The eroded sediment was partly redistributed towards the adjacent bars and intertidal areas and partly exported towards the sea (Figure 4c). Hence, CV contributed to an increase of mud fractions and bed elevations offshore the mouth area. In comparison, AM led to limited effects on bed elevation change. In contrast, MPB promoted higher bed elevations on the floodplain and along the shores close to the main channel. Since the biologically mediated rate in local sediment erosion led to varying sediment export rates towards the sea, benthic organisms can be important determinants in the retention and export of sediment in and from estuaries.

In Figure 5 we compare the results from the dynamic model (E) with the control runs that excluded either eco‐engineering effects (N) or both eco‐engineering effects and species interactions (N2). This allowed us to quantify the impact of eco‐engineering activity and species interactions on the species abundances in multi‐species runs. Comparison between the static and dynamic model predictions reveals that the static model results in larger species extent and fractions (Figure 5a). AM shows different species abundances between the static and dynamic scenarios. The inclusion of eco‐engineering effects significantly reduces species cover predicted by the dynamic model for scenarios including CV, especially in the centre estuary. MPB predictions are independent of eco‐engineering effects leading to similar predictions between the dynamic and static models. Eco‐engineering and species interactions benefited the habitat extent of AM and negatively affected CV and MPB in all scenarios (Figure 5b). This led to a strong reduction of habitat area for CV and MPB in the multi‐species scenarios. Hereby, the comparison of the N and N2 models shows that competition and grazing had a lower impact on species abundance than the eco‐engineering effects. In contrast, AM occupies less habitat in Figure 5b, showing that eco‐engineering induces a positive effect on AM. The mean species fraction, describing the density of the species within their habitat, was negatively affected by eco‐engineering (Figure 5c). This was true for both AM and CV, whereas MPB remained at a fraction of unity in the single‐species scenario which was related to the binary formulation in the microphytobenthos model. In the multiple‐species runs, however, MPB fractions were reduced both due to grazing and eco‐engineering. Surprisingly, CV showed largest fractions in the static model that included species interactions (N). Possibly, the reduction in their habitat area led to a concentration of higher fractions in the remaining habitat.

FIGURE 5

(a) The maps show the difference in biota fraction between the dynamic (E) and the static (N, N2) scenarios, with predicted larger (grey colours) or smaller fractions (red colours) by the static scenario. (b) Habitat extent of all species predicted by the three models E, N and N2 (in shades) in all scenarios showing an increase of area for AM through eco‐engineering and competition and the opposite trend for CV and MPB. (c) Mean species fraction of colonized cells for all species in E, N and N2 models (in shades) showing a decrease in mean fraction with eco‐engineering effects [Colour figure can be viewed at wileyonlinelibrary.com]

The control runs predicted larger biota abundance than the dynamic model with eco‐engineering for all scenarios but AM. This was due to the different eco‐engineering effects changing the morphology: CV resuspends mud and reduces its own habitat, whereas AM is able to maintain and expand its habitat through reduction of the local mud fraction. Thus, in multi‐species runs, the expansion of the intertidal area and a redistribution of the mud by the presence of CV resulted in newly created habitat that could be occupied by associated species. However, the species density, here described as species fraction, was negatively affected by the bioturbating activity of both AM and CV. Consequently, the results show that species interactions determine where species can settle, whereas eco‐engineering controls habitat suitability.

We compare estuary width, mean mud thickness and relative colonized area along the estuary and along the inundation gradient to study eco‐engineering effects and related trends in species abundance (Figure 6a‐j). The presence of bioturbators leads to a widening of the estuary while MPB locally reduces width compared with the reference scenario (Figure 6a). The effects of the species increase away from the mouth of the estuary. Especially under the presence of CV in the scenarios CV, AMCV, and AMCVMPB, the estuary width increases, which can be related to the erosion of intermediate intertidal elevations between 0 and 0.4 of inundation period (Figure 6a and b). The above trend is also reflected in the mean mud fraction (Figure 6b and g): Along estuary, the reference scenario showed strong variations in mud thickness from a few centimetres up to 0.5 m in the centre of the estuary. In line with the previous results, AM showed comparably limited variations in mud thickness as opposed to CV, AMCV and AMCVMPB that reduced mud thickness along the entire length of the estuary (Figure 6b). CV and AMCV induced strongest effects on mud thickness in the mouth area of the estuary in contrast to scenario AMCVMPB (dark purple) that mainly altered mud in the centre and upstream reach. These variations were likely induced by the presence of MPB that stabilized the mud along the shores in the centre and upstream (Figure 4b). These locations coincided with the alterations induced in the MPB scenario at kilometre 15. For the single‐species runs (AM, CV, MPB) the along estuarine mud thickness could directly be related to relative colonized area (Figure 6c‐e). For multi‐species scenarios, no clear pattern was visible.

FIGURE 6

Estuary width, mean mud layer thickness and relative species cover along the estuary (a‐e) and along the inundation gradient (f‐j) at year 50 for all scenarios showing a widening coinciding with mud erosion along the entire estuary for CV, AMCV and AMCVMPB and accretion by MPB while AM has marginal effects. When multiple species are present, species abundances of AM and MPB are locally enhanced while CV abundances reduce. The new abundances result from an additional feedback between eco‐engineering effects and species interactions and lead to a redistribution of species along the inundation gradient. Please note the logarithmic y‐axis in panel A and F [Colour figure can be viewed at wileyonlinelibrary.com]

Along the inundation gradient, mud thickness reduced from high to low intertidal areas. CV, AMCV and AMCVMPB significantly reduced mud thickness at higher intertidal elevations at inundation periods below 0.4 (Figure 6g). This trend is proportional with an increasing fraction of AM from single‐ to multi‐species runs, whereas CV is decreasing (Figure 6h‐i). In contrast, MPB affects intermediate and high elevations even though its occurrence is limited to the lower intertidal, suggesting that eco‐engineering mediates cross‐shore transport and leads to protection or exposure of the higher intertidal. There is no significant change in MPB between single‐ and multi‐species scenarios (Figure 6j).

Interestingly, the combination of several species (AMCV and AMCVMPB) led to a shift in biota abundance (Figure 6h‐j). While CV coverage was significantly reduced, AM migrated to higher bed elevations, probably because of reduced mud fractions at these elevations induced by CV. Similarly, the presence of MPB in AMCVMPB enhanced both AM and CV presence through mediation of the mud content. Astonishingly, these feedbacks induced a shift in the species curves with AM being dominant at higher elevations and CV at lower elevations in the multi‐species run (AMCVMPB). Consequently, when multiple species were present, species abundances were affected by individual eco‐engineering effects. This means that the emerging species pattern was defined by the new habitat properties emerging from the feedback‐loop between eco‐engineering effects and species interactions.

The biota‐induced changes observed in Figure 4 led to varying effects on the hypsometry across the three reaches of the estuary (Figure 7a). Large bed level changes occurred under CV, AMCV and AMCVMPB. Especially, the upstream and centre reach in these scenarios showed increased erosion of the higher intertidal elevations, which led to a smoothened hypsometric curve. AM induced a similar trend but to a lesser extent. On the other hand, MPB increased the bed levels in the upstream and centre reach with largest accretion around and above mean water (0 m). In the mouth, the concavity of the hypsometric distribution was enhanced especially by CV and AMCV. Overall, the effects on bed elevation distribution were greatest in the upstream and centre reach compared with the mouth area, possibly linked with increasing sediment availability provided from erosion in adjacent areas.

FIGURE 7

Hypsometric curves per reach for bins of 0.2 m (right axis) as mean values and biota probability function (left axis) for each estuary reach show largest biota‐induced effects in the upstream and centre reach of the estuary. Dotted lines represent high and low water. (b) Mud fraction, thickness, tidal area extent and colonized areas per reach. Bars are relative to the values in the reference run [Colour figure can be viewed at wileyonlinelibrary.com]

The bar plot (Figure 7b) quantifies morphological changes for the three reaches of the estuary between the biota scenarios and the reference scenario. Comparison of the mud thickness and tidal area extent revealed how strongly sediment was resuspended through eco‐engineering along the tidal gradient. In general, the effects of AM and MPB were small with greatest effects in the intertidal and limited effects in the supratidal area. In contrast, scenarios including CV strongly reduced mud in the supra‐ and intertidal domain by close to 50% in the upstream and the centre reach. As a result, supratidal area extent decreased, leading to an increase in intertidal and subtidal area extent. The eroded mud partly settled in the subtidal area raising local mud fractions or was transported downstream. In the mouth, more than 50% of the mud was eroded in the intertidal by CV and AMCV compared to the control run, leading to larger mud thickness in the subtidal for AMCV. Interestingly, MPB in the AMCVMPB scenario protected the intertidal mud from erosion with increasing effects towards the centre and mouth of the estuary. Our findings (Figure 7) show that the efficient bioturbator CV promoted regional effects across the estuary by exposing the supratidal domain and promoting lateral expansion of the estuary. On the other hand, the moderate bioturbator AM and microphytobenthos mostly affected local sediments, having a smaller impact on estuarine evolution.

DISCUSSION

Eco‐engineering and species interactions determine species abundance

The results of the multi‐species scenarios revealed new emerging species patterns compared to the single‐species scenarios. When CV is present, the dominant species AM expands its habitat (Figure 4a, 5 and 6): the eco‐engineering effect of CV promotes an advantageous modification of the morphology for AM, allowing for migration of AM towards higher intertidal elevations. Similar successive behaviour can be observed in marine vegetation with pioneer species colonizing harsh environments and through their habitat modification allowng colonization by successive plant species (Bertness & Pennings, 2002; Townend et al., 2011). This mechanism can lead to the disappearance or migration of pioneer vegetation and alters species composition. On the other hand, competitive effects reduce the abundance of CV along the entire tidal gradient. Consequently, competition effects do not only constrain recessive species growth but affect the feedback loop between eco‐engineering and hydromorphology, which governs species community structure. Surprisingly, when MPB is introduced to the model, both AM and CV are able to expand as MPB locally stabilizes the mud, which enlarges the suitable habitat for both bioturbators (Figure 6). New zonation patterns emerge, suggesting that overall species distribution is the result of four combined factors: species‐specific habitat preferences, habitat modification, competition and grazing. Hereby, eco‐engineering effects appear to be especially important for species abundance and ecosystem dynamics (Figure 5). The feedback loop between species interaction, eco‐engineering effects and environment has wide‐ranging implications for species growth and morphology. Eco‐engineering effects are important for species diversity in ecosystems as they promote a variety of habitat and facilitate niche development (Crooks, 2002). These results imply that species interactions have measurable effects on the large‐scale morphology of tidal environments.

The model of the Western Scheldt Estuary predicts trends in bioturbator abundance observed in the data: The model and field data both show the presence of A. marina on the lower elevations closer to the channels, whereas C. volutator occurs along higher elevations at the shores and on top of the tidal bars (Figure 3). The larger abundances in the model result from omitting the salinity gradient in species predictions that in the field data constricts CV abundances closer to the mouth and AM upstream. Moreover, the differences in data sampling locations and timing explain the reduced abundances in the field data compared to the model results. Despite the divergence between the model and field data regarding the overall abundance of species, the model correctly predicted the species distribution trends as observed in the data.

In the idealized estuary, the model can represent species‐specific occurrence along the inundation gradient as reported in Ysebaert et al. (2002), Cozzoli et al. (2013), and Fang et al. (2019): CV grows in high fractions in hydrodynamically calm locations on top of the bars and at the higher elevations along the shores of the estuary; AM occurs at intermediate intertidal elevations in the downstream estuary (Figure 4 and 6h). This trend is in line with species descriptions in literature that report a medium density along large parts of the cross‐sectional gradient of A. marina and high densities of C. volutator (Beukema & De Vlas, 1979; De Backer et al., 2011; Fang et al., 2019). As a result, the model allows us to explore the isolated and combined effects of the eco‐engineering species on the morphology with our idealized scenarios.

Eco‐engineering affects reach‐scale morphology

The eco‐engineering effect of the species is to alter mud distribution in the estuary that feeds back on species growth: a negative eco‐engineering effect of CV causes potentially restricted occurrences through high resuspension, whereas a positive eco‐engineering effect by AM maintains its habitat by resuspending superfluous mud from the top layer (Figure 4 and 6). This mechanism can also be observed in marine vegetation, where dense vegetation cover leads to higher inundation times or ponding, which eventually restricts vegetation to expand or enhances their mortality (Brückner et al., 2019). The bioturbation scenarios represent combined grazing activity and resuspension by movement and reproduction behaviour (Le Hir et al., 2007), hence, a potential maximum effect induced by the bioturbators. As a result, the bioturbation represents ideal conditions for bioturbator activity, assuming sufficient food, absent seasonal variability and neglecting predators. The model results are informative for estuaries that are subjected to maximum biostabilization and biodestabilization effects under ideal conditions, but also elucidate the trends in estuaries under less ideal conditions for benthic species through competition and grazing.

We observe regional variations between mud erosion and bed level change magnitudes in the reaches of the estuary (Figure 6 and 7). Previous studies have shown that the equilibria between bioturbation and biostabilization determine seasonal mud content and bed elevations on intertidal flats (Herman et al., 2001; Widdows et al., 2004; van de Koppel et al., 2001). The presence of biogenic structures, such as shells or protruding structures from burrows, tracks and mounds, can alter the roughness of the sediment and create local turbulence. These effects are highly density and species dependent, possibly creating local skimming flow that reduces near‐bed roughness and therefore exceed the scope of this paper (see Friedrichs et al., 2000). We observe the strongest effects of eco‐engineering in the upstream and centre estuary, resulting from the combination of local species abundances, eco‐engineering efficiency and hydrodynamic stresses.

In the mouth, mainly AM occurs, which is characterized by low bioturbation efficiency and mainly sandy habitat. Furthermore, relatively higher hydrodynamic stresses through a large tidal prism have a larger effect on the morphology of the mouth (Corenblit et al., 2007). Although Volkenborn et al. (2007) showed that the presence of A. marina can lead to system changes between muddy and sandy flats, suggesting that bioturbators can have wide‐ranging effects on system state, the model predicts strong localized effects. This trend evolves from the low mud availability in its mostly sandy habitat combined with low bioturbation efficiency. Consequently, AM exposes a smaller mud fraction and leads to local morphological changes constrained to the intertidal zone. Overall, the effects of AM are smaller compared to the pronounced erosion induced by CV. However, resuspension of sandy sediments under bioturbation was only indirectly accounted for in the model with erosion of sand‐mud mixture for mud percentages larger than 40% (Van Ledden et al., 2004). With this assumption we exclude the morphological effects of species bioturbating sand in the dynamic parts of the estuary that possibly enhance changes in the morphology of the mouth area.

The largest morphological effects occur in the upstream and centre estuary, where CV occurs in calm, muddy areas. Resulting from the upstream erosion, mud is transported towards the mouth but does not accrete on intertidal bars and shores (Figure 7). Instead, the mobilized mud is transported offshore. The upstream erosion also increases intertidal habitat. However, the newly created habitat is not colonized by CV since the reduced mud content constitutes uninhabitable habitat for CV (Figure 7 and S2). On the other hand, MPB stabilizes the intertidal domain upstream and in the centre (Figure 7). As a result, a reduction in mud transport towards the channels leads to reduced mud export while bed levels can be maintained or locally increased. Consequently, large‐scale morphological response is governed by local effects that depend on both the efficiency of the bioturbator in resuspending mud as well as the available mud in their habitat. Hence, the overall potential of the bioturbator to induce modification of system‐scale morphology is determined by a combined effect of species density, species‐specific bioturbation potential and habitat characteristics.

When multiple species are present, morphology is determined by the efficient bioturbator CV. The morphology of the estuaries evolving from scenario CV and AMCV have similar bar pattern (Figure 4), bed level distribution (Figure 7) and mud content (Figure 6). Even with low relative abundances, as AM is the main competitor but with lower bioturbation effect, the presence of CV determines morphology. Consequently, species dominance seems to be secondary in multi‐species environments, whereas presence of species with strong eco‐engineering efficiency is the decisive factor on how the overall morphology evolves. This is similar to the concept of 'keystone species' that have a disproportionately large effect on environmental modification (Power et al., 1996). We therefore confirm that bioturbation efficiency defines the potential of a species to alter their habitat and, independently of abundance, strongly contributes to the morphological evolution of estuaries.

Eco‐engineering affects large‐scale morphology

The morphology of the estuary is altered by eco‐engineering such that the estuary shifts from an exporting system under bioturbation to a depositional system under biostabilization (Figure 4 and 7). Local bioturbators determine total mud content of the estuary through controlling resuspension of mud deposits that are exported towards the sea. This results in an overall sandier environment under the two contrasting bioturbators as well as their combination. On the other hand, biostabilization retains mud in the estuary by protecting the intertidal mud from erosion (Figure 4b).

Consequently, species‐specific resuspension and stabilization lead to contrasting erosion patterns. These modified sediment erosion rates in turn determine bed elevations, affect bank slopes and therefore control slope failure. If bioturbator‐induced erosion is large, i.e. in scenarios with CV, supratidal erosion increases intertidal floodplain extent, whereas lower floodplains transform into subtidal environment. As a consequence the estuary widens laterally (Figure 6a and f). In contrast, lower resuspension by AM transforms the intertidal into subtidal area, reducing floodplain and bar extent without promoting lateral erosion of the estuary. Biostabilization on the other hand protects the banks and promotes local mud accretion, which raises lower intertidal bed elevations. As a result, the trends in large‐scale morphological evolution are governed by the dominance of either biostabilizers or bioturbators. These findings confirm that potential shifts from bioturbator‐ to biostabilizer‐dominance, for example induced by global warming, can affect estuarine evolution, not only locally but also at the large‐scale (Widdows & Brinsley, 2002).

When biostabilization is governing the system, we observe channel confinement and increasing mud content, whereas the presence of destabilizers leads to channel erosion and lateral expansion (Figure 8). Under bioturbation, overall mud content in the system is reduced. When multiple species are present overall channel erosion prevails, even though biostabilizers locally promote bed accretion. What is more, eco‐engineering effects lead to restructuring of species zonation across the system.

FIGURE 8

Conceptual channel adaptation and mud content in an estuary dominated by biostabilization (a), biodestabilization (b), and a combination of biostabilization and ‐destabilization (c). A) a depositional system evolves that develops confining channels with overall increasing mud content and steeper slopes; (b) and (c) promote an erosional systems with lower mud content and gentle slopes, allowing for lateral channel migration. Local biostabilization in (c) partly counteracts the erosional trend. Adapting morphology and species interactions lead to a restructured species distribution [Colour figure can be viewed at wileyonlinelibrary.com]

A positive feedback‐loop emerged between increasing MPB occurrence on mud and subsequent enhanced mud stabilization. This causes seasonal stabilization of mudflats (Le Hir et al., 2007; Orvain et al., 2012; van de Koppel et al., 2001). The mean stabilising effect by MPB leads to a long‐term effect on mud and morphological stabilization. Previously reported modelling results showed that seasonally growing MPB has limited effects on long term bed elevation change as mud gets washed away in winter (Le Hir et al., 2007). However, in subtropical regions or warmer climate, persisting biostabilization possibly allows for long‐term mud accretion in dynamic estuaries (Day Jr et al., 2013). Furthermore, waves also affect the distribution of benthos. Here our models represent sheltered systems by seaward spits or barrier islands with limited exposure to offshore and wind generated waves. The absence of waves in the model allows for MPB establishment on the bars in the mouth where mud can settle. However, in exposed systems those seaward reaches are affected by waves and possibly limit both mud content and presence of MPB.

Previous modelling studies investigating eco‐engineering effects through vegetation on the large scale observed infilling and stabilization in the centre of the estuary which led to confinement of channels and impeded landward erosion (Kleinhans et al., 2018; Lokhorst et al., 2018). MPB shows a similar trend, whereas bioturbation enhances landward erosion. Possibly, systems that include both vegetation and benthic organisms experience mud accretion by vegetation in the higher intertidal areas while low to intermediate intertidal elevations evolve erosion‐deposition patterns induced by alternating microphytobenthos and bioturbator cover. However, species interactions between vegetation and macrobenthos can have additional effects on local species dominance. For instance, van Wesenbeeck et al. (2008) found that presence of A. marina can inhibit saltmarsh establishment linked with the high erodibility of the sediment while A. marina growth is prevented by roots and compaction within the salt marsh sediments. Similarly, seed predation of pioneer vegetation can inhibit vegetation establishment on intertidal flats (Zhu et al., 2016). This wide range of interactions impairs predictions of species effects in multi‐species environments and requires further research.

Species shifts affect morphological adaptation

Especially in the face of a warming climate, we need a more detailed understanding of the morphological response of estuaries to varying eco‐engineers as we expect future species shifts due to habitat degradation, biodiversity changes or species invasion (Dippner et al., 2010). Climate‐ and human‐induced changes will lead to significant effects on local ecosystems and consequently estuarine morphology. Especially macrobenthic organisms were reported to respond to changing environmental conditions. At the same time, field studies have reported a high sensitivity of macrobenthic organisms to climatic changes and external disturbances (Kristensen et al., 2014; Kröncke et al., 2013; Reise et al., 2007; Singer et al., 2016; Van der Wal et al., 2008). For example, increasing temperatures or sediment input can cause significant changes in species communities (Dippner et al., 2010; Pillay et al., 2008). We show that species presence affects the local and the large‐scale morphological development of estuaries (Figure 7). Hence, our results suggest that the response of benthic eco‐engineering species to climate change can be critical to the evolution of estuaries.

Microplastic pollution in estuaries can be mediated by eco‐engineers

The dynamics of estuarine morphology determine the settling and transport of pollutants, such as heavy metals and microplastics that are usually associated with mud deposits. Especially in recent years, the deposition of microplastics (MP) (diameter several millimeters or smaller) was identified as a key component of soil contamination, posing a threat to marine life and human health (do Sul & Costa, 2014). Estuaries are suspected to be one of the key sinks of microplastic deposition owing to the proximity to the river input and tendency to trap lightweight sediments in the intertidal and supratidal zones (Corcoran, 2015). Recent studies showed that accumulation of MP's can statistically be attributed to the distance from a source as well as the properties of the local sediment, such as mud content and grain size (Ballent et al., 2016; Enders et al., 2019; Maes et al., 2018; Vianello et al., 2013), which would predict that the largest MP deposition may occur close to rivers and in the sheltered intertidal sediments. Relations between MP abundance and mud content give a general indication for MP occurrence independently of the MP characteristics, such as density, material or size (Enders et al., 2019). Associated with eco‐engineering activity, sediment resuspension and stabilization possibly determines MP preservation of sediment‐microplastic mixtures. Flushing time of the estuary can be affected by macrobenthos presence, which acts as an active ‘filter’ that affects pollutant export (Dürr et al., 2011). As a result, MP will be exported when strong bioturbation is present. In contrast, biostabilizers possibly retain more MP in muddy sediments and enhance residence times on the mudflats on the tidal bars. Moreover, the ingestion of MP's by marine organisms can lead to reduced productivity or be lethal to benthic species (Besseling et al., 2013; Hope et al., 2020; Wright et al., 2013). Furthermore, MP can be bio‐accumulated along the trophic chain up to humans. Consequently, MP presence can reduce bioturbation efficiency and could change overall morphological response.

Our modelling in conjunction with the above arguments imply that the presence of eco‐engineers has a potentially significant effect on where microplastics settle and accumulate in the system. Based on the observations in the literature and the reasoning above, we expect that there will be a delayed response of the morphology to the supply of MP in estuaries. At first, bioturbators resuspend the MP from the upstream reaches such that they are distributed in seaward direction and possibly captured in the high intertidal zones. As the MP accumulate, the toxicity may reduce the macrobenthic population, which could lead to more trapping at the cost of reducing macrobenthos abundance. This effect may propagate to higher trophic levels leading to ecosystem deterioration. This response has important implications for identifying main areas of pollution in estuaries that can threat marine life and human health. On the other hand, future research on microplastic occurrence needs to take into account resuspension of benthic species to quantify pollutant retention in estuaries and how they are exported towards the marine realm.

CONCLUSIONS

This research shows that benthic eco‐engineers can have large‐scale effects on the morphology of estuaries. To understand the eco‐engineering effects of two macrobenthic species and microphytobenthos, we built a novel eco‐morphodynamic model that combines a hydro‐morphodynamic model in Delft3D with a dynamic species growth model. We ran model scenarios in an idealized estuary domain with a dominant sand‐prone moderate bioturbator, a recessive mud‐prone efficient bioturbator and microphytobenthos in isolation and in combination.

The model showed that the presence of an efficient bioturbator reduced mud thickness and bed elevations in the upstream and centre reach of the estuary through erosion of the supra‐ and intertidal areas. Consequently, the estuary widened laterally. In contrast, a moderate bioturbator mainly affected the intertidal areas at the mouth of the estuary. We found that this contrast resulted from a combination of the bioturbation efficiency of the species and the mud available within their habitat. In runs with multiple interacting species, the morphology is mainly controlled by the efficient mud‐prone bioturbator, even though their abundance is small due to habitat constraints and species interactions. On the other hand, biostabilization by microphytobenthos confined the estuary by protecting the banks from erosion and, hence, facilitating mud accretion in the estuary. Overall, bioturbation led to mud export towards the sea whereas biostabilization retained mud within the estuary.

Eco‐engineering creates suitable habitat for co‐existing species, whereas species interactions determine their abundance. The mud‐prone bioturbator reduces the mud content in the higher intertidal area, which increases habitat quality for the sand‐prone species. Subsequently, once the sand‐prone species colonizes the modified areas, it outcompetes the mud‐prone bioturbator. As a result, a species shift is induced in the upper intertidal where the mud‐prone bioturbator is replaced by the dominant sand‐prone species. On the other hand, microphytobenthos enhances mud fractions in multi‐species runs at intermediate elevations, which in turn facilitates the expansion of the mud‐prone bioturbator.

Eco‐engineers modify the large‐scale morphology of estuaries, whose abundance is determined by both physical habitat modification and by species interactions. As a result, changes in the species community may have wide‐ranging effects on the response of estuarine morphology in terms of mud export to coastal seas and lateral erosion of estuaries under changing conditions. Climate change, human usage and pollutants, such as microplastics, exert pressures on the ecosystem with still unknown consequences and feedbacks that are imperative to future investigation.

Supporting Information

Click here for additional data file.