Effect of hydro-climate variation on biofilm dynamics and its impact in intertidal environments

Shallow tidal environments are very productive ecosystems but are sensitive to environmental changes and sea level rise. Bio-morphodynamic control of these environments is therefore a crucial consideration; however, the effect of small-scale biological activity on large-scale cohesive sediment dynamics like tidal basins and estuaries is still largely unquantified. This study advances our understanding by assessing the influence of biotic and abiotic factors on biologically cohesive sediment transport and morphology. An idealised benthic biofilm model is incorporated in a 1D morphodynamic model of tide-dominated channels. This study investigates the effect of a range of environmental and biological conditions on biofilm growth and their feedback on the morphological evolution of the entire intertidal channel. By carrying out a sensitivity analysis of the bio-morphodynamic model, parameters like (i) hydrodynamic disturbances, (ii) seasonality, (iii) biofilm growth rate, (iv) temperature variation and (v) bio-cohesivity of the sediment are systematically changed. Results reveal that key parameters such as growth rate and temperature strongly influence the development of biofilm and are key determinants of equilibrium biofilm configuration and development under a range of disturbance periodicities and intensities. Long-term simulations of intertidal channel development demonstrate that the hydrodynamic disturbances induced by tides play a key role in shaping the morphology of the bed and that the presence of surface biofilm increases the time to reach morphological equilibrium. In locations characterised by low hydrodynamic forces, the biofilm grows and stabilises the bed, inhibiting the transport of coarse sediment (medium and fine sand). These findings suggest biofilm presence in channel beds results in intertidal channels that have significantly different characteristics in terms of morphology and stratigraphy compared abiotic sediments. It is concluded that inclusion of bio-cohesion in morphodynamic models is essential to predict estuary development and mitigate coastal erosion.

1 Introduction

Tidal inlets are some of the most sensitive systems to sea-level rise and environmental change. Their morphology is shaped and influenced by tides, waves, river discharge, and associated sediment supply of marine and riverine sands and muds (Corenblit et al., 2007; De Haas et al., 2018; Noffke, 2008). The availability of nutrients and sediment from the surrounding area in combination with dynamic environmental conditions provide a favourable setting for numerous aquatic species, making them one of the most ecologically important environments (Meire et al., 2005). Even though strongly driven by abiotic processes, biotic processes can determine the geomorphological evolution of intertidal areas (Defew et al., 2002; Malarkey et al., 2015; Parsons et al., 2016; Vignaga et al., 2013). In order to manage these systems and adapt for future changes, there is the need for models that are able to incorporate the role of bio-cohesion on geomorphology. Those currently available are not yet robust enough to predict with confidence very far into the future. Consequently, understanding the interactions between hydrodynamics, sediment erosion and deposition, and biological communities becomes crucial for the sustainable management of estuaries and intertidal environments.

Biological activity on the seabed is known to have a significant influence on the bed composition and dynamics of cohesive and non-cohesive sediment at both small spatial and temporal scales (Decho, 2000). The presence of benthic microorganisms and the substances that they secrete strongly mediate the physical behaviour and functionality of the depositional system, influencing the structure and behaviour of sedimentary habitats, acting as ecosystem engineers (Paterson, 1997; Paterson et al., 2018). Microphytobenthos (MPB) is an assemblage of microbial cells, e.g. diatoms, cyanobacteria and heterotrophic bacteria, aggregated within a gel matrix composed of a mixture of lipids, proteins and polysaccharides, known as extracellular polymeric substances (EPS), that form benthic biofilms in intertidal and subtidal areas (Austin et al., 1999; Decho, 2000; Paterson et al., 1994; Tolhurst et al., 2002; Underwood and Paterson, 1993). Biofilms composed of MPB and EPS are ubiquitous in aquatic sediments (sand and mud) from shallow fluvial systems to continental shelves within the photic zone (Cahoon, 1999), even under physical disturbance from flow (Hope et al., 2020; Pinckney et al., 2018). While prevalence and patchiness can be greater on intertidal muddy flats, biofilm distribution in sandier intertidal and subtidal channels can be more homogenous, as seen in the Western Scheldt (Daggers et al., 2020).

It has been shown that secreted EPS is crucial in the adhesion and cohesion of the substratum and sediment particles, and it can act as a protective layer at the bed surface, reducing the bed roughness and significantly influencing the erosion and deposition of sediment particles by raising the sediment erosion threshold due to cohesion (Tolhurst et al., 2002, 2006, 2009; Paterson et al., 2018; Hope et al., 2020). This promotes the sedimentation of fine-grained particles and subsequently stimulates biofilm growth (Weerman et al., 2010) as nutrients are supplied to the bed. Microbial production of EPS is not only influenced by nutrient availability but can also be stimulated by exposure to contaminants such as heavy metals and nanoparticles (Ruddy et al., 1998; Lubarsky et al., 2010). Even at low EPS content (Tolhurst et al., 2002), both EPS concentrations (quantity) and EPS components (quality) play important roles in the binding effect of sediment particles, increasing the critical threshold for erosion and “biostabilisation” (Paterson et al., 1989; Tolhurst et al., 2002; Widdows et al., 2000) and thereby reducing sediment resuspension and bed erosion (Lubarsky et al., 2010; Malarkey et al., 2015; Parsons et al., 2016). This allows for the spatial development of biofilms and stabilisation across large geomorphological features (Weerman et al., 2010; Friend et al., 2008). By reducing the concentration of resuspended fine sediment and consequently the turbidity of the water column, biostabilisation improves light penetration to the sediment surface, creating a positive feedback to the biofilm community and more growth. Biostabilisation also limits the resuspension of coarse particles that, by moving, could cause abrasion to the biofilm layer and the removal of large sections of biofilm from the bed (Lanuru et al., 2007). Further, the stabilisation of the water–sediment interface by benthic biofilm is important for the regulation and bentho-pelagic exchange of carbon, nitrogen and oxygen with the substrate (Cahoon, 1999) and subsequently the transfer of energy and resource to adjacent habitats (Savage et al., 2012).

These processes are complicated by the presence of benthic bioturbators that disrupt and graze on MPB, and they can have a high impact on mudflat morphology because they can physically destabilise the bed (e.g. de Deckere et al., 2001; Brückner et al., 2021) and trigger sediment resuspension that is otherwise stabilised by diatoms. Furthermore, the establishment of biostabilisers might be affected by sediment destabilisation and seed predation from bioturbators (Cozzoli et al., 2019). In turn, bioturbators organically enrich the sediment via biodeposition, which can promote the MPB growth (e.g. Andersen et al., 2010; Donadi et al., 2013), and biostabilisers can modify the hydrodynamics and sediment properties around them (Brückner et al., 2020), impacting the size and density of bioturbators communities (Walles et al., 2015).

While microbially produced EPS is more abundant in cohesive sediment (muddy bed), studies have shown that EPS production by bacteria and microphytes can also play a significant role in non-cohesive and mixed sediment substrates by hindering bedform development and inhibiting erosion (Malarkey et al., 2015; Parsons et al., 2016; Chen et al., 2017; Hope et al., 2020). The influence of benthic biofilms and EPS on sediment erosion is widely recognised and characterised across different sedimentary habitats (e.g. Paterson, 1989; MacIntyre et al., 1996; Marani et al., 2010; Malarkey et al., 2015; Hope et al., 2020; Chen et al., 2021), but few numerical studies account for these processes. The exclusion of bio-cohesion and biostabilisation effects makes it difficult for predictive models of sediment stability to be sufficiently accurate. This is primarily due to the difficulty of simultaneously coupling the physical, biological and biodiversity components. Seasonal changes in environmental conditions and grazer communities can mediate biofilm grow rate (Underwood and Paterson, 1993; Montani et al., 2003; Zhang et al., 2021; Daggers et al., 2020; Brückner et al., 2021), but interannual changes in key biota, through their influence on sediment erosion, and the consequences for intertidal ecology and morphology, can also be driven by climatic factors such as changes in water and sediment temperature (Marani et al., 2007, 2010; Mariotti and Fagherazzi, 2012), which is strongly regulated by the light availability due to the turbidity of the water column. Quantifying and understanding these benthic processes in order to parameterise them into mathematical models is critical for providing insights into the relative importance of biological and physical factors in sediment erosion and accretion in the intertidal zone.

A range of hydro-morphodynamic models have attempted to parameterise eco-engineering processes on varying spatial and temporal scales (Brückner et al., 2020, 2021; Coco et al., 2013; Le Hir et al., 2007; Mariotti and Canestrelli, 2017). While extensive field and flume studies are available in the literature on the effect of MPB and faunal on sediment erosion (Le Hir et al., 2007; Cozzoli et al., 2019), the main challenges in modelling these types of environments are the complexity of the interaction between the different biotic and abiotic contributors, the time and spatial scales, and the fact that variation in sediment stability might reflect site-specific differences (Le Hir et al., 2007; Pivato et al., 2019). In fact, the interactions between these processes are strongly regulated by spatio-temporal conditions (e.g. Widdows et al., 2000; van de Lageweg et al., 2018; Paterson et al., 2018; Best et al., 2018; Cozzoli et al., 2019), making it difficult for predictive models of sediment stability to make generalities from site-specific findings and to be sufficiently accurate.

For the first time, this study investigates the effect of the environmental conditions, such as temperature, seasonality and sediment rheology, on biofilm growth and its feedback to the bed stability and morphological evolution over an entire intertidal channel. The main objective was to investigate and define the key parameters of the biofilm development model that influence the intertidal channel morphology. The combined effect of temperature, biofilm growth rate and surface biofilm removal due to tidal dynamics is investigated for different scenarios.

A one-dimensional eco-morphodynamic shallow-water model is implemented and tested in this study to investigate the effect of biostabilisation due to the presence of surface biofilm. The model accounts for the effect of tidal oscillation on a non-uniform non-cohesive sediment channel subject to tidal fluctuations at the ocean boundary, and it allows us to store the information of the stratigraphy of the deposit emplaced. The biofilm logistic growth model accounts for the effect of hydro-climate variation on the biofilm development, such as temperature changes and carpet-like erosion, as these are key factors controlling biofilm development (Pivato et al., 2019). The model is tested for different benthic biofilm growth rates. Biostabilisation from presence of surface biofilms is implemented in the 1D morphodynamic shallow-water model assuming a linear relationship that correlates the amount of biofilm biomass with the increase in the sediment critical shear stress for erosion (Le Hir et al., 2007). The model is applied to an initial flat bed to investigate the implications of different sediment temperatures, representative of different climate scenarios, and different sediment rheology on the channel development.

Bio-sedimentology summary of processes and controls

Since the living and abiotic elements vary temporally and spatially, it is not surprising that the functions and importance of these various factors in determining sediment stability also vary (Defew et al., 2002; Friend et al., 2003; Paterson et al., 1994; Riethmüller et al., 2000; Underwood et al., 1995; Yallop et al., 1994b). Benthic biofilms change the fundamental properties of sediment and bed substrate: when biofilm develops on the bed surface, it acts as a protective skin on the sediment surface inhibiting entrainment (Paterson et al., 2000) with greater volumes of biofilm required to stabilise sandier beds (Hope et al., 2020).

Numerous studies in marine intertidal environments show a positive correlation between sediment stability in terms of critical shear stress for erosion (τ_bc) and EPS components of biofilm. Although it is EPS that stabilises the bed, not the MPB per se, chlorophyll-a (Chl-a), a proxy of living MPB biomass, provides a good approximation of biostabilisation potential (Defew et al., 2002; Paterson et al., 2000; Riethmüller et al., 2000; Haro et al., 2022; Le Hir et al., 2007; Yallop et al., 1994a). Chl-a is often the preferred measurement, due to its ecological significance and the fact that it is easy to evaluate (both in the field and by optical remote sensing) (Andersen, 2001; Le Hir et al., 2007), but Chl-a–stability relationships can often be weak, emphasising the complexity of this phenomenon and that important interactions are being missed. Hydro-sedimentary processes, modulated by the shear stresses at the bed due to tidal and waves, regulate the biofilm resuspension process and its flux in the water column. The erosion fluxes depend on the bed erodibility, described by the resistance of the sediment to be eroded (Orvain et al., 2014). Changes in bed erodibility, which vary largely in space and time, is the result of a complex interaction between sediment properties, bioturbation activities, grazing, biofilms deposition, reseeding and growth rate (Wood and Widdows, 2002; Thrush et al., 2012; Cozzoli et al., 2019). Due to the complexity of these systems, multiple factors play a relevant role in defining a relationship between critical shear stress for erosion and Chl-a or EPS. There is non-standard relationship but a general tendency for shear stress to increase with Chl-a content (Paterson et al., 1994; Yallop et al., 1994a; Underwood et al., 1995; Riethmüller et al., 2000; Defew et al., 2002; Friend et al., 2003; Le Hir et al., 2007; Righetti and Lucarelli, 2007; Fang et al., 2014), and often results are site specific (Riethmüller et al., 2000; Le Hir et al., 2007; Katz et al., 2018). There is thus a fundamental need for a broad-scale bio-morphodynamic approach to synthesise the general effects across habitats modulated, for example, by the distribution of benthic macrofauna, the sediment types, the water content or the tidal range.

The development of biofilm is controlled by various sedimentary characteristics, biogeochemical drivers and light-related photosynthesis parameters (e.g. optimum and maximum temperature for MPB photosynthesis, light saturation parameter) and their spatio-temporal variability (MacIntyre et al., 1996; Pivato et al., 2019, Savelli et al., 2020), the availability of nutrients (Hillebrand and Sommer, 1997), hydrodynamic disturbances such as currents and waves (Mariotti and Fagherazzi, 2012; Tolhurst et al., 2009, 2006), and grazing benthic macrofauna (Hillebrand et al., 2000; Montserrat et al., 2008). Even when biofilms are removed during tidal inundation, the remaining MPB community can quickly re-establish itself, depending on the prevailing conditions, with a subsequent increase in biostability, as cell numbers increase and EPS secretions once again build up (Valentine et al., 2014; Chen et al., 2019; Hope et al., 2020). The prevailing environmental conditions can significantly influence biostabilisation processes, with the temperature of the water and underlying sediment layers exerting a major influence on chemical and biological processes and kinetics including benthic nutrient cycling (Smith, 2002; Pivato et al., 2019). In shallow-water environments, the energy exchange at the water–sediment interface, the turbidity of the water column and the light reaching the bed surface are crucial to appropriately describe the sediment temperature (Pivato et al., 2018, 2019). Experimental studies of the response of biofilm communities to water warming have shown faster biofilm growth with increases in temperature (Majdi et al., 2020). Therefore, seasonal temperature changes influence the resistance to erosion (Thom et al., 2015); for example, increases in temperature during the spring promote photosynthesis, leading to higher Chl-a concentrations and biofilm growth and greater biostabilisation (Underwood and Paterson, 1993; Savelli et al., 2018; Pivato et al., 2019; Haro et al., 2022).

When sediments are covered by biofilm, the entrainment process can occur as sediment–biofilm coated particles (flocs) or via the resuspension of sediment–biofilm aggregates (biofilm failure due to carpet-like erosion) (Shang et al., 2014; Fang et al., 2016; Fang et al., 2017). Resuspended biofilm-coated particles can be transported as bed load, and deposited under the different settling velocities, governed by the sediment shape and size, amount of biofilm and density of the particles or flocs (Koh et al., 2007). Hydrodynamic disturbances from currents, tides and waves play a cardinal role, eroding the biofilm and eventually detaching it from the sediment surface. Once the protective biofilm is broken or removed, the underlying clean sediment is exposed, the erodibility of which is regulated by the characteristic sediment grain size of the substrate (Defew et al., 2002; Le Hir et al., 2007; Sutherland et al., 1998).

2 Methodology

A 1D morphodynamic model for tide-dominated channels implemented with a function that describes the surface biofilm growth was used to determine the relative importance of different bio-physical factors on the development of an intertidal channel longitudinal profile and stratigraphy. The abiotic physical processes included in this study are tidal currents, sediment erosion, transport and deposition. The model takes into account the dynamics of biofilm development and its feedback on the erosional and depositional sediment transport processes. The model is based on the one-dimensional shallow-water equations (1D-SWE) for the flow mass, sediment and momentum conservation, modified according to Defina (2000) to account for partially dry areas, such as the beach that can be formed at the landward boundary of the model domain (Fig. 1). The model is implemented with a procedure that stores and accesses the information of the grain size of the stratigraphy of the deposit.

Figure 1 Schematic representation of the model geometry. The ξ and η are the water surface elevation and channel bed elevation. At the beginning of the simulation the bed is assumed horizontal (η₀), the water depth is H and the mean water surface elevation is ξ₀. ω_t represents the tidal amplitude, and the water surface elevation at the ocean boundary (x=L) is ξ_d.

2.1 A 1D SWE model for tidal channel accounting for partially dry areas

The shallow-water equations model (Chaudhry, 2008) is used to describe temporal and 1D spatial variation in idealised tidally dominated channel reaches (Fig. 1). The domain is bounded by the ocean, where the tidal oscillations are modelled as a sine curve with amplitude α_t and period ω_t. Input of riverine water and sediment at the landward boundary (Lanzoni and Seminara, 2002) and interaction of the channel with tidal flats and intertidal areas (Todeschini et al., 2008) are assumed to be negligible.

The shallow-water equations, modified by Viparelli et al. (2019) according to Defina (2000), account for the partially dry areas such as when the channel bed is only periodically submerged. Defina (2000) derived the two-dimensional shallow-water equations by averaging the Reynolds equations over the bottom irregularities and then integrated them for mass and momentum conservation in the direction normal to the channel bed. The one-dimensional form is obtained by integrating the equations in the transverse direction (Viparelli et al., 2019), giving the following equation:

$$\begin{array}{}\text{(1)}& \left\{\begin{array}{l}{F}_H\frac{\partial A_{\mathrm{i}}}{\partial t}+\frac{\partial Q}{\partial x}=\mathrm{0},\\ \frac{\partial Q}{\partial t}+\frac{\partial }{\partial x}\left(\frac{Q^{\mathrm{2}}}{A_{\mathrm{c}}}\right)+\mathit{\rho }g{A}_{\mathrm{c}}\frac{\partial \mathit{\xi }}{\partial x}+\frac{{\mathit{\tau }}_{\mathrm{b}}}{\mathit{\rho }}\mathit{\chi }=\mathrm{0},\end{array}\right.\end{array}$$

where Q is the volumetric flow discharge, A_c is the cross-sectional area and ρ is the water density. The cross-sectional area averaged over bed irregularities A_i is equal to (W⋅ξ). The wet fraction of the channel bed (F_H) is computed as function of the characteristic length scale of the bed irregularities (a_r, assumed equal to 1 cm), the effective flow depth (Y) and the average bed shear stress (τ_b) acting over the wetted perimeter χ (see Viparelli et al., 2019, for further details of the 1D morphodynamic model). The model validation is presented in Appendix A and shows that the model can reasonably capture the magnitude and timing of the bed changes. The numerical model is demonstrated to be second-order accurate, and model parameters are reported in Table 1.

Table 1 Parameters for the shallow-water model.

2.2 Sediment transport model

A sediment transport model is incorporated to describe well-mixed, non-cohesive sediment transport and the coupled morphodynamics (Viparelli et al., 2019). The total volumetric bed material load (Q_b) is calculated as the contribution of bed load and suspended load. The equations to compute the bed load and the suspended load implemented in the model have been selected to let the direct correlation between the amount of biofilm biomass on the bed and the updated critical shear stress for sediment motion that results in biostabilisation.

The bed load is computed using the Ashida and Michiue relation, while the McLean formulation is used to model the entrainment of sediment in suspension. The total material load (Q_b) is the sum of the contribution of bed load and suspended load, summed over all the grain sizes, and the volume fraction content of sediment with characteristic diameter D_i can be computed as $(Q_{\mathrm{b},\mathrm{bi}}+Q_{\mathrm{b},\mathrm{si}})/Q_{\mathrm{b}}$.

The equation for the conservation of the sediment material, coupled with a procedure to store the information of the stratigraphy of the deposit, is solved to compute the temporal evolution of the bed profile (η) and the spatial distribution of the sediment size (Viparelli et al., 2010). To solve this equation, according to the Hirano active-layer approximation (Hirano, 1971), the deposit can be divided into two regions, the active layer and the substrate. The active layer (L_a) is the topmost part of the deposit where the sediment particles can interact with the flow and it is assumed to be well mixed, meaning that the grain size distribution of the sediment on the active layer can change in space and time but is assumed constant in the vertical direction. The substrate (η−L_a) is located below the active layer and does not interact with the flow; the sediment fraction in the substrate varies in space but not in time. Exchange between the substrate and the active layer occurs in the case of aggradation and degradation. During aggradation the distance between the substrate and the active layer increases, and layers can be added to the grid for the storage of the newly deposited sediment. The grain size distribution of the antecedent storage layer is computed as a weighted average, while the sediment composition of the new storage layers has the same grain size distribution as the newly deposited material. Interested readers may refer to Viparelli et al. (2010) for further details about the deposit storage procedure.

2.3 Biofilm-dependent erodibility

The novelty of this work is the implementation of a 1D morphodynamic model for intertidal channels with a biofilm growth model that accounts for the effect of seasonality on sediment temperature and light. This study aims to understand the general behaviour of the system and investigate the sensitivity of the biofilm model parameters on the channel development process; hence, the assumption of spatially homogenous biofilm or constraint on the development of biofilm only in the cells where the water depth is smaller than 0.05 m are reasonable. Once the biofilm biomass is estimated according to the biofilm growth model, the critical shear stress for erosion is updated to account for the biostabilisation. According to Le Hir et al. (2007), the increase in critical shear stress is assumed to be proportional to the biofilm biomass available on the bed (B, measured in mg Chl-a m⁻²):

$$\begin{array}{}\text{(2)}& {\mathit{\tau }}_{\mathrm{bc}}={\mathit{\tau }}_{\mathrm{bc},\mathrm{0}}+\mathit{\alpha }B.\end{array}$$

Here (τ_bc,0) is the critical shear stress for clean sediment. The updated value for the critical shear stress is used in the bed load and suspended load equations to correlate the sediment mobility with the amount of surface biofilm. The time evolution of biofilm biomass (B) is estimated by a simplified model proposed by Mariotti and Fagherazzi (2012) that assumes a logistic growth function for the biofilm biomass:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{d}B}{\mathrm{d}t}}=P^{\mathrm{B}}B{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+K_{\mathrm{B}}B}}-\mathit{\epsilon }(B-B_{\min })-E,\end{array}$$

where P^B is the effective maximum growth rate, and K_B is the half-saturation constant, which represents the biofilm concentration at which it has reached half of the maximum growth rate, and this term accounts for the effect of density limitation. The second term of the equation accounts for the chronic and self-generated biofilm detachment (ε being the global decay parameter), which is not associated with the simulated hydrodynamics (e.g. senescence, heterotrophic processes, benthic macrofauna grazing), and B_min is the amount of background biofilm biomass that allows the recolonisation after removal. Starting from a background value for the surface biofilm (B_min), the biofilm grows only if there are no disturbances limiting the establishment of biofilm. The last term of Eq. (3) takes into account the effect of extremely high-intensity flow events (E) that are able to mobilise the bed and completely remove the surface biofilm, exposing the clean sediment underneath. The reference values for the parameters of the biofilm growth function (Table 2) are based on field observations, assuming that in equilibrium conditions the surface biofilm biomass is equal to 200 mg Chl-a m⁻² (Mariotti and Fagherazzi, 2012; Le Hir et al., 2007), which is a value commonly found in intertidal environments in temperate areas (Mariotti and Fagherazzi, 2012; Le Hir et al., 2007). At the initial stages, the growth of undisturbed biofilm is approximately exponential. As the saturation begins, it slows to linear until it reaches maturity when the growth stops. The amount of biofilm on the bed surface remains constant for the entire duration of the simulation, asymptotically reaching an equilibrium condition (Fig. 2a).

Figure 2 Biofilm development over time. Biofilm growth patterns in cases of (a) annual undisturbed growth following the logistic growth function, (b) growth affected by the variation in sediment temperature due to seasonality over a year-long simulation, and (c) growth affected by carpet-like erosion.

The biofilm model has been implemented to account for the seasonal cycle of temperature and light as proposed by Pivato et al. (2019), based on the vertical energy transfer within the water–sediment continuum. This sediment temperature model simulates natural conditions that regulate the development of biofilm, such as the effect of winter conditions that limit the growth of MPB, leading to lower surface sediment biostabilisation and resistance to erosion compared to late spring, summer and early autumn (Fig. 2b), as also confirmed by in situ observations (Friend et al., 2003). The MPB photosynthesis and biofilm development are both strongly influenced by the seasonal changes in sediment temperature and light availability, which are controlled by the water depth and turbidity. The sediment temperature model implemented in this study accounts for the effect of seasonality, and it is based on sediment temperature parameters that are based on temperate environments as proposed by Guarini et al. (2000), Pratt et al. (2014) and Pivato et al. (2019) (Table 2). The maximum growth rate of MPB (P^B) is computed according to Guarini et al. (2000):

$$\begin{array}{}\text{(4)}& P^{\mathrm{B}}=P_{\max }^{\mathrm{B}}\tanh (H_{\mathrm{res}}/E_k).\end{array}$$

The light saturation parameter E_k (W m⁻²) is assumed to be constant. The light availability (H_res) is represented by the residual solar radiation reaching the bed and not reflected by the water surface albedo (A= 0.04), and it is computed as follows:

$$\begin{array}{}\text{(5)}& H_{\mathrm{res}}=R_{\mathrm{0}}{\mathrm{e}}^{-\mathit{\lambda }Y};R_{\mathrm{0}}=\left(\mathrm{1}-A\right)R_{\mathrm{sun}}.\end{array}$$

The extinction coefficient λ represents the capability of the water column to absorb the solar radiance, describing the average effect of the turbidity in the water column (Y being water depth) on radiative transfer, and R_sun is the solar radiation. $P_{\max }^B$ (h⁻¹) represents the growth rate under light saturation conditions. This parameter varies in time and depends on the surface sediment temperature (T_s0) according to:

$$\begin{array}{}\text{(6)}& \left\{\begin{array}{ll}\mathrm{if}{T}_{\mathrm{s}\mathrm{0}}<T_{\max }:& P_{\max }^B=P_{\max }{\left(\frac{T_{\max }-T_{\mathrm{s}\mathrm{0}}}{T_{\max }-T_{\mathrm{opt}}}\right)}^{\mathit{\beta }}\\ & \exp \left[\mathit{\beta }\left(\mathrm{1}-\frac{T_{\max }-T_{\mathrm{s}\mathrm{0}}}{T_{\max }-T_{\mathrm{opt}}}\right)\right]\\ \mathrm{if}{T}_{\mathrm{s}\mathrm{0}}\ge T_{\max }:& P_{\max }^B=\mathrm{0}.\end{array}\right.\end{array}$$

This is the function of the optimal and maximum temperature for photosynthesis (T_opt= 25 ^∘C, T_max= 38^∘C), where the shape factor (β) is site dependent. The parameter P_max represents the maximum value for $P_{\max }^B$ and is site and time dependent. The seasonal changes in the sediment temperature modulate the amount of biofilm biomass and, as a consequence, the biostabilisation of the bed (Fig. 2b). For simplification purposes, in this study the sediment temperature will be assumed following a parabolic trend during the 1-year interval (continuous blue line in Fig. 2b, Pivato et al., 2019). Biomass increases exponentially at the beginning of the year, reaching its maximum when the sediment temperature is equal to the optimal temperature for photosynthesis (dotted orange line in Fig. 2b, T_opt) during spring and autumn. As the sediment temperature increases during the summer months (continuous orange line in Fig. 2b), photoinhibition can occur and the biofilm biomass decreases (blue line in Fig. 2b), reaching a local minimum when the sediment temperature is at its maximum and close to the maximum temperature for photosynthesis (dashed orange line in Fig. 2b, T_max). The growth rate during these months is still sufficient to enable a fast recovery of the biofilm. As light and sediment temperature decrease during the winter season, the environmental conditions are less favourable for the growth of biofilm. In cases when the availability of light at the bed is limited and the sediment temperature is lower than the optimal temperature for photosynthesis, surface biomass decreases.

The quantification of the removal of the surface biofilm by intense hydrodynamic forces (carpet-like erosion) occurs in a very short period of time, and it can thus be considered instantaneous, and the catastrophic erosion (E) is as follows:

$$\begin{array}{}\text{(7)}& E\left(B,t\right)=E_{\mathrm{0}}\left(B\right){\sum }_i\mathit{\delta }\left(t-t_i\right),\end{array}$$

where δ is the Dirac function, t_i is the time of the detachment and E₀ is the intensity of the extreme event, assumed to be an “all-or-nothing” process that can be described as a function of the shear stresses acting on the bed (τ):

$$\begin{array}{}\text{(8)}& E_{\mathrm{0}}=\left\{\begin{array}{ll}\mathrm{0}\mathit{\tau }\le & {\mathit{\tau }}_{\mathrm{bc}}\\ B-B_{\min }& \mathit{\tau }>{\mathit{\tau }}_{\mathrm{bc}}\end{array}\right..\end{array}$$

In the case that shear stresses due to the hydrodynamic forces (τ) are smaller than or equal in value to the sediment critical shear stress for erosion, there is no disruption of the surface biofilm. In the case that the stress on the bed exceeds the critical value for erosion (τ_bc), the biofilm is eroded and is reduced to the background value B_min, which allows for the establishment and growth of biofilm (Fig. 2c). When biofilm is removed from the bed surface as a carpet-like erosion, the resistance of the bed reduces to a minimal value (Fig. 2c) under the assumption of linear relationship between surface biofilm biomass and critical shear stress for erosion (Le Hir et al., 2007). This simplified model assumes that the erosion is on the order of millimetres to centimetres in the case of extreme hydrodynamic events, which is much larger than the thickness of the biofilm (micrometres to millimetres). The range of values found in the literature and the reference values selected here are reported in Table 2.

Table 2 Parameter ranges found in the literature (Mariotti and Fagherazzi, 2012; Pivato et al., 2019) and parameter values used as a reference in the model.

The sediment mixture used for the simulations is characterised by median diameter D₅₀= 0.323 mm and geometric mean sediment grain size D_g = 0.303 mm.

By changing the biofilm model parameters within the range found in the literature, this study investigates the sensitivity of the key biofilm model parameters to the morphological evolution of an intertidal channel (Table 3).

Firstly, a sensitivity analysis of the biofilm stability is presented under different hydrodynamic disturbances characterised by periodicity (T) and intensity (τ₀, shear stress). The sensitivity analysis is performed for all model parameters and within the range of values suggested in the literature (between 0.6 and 1.4 times the reference value suggested in the literature, Table 2) by systematically changing the periodicity and intensity of the hydrodynamic disturbances to evaluate under which conditions the biofilm is stable. Then, with the objective to test the biofilm stability modulated by the effect of seasonality, a sensitivity analysis is performed for the maximum and minimum values of the parameters that are strongly affected by seasonality, like the biofilm growth rate (P^B) and the sediment temperature (T+5 ^∘C, T −5 ^∘C).

A second set of sensitivity analysis tests aim to understand the bio-modulation of channel morphodynamic evolution for an idealised channel characterised by a semidiurnal tidal environment. The results show the comparison between the impact on biostabilisation for uniform and spatially varied surface biofilm under the effect of carpet-like erosion, which is regulated by the periodic hydrodynamic disturbance changes in the water level at the seaward boundary due to the tidal forces, and seasonality. In the case of a combined effect of these conditions (seasonality and carpet-like erosion), it has been investigated how channel morphology changes with the variation in (i) biofilm growth rate, (ii) sediment temperature and (iii) sediment bio-cohesivity.

The numerical simulations have been performed in the absence of an imposed input of sand from the ocean and without riverine water and sediment at the landward boundary. An idealised 25 m long channel with constant width equal to 30 cm has been assumed.

Table 3 Summary of the simulations performed in this study.

3 Results

3.1 The control of hydrodynamic disturbances and biofilm model parameters on biofilm stability

The sensitivity analysis of the biofilm model parameters (Eqs. 2–6) have been investigated by systematically changing the intensity and the periodicity of the disturbances to find the hydrodynamic conditions at which the status of the biofilm changes from stable to detached (Fig. 3). Mariotti and Fagherazzi (2012) have shown that biological biofilm growth parameters (P_B, K_b, ε) can affect the stability of surface biofilm in terms of the resistance of biofilm to being eroded from the bed by high-intensity hydrodynamic forces, i.e. tides. The bio-cohesivity parameter that correlates the presence of surface biofilm with the increase in the bed resistance (α) ends up being important in the determination of the equilibrium configuration (steady biofilm). The parameters that describe the biofilm maximum growth rate such as the dimensionless shape factor (β), the surface albedo (A), the extinction coefficient, which is proxy for the water column turbidity (λ), and the light saturation parameter (E_k), do not influence the growth of biofilm under the effect of different hydrodynamic disturbances.

Figure 3 Sensitivity analysis of the model parameters for the determination of the equilibrium configuration. The dashed black line represents the simulation considering the reference value (RV, Table 2), dark green lines represent conditions of the parameters under examination above the reference value (1.2 and 1.4 times the reference value, respectively), and light green lines represent conditions below the reference value (0.8 and 0.6 times the reference value). The area below the curves represents conditions of stable biofilm.

The combined effect of seasonality of sediment temperature and hydrodynamic events are reported in Fig. 4 under a set of different temperature-influenced scenarios that are intended to simulate the changes in nutrient availability in the water (growth rate parameter) and the long-term variation in temperature. The reference profile for the development of biofilm is reported in Fig. 2b ($T_{\mathrm{s}\mathrm{0},\max }=$ 32 ^∘C, P_B= 1.068 d⁻¹). The sensitivity analysis is carried out for a year-long cycle. The intensity and periodicity of the hydrodynamic conditions are selected from the previous analysis according to what has been observed by Mariotti and Fagherazzi (2012). High-intensity and low-disturbance periodicity events (case 1 in Fig. 4, T = 15 d, τ₀= 1.5 Pa) are assumed to allow for the growth of biofilm under reference values for the biofilm model parameters, while under frequent and weak disturbances (case 2 in Fig. 4, T = 5 d, τ₀= 0.5 Pa) the biofilm is not fully established on the bed.

Figure 4 Effect of seasonality and hydrodynamic forces on the evolution of surface biofilm biomass. The evolution of the temperature of the sediment at the bed is simulated for a period of 1 year under rare (every 15 d) and strong (1.5 Pa) hydrodynamic disturbances. The effect of different values of the growth rate parameter (a, b c) and sediment temperature are investigated (d, e). The evolution of the temperature of the sediment at the bed is simulated for a period of 1 year under frequent (every 5 d) and weak (0.5 Pa) hydrodynamic disturbances. The effect of different values of the growth rate parameter (f, g, h) and sediment temperature are investigated (i, j).

With high-intensity and rare events (case 1, T= 15 d, τ₀ = 1.5 Pa) and small values of the grow rate parameter ($P_{\max }^{\mathrm{B}}$ = 0.0078 and 0.5617 d⁻¹), the new settled biofilm is periodically detached by the disturbances (Fig. 4a and b). Biofilm grows during the time span between two consecutive events (Fig. 4a–b), but it is destroyed every time a significant hydrodynamic event occurs. The increase in sediment resistance is not enough to prevent the erosion caused by the high-intensity events. A further increase in the maximum growth parameter ($P_{\max }^{\mathrm{B}}=$ 1.068 d⁻¹) results in a more rapid growth and establishment of biofilm. During the initial and final months of the simulated year (January to mid-March and after mid-November) the biofilm is periodically removed because the temperature of the sediment inhibits the development of biomass (Fig. 3c). During the spring and summer months the combination of the temperature conditions and high growth rate promote the development of stable biofilm that is able to resist the periodic disturbances.

The biofilm biomass profile under rare and intense hydrodynamic disturbances, with a variation in the annual sediment temperature of ±5 ^∘C compared with the previous simulation ($P_{\max }^{\mathrm{B}}=$ 1.068 d⁻¹), is reported in Fig. 4d–e. Cases with an increase in the sediment temperature of the profile, analogously to the previous cases, show a slowdown in the development of biofilm in winter and autumn. Furthermore, during the summer period (June to September) the sediment temperature increases above the optimal temperature for photosynthesis (T_opt= 25 ^∘C), resulting in a drop in EPS production, a reduction in the bed stabilisation (linear decrease in the critical shear stress for erosion) and an increase in vulnerability to hydrodynamic disturbances (Fig. 3d). The effect of an overall annual reduction of the sediment temperature on biofilm is shown in Fig. 4e. In this case the biofilm is more vulnerable to the disturbances at the beginning and end of the simulated year compared with the profile in Fig. 4c, due to the fact that temperature conditions further from the optimal temperature for photosynthesis reduce the rate of development of biofilm.

When the biofilm growth rate parameter is low (Figs. 4f and 3g) under frequent and weak disturbances (case 2, T= 5 d, τ₀= 0.5 Pa), biofilm is periodically detached and it cannot establish during the entire simulated year. An increase in the growth rate parameter shows that in summer the biofilm can establish and cover the bed surface until the end of the year, even though the biomass decreases in autumn and winter ($P_{\max }^{\mathrm{B}}=$ 1.068 d⁻¹, Fig. 4h), unlike in cases of strong disturbances (Fig. 4c). The increased amount of biofilm enhanced the bed stabilisation, also inhibiting the erosional behaviour under further disturbances. Comparing Fig. 4h with the case in which the annual sediment temperature is increased (Fig. 4i) or decreased (Fig. 4j) by 5 ^∘C shows that an increase in temperature would decrease the amount of biomass at the bed, whereas a decrease in sediment temperature would not allow biofilm to establish because it would be constantly destroyed by the frequent disturbances.

It is reasonable to conclude that the presence of consolidated biofilm able to stabilise the bed not only depends on the intensity and the frequency of the disturbing events but also on sediment temperature and seasonal parameters. The amount of biofilm biomass on the bed surface plays a significant role in defining areas of erosion, even under the same hydrodynamic conditions (Mariotti and Fagherazzi, 2012; Hope et al., 2020).

3.2 Effect of seasonality and carpet-like erosion

This section explores the effects of seasonality and carpet-like erosion on the morphological evolution of an intertidal channel. Two main cases are considered in this study for the spatial distribution of biofilm. One case assumes that biofilm is uniformly distributed in the entire computational domain (centre column in Fig. 5) to also explore the case of biofilm development in the deepest portion of the channel. In fact, biostabilising organisms are found along the entire tidal range, from intertidal and subtidal areas to shellfish reefs and the continental shelf, as has been suggested in the literature (Cahoon, 1999; Pinckney, 2018; van de Vijsel et al., 2020). For the second set of simulations, the biofilm is assumed to grow in turbid systems, where light attenuation would prevent substantial growth of surface biofilm due to the limited availability of light for photosynthesis processes. The development of biofilm is therefore limited in locations where the water depth is below 0.05 m, which corresponds to the portion of the channel that experiences the wet–dry transition according to the tidal amplitude range used for these simulations (right column in Fig. 5).

The model is applied to investigate the separate and combined effects of carpet-like erosion and seasonality. The effect of carpet-like erosion is modulated by changes in the water level at the seaward boundary due to the tidal forces creating periodic hydrodynamic disturbances. The bed evolution after 30 000 tidal cycles is reported in Fig. 5, where for the bed evolution profiles the dashed green line represents the initial bed and the dashed blue line the initial mean water surface. The profiles are compared to the final equilibrium bed profile in the case of clean sediment (Fig. 5a, red dashed line). In Fig. A3 the spatial distributions of the geometric mean diameter of the deposit at the end of each simulation are reported. The initial mean diameter of the transported sediment and the bed is 0.3 mm.

The reference case is characterised by clean sediment (Fig. 5a). The model initial conditions assume a flat bed, and there is a formation of an upstream-migrating shore at the landward boundary due to the effect of tides at the ocean boundary, creating an alluvial deposit characterised by sediment erosion at the ocean boundary (Lanzoni and Seminara, 2002; Tambroni et al., 2005; Todeschini et al., 2008; Viparelli et al., 2019). As the shoal reaches and is impeded at the landward boundary, a beach forms and grows until conditions of morphodynamic equilibrium are met (approximately after 20 000 tidal cycles).

First, the case of spatially uniform and stable biofilm on the bed surface for the entire duration of the simulation (Fig. 5b) is modelled. This resulted in less sediment mobility compared to the clear sediment scenario (Fig. 5a). The bed exhibits minor erosional behaviour at the ocean boundary, while at the land boundary the bed profile does not change in time and the bed is horizontal and stable (Fig. 5c). Figure 5d shows the bed elevation in the case of water depth being a constraint for the development of biofilm (H < 0.05 m). The bed is more mobile both at the ocean and in landward areas (Fig. 5d).

Figure 5 Bed evolution after 30 000 tidal cycles under different biofilm conditions. Panel (a) shows the bed evolution in the reference case with clean sediment. The panels in the left column represent the time evolution of biofilm during the simulation, showing stable biofilm (b) and the effect of seasonality (e) and carpet-like erosion (h). The rows represent, respectively, the bed evolution under these biofilm conditions. The panels in the central column show the bed evolution profile in the case of biofilm being uniformly distributed spatially. The panels in the right column show the bed evolution profile in the case of biofilm being developed only in locations where the water depth is smaller than 0.05 m. The bed profiles are compared with the final bed elevation of clean sediment at equilibrium (dashed red lines). The dashed blue and green lines represent the initial water surface elevation and the initial bed profile, respectively.

Considering the effect of seasonality (Fig. 5e), there are slight increases in sediment mobility seaward (Fig. 5f). However, assuming that biofilm is present only in shallow-water conditions (H < 0.05 m) results in an increase in bed mobility (Fig. 5g). Additionally, the bed needs a longer time to reach an equilibrium state, so it is reasonable to conclude that, even after 30 000 tidal cycles, the bed profile is still evolving.

The bed evolution in the case of surface biofilm being periodically removed by the tidal-induced stresses on the bed (carpet-like erosion, Fig. 5h) is shown in Fig. 5i. Hydrodynamic forces play a relevant role in shaping the bed, and the final profile is similar to the benchmark case, with erosion at the ocean boundary and deposition at the land boundary (Fig. 5a). The presence of biofilm hinders bed evolution, and more time is required to reach the equilibrium state. This is due to the periodic removal of surface biofilm due to the tidal forces, which causes periodic decreases in the bed critical shear stress for erosion and therefore biostabilisation. Assuming that biofilm is developed only under shallow-water conditions (water depth smaller than 0.05 m), the channel needs even more time to reach equilibrium (Fig. 5j).

For both the biofilm spatial distribution conditions investigated in this study (uniform biofilm and biofilm only in shallow-water conditions, H < 0.05 m), the combined effect of seasonality and carpet-like erosion results in a similar bed profile. In fact, the bed morphology is comparable in terms of bed elevation (Fig. 5k and l).

3.3 Effect of maximum biofilm growth rate parameter

Changes in the maximum growth rate in the biofilm development model (P_max, Eq. 6) result in a faster development of the biofilm, furthermore the peak of biofilm biomass appears in early stage of biofilm development (Fig. 6a). In the simulations showed above (Fig. 4), the maximum growth rate parameter has been assumed equal to 1.07 d⁻¹, which is a reference value that would give a biofilm biomass of 200 mg Chl-a m⁻² in steady-state conditions (model value, Table 2). Figure 6 shows the morphology and stratigraphy of the final bed after 30 000 tidal cycles under different values of the biofilm grow rate parameter: small (Fig. 6b and e, P_max= 0.0078 d⁻¹), medium (Fig. 6c and f, P_max= 0.56 d⁻¹) and large (Fig. 6d and g, P_max= 1.10 d⁻¹). Surface biofilm in these simulations has been assumed to develop only in locations where the water depth is smaller than 0.05 m.

Small or medium values for the maximum growth parameter for biofilm create a similar final longitudinal bed profile, while for large values of P_max the morphology of the bed is significantly influenced by the presence of surface biofilm. In the case of small (Fig. 6b) and medium (Fig. 6c) values of P_max the final bed profiles are similar, even if the smaller growth rate parameter results in a slightly higher bed mobility and the bed reaches the final bed equilibrium condition sooner. A large maximum growth rate parameter influences the morphological evolution of the channel by promoting the development of surface biofilm from the early stages of the simulation and reducing the sediment mobility (Fig. 6a, green line). Under this condition, the bed after 30 000 tidal cycles is still dynamic at both the landward and seaward boundary. Overall, the grain size distribution of the channel bed is preferentially coarse in the seaward boundary and fine at the landward boundary. Simulations with small values of the growth parameter result in higher sediment mobility, and the deposit at the landward side is relatively coarse (Fig. 6e) compared with the stratigraphy of the deposit created in the case with a large growth rate parameter. In this case, due to the high stabilisation, the coarse fraction characterises the bed surface at the sea boundary (Fig. 6g).

Figure 6 Effect of different values of maximum growth rate (P_max) on the surface biofilm biomass (a). The top row represents the bed evolution profile with small (b), medium (c) and large (d) P_max after 30 000 simulated tidal cycles. The bottom row represents the geometric mean diameter of the final deposit in the case of small (e), medium (f) and large (g) P_max. The bed profiles are compared with the final bed elevation of clean sediment at equilibrium (dashed red lines). The dashed blue and green lines represent the initial water surface elevation and the initial bed profile, respectively. The initial geometric mean size of the bed is 0.30 mm.

3.4 Effect of temperature variation

Biofilm growth differs during the course of a year due to environmental conditions, with a higher growth rate during spring and beginning of summer (Thom et al., 2015; Widdows et al., 2000). The seasonality and the variation in the sediment temperature affect the development of biofilm and the consequent morphological evolution of the channel, as shown in Fig. 7a. The simulations presented here focus on the effect of seasonality and changes in temperature; therefore, the effect of carpet-like erosion is neglected.

The variation in the sediment temperature is function of the light availability and the turbidity of the water column. Here a sediment temperature variation of ±5 ^∘C is assumed compared to the previously simulated temperature profile, to simulate the possible scenarios in shallow-water environments (Pivato et al., 2019). As mentioned before, the amount of biofilm biomass developed on the bed surface is strongly regulated by the sediment temperature (Fig. 7a). Compared with the reference case (orange line), a decrease in the annual sediment temperature (T −5 ^∘C) results in an overall slower development of biofilm; in other words, it takes longer for the biofilm to reach the maximum amount of biofilm biomass at the bed (approximately 150 simulated days, blue line Fig. 7a). In this scenario the sediment temperature does not reach the maximum temperature for photosynthesis (T_max), resulting in a stable biofilm biomass over a relatively long period (∼ between 150 and 250 d). While an increase in the annual sediment temperature (T+5 ^∘C) would result in a more rapid development of biofilm compared with the reference case, reaching the maximum amount of surface biomass after approximately 80 simulated days (green line in Fig. 7a). The sediment temperature reaches and surpasses the maximum temperature for photosynthesis (T_max), resulting in a decrease in surface biofilm.

In the reference scenario (Fig. 7a, orange line), the total amount of biofilm biomass covering the bed over the year interval is comparable to the case of low sediment temperature (Fig. 7a, blue line), but these two scenarios result in a slightly different final bed profile. In the case of low sediment temperature, the bed is covered by biofilm for the period between 120 and 240 d (May–August), resulting in a more mobile bed (Fig. 7b). In the reference case the bed shows the presence of biofilm for a longer period of time, even if it is not always at its maximum value (between 90 and 280 d), and the bed evolves more slowly (Fig. 7c). An increase in temperature (Fig. 7a, green line) results in high biofilm biomass around day 70 and 300, while for the rest of the year the presence of biomass on the bed is low; therefore, the bed results are more mobile (Fig. 7d).

The stratigraphy of the deposit emplaced is coarser at the ocean boundary and finer at the landward boundary after reflection for the reference case (Fig. 7f) compared to the other simulated temperature conditions.

Figure 7 Effect of sediment temperature on the development of biofilm (a). Bed evolution profile (b, c, d) and final stratigraphy of the deposit (e, f, g) after 30 000 tidal cycles in the case of low sediment temperature (−5 ^∘C, left panels), the reference case (centre panels) and high sediment temperature (+5 ^∘C, right panels). The bed profiles are compared with the clean sediment final bed elevation at equilibrium (dashed red lines). The dashed blue and green lines represent the initial water surface elevation and the initial bed profile, respectively. The initial geometric mean size of the bed is 0.30 mm.

3.5 Effect of the sediment bio-cohesivity parameter (α)

There is no universal relationship available in the literature between critical shear stress for erosion (τ_bc) and the amount of Chl-a, which is considered an approximation of biostabilisation potential. This uncertainty can be explained by sediment rheology and different sampling techniques. Furthermore, the distribution of Chl-a content can vary spatially due to the small-scale morphology of the bed (Le Hir et al., 2007).

The effect of this variability has been investigated by changing the parameter (α) used to correlate the critical shear stress for erosion (τ_bc) with the amount of Chl-a on the bed (Eq. 2). The results of the channel morphology and stratigraphy obtained by assuming α is equal to 0.01, as suggested in the literature (Le Hir et al., 2007; Mariotti and Fagherazzi, 2012), are compared with scenario that account for the variability of the sediment bio-cohesivity (Fig. 8). The values of the bio-cohesivity parameters tested here (α= 0.001 and 0.02) have been suggested by previous studies (Le Hir et al., 2007).

Small values of the bio-cohesivity parameter (α= 0.001) result in higher channel mobility, which is able to reach equilibrium by the end of the simulation (Fig. 8a). A further increase in the bio-cohesivity parameter results in a slower morphological evolution of the channel. After 30 000 simulated tidal cycles, the channel is still evolving in time (Fig. 8c).

For small values of the bio-cohesivity parameter, the final surface grain size distribution (Fig. 8d) is mostly characterised by fine sediment (D < 0.25 mm), with the coarse fraction covering the landward boundary of the domain (x < 10 m), which differs from the the other two cases where a larger surface area is increasingly covered by coarse material, inhibiting sediment motion (Fig. 8e and f).

Figure 8 Effect of the bio-cohesivity parameter (α) that relates the biofilm with critical shear stress for erosion. Bed evolution profile (a, b, c) and final stratigraphy of the deposit after 30 000 tidal cycles (d, e, f) in the case of α= 0.001 (left), the reference case (α= 0.01, centre) and α= 0.02 (right). The bed profiles are compared with the clean sediment final bed elevation at equilibrium (dashed red lines). The dashed blue and green lines represent the initial water surface elevation and the initial bed profile, respectively. The initial geometric mean size of the bed is 0.30 mm.

4 Discussion

The complex interaction between physical, chemical and biological processes and properties that govern sediment transport mechanisms are still poorly understood and quantified. It is therefore difficult for morphodynamic models to be accurate and predict into the future. Whilst some factors will be similar between estuaries, our findings confirm the need for site-specific calibration of morphodynamic models. These models must account for the contribution of different eco-engineers to tidal flat development. Nonetheless, our investigation offers both fundamental qualitative and quantitative information regarding the role of key environmental parameters in sediment stability and morphological evolution in a simplified intertidal channel.

Local hydrodynamic conditions (e.g. tides, waves) affect not only the establishment of biofilms but also their recovery processes (Defew et al., 2002). Small but frequent disturbances hinder the early stages of biofilm development, while strong disturbances can detach established biofilm (Fig. 4). It is reasonable to conclude that local hydrodynamics play a crucial role in mediating the presence of biofilm, with carpet-like being erosion possible when disturbance is high. Results presented in this study show that in low dynamic environments where carpet-like erosion is not dominant (e.g. on bars, in central areas of tidal flats), biofilm growth is prominent (inserts of Fig. 5c), resulting in a strong biostabilising effect on the bed (Fig. 5c). This supports field investigations where higher bed stability is observed in central tidal flats compared to the edges (Widdows at al., 2000) or there are differences between channels and flats (Daggers et al., 2020). Biofilm presence inhibited the sediment movement for all shallow-water habitats with low tidal forces, as demonstrated by a lack of significant changes after 30 000 tidal cycles (Fig. 5c, d, f, g). Furthermore, deposited sediment was coarse (Fig. 3Ac, d, f, g). In high dynamic environments, carpet-like erosion can remove surface biofilm exposing the clean sediment underneath and reducing biostabilisation (e.g. close to the channel, at the edge of the tidal flat), resulting in a more mobile bed profile (Fig. 5i, j, k, l). Moreover, high bed shear stresses due to hydrodynamic forces (tides) can cause a general delay in biofilm formation and biostabilisation (Fig. 4) and a significant decrease in the biofilm stability (Schmidt et al., 2018). This study does not incorporate the combined hydrodynamic effect on surface biofilm mass of occasional storms and periodic tidal forces. Morphology and sedimentary processes on tidal flats can be strongly affected by storms and associated high-energy activities over a short time. The simulations presented in Fig. 4 assume periodic disturbances to investigate the effect of changes in biofilm model parameters on biofilm establishment and growth. The frequency and intensity of storms is likely to increase in the future due to climate change, and the resulting drastic morphological changes on tidal flats can occur over short durations. This will affect biofilm evolution and establishment and therefore the degree of biological stabilisation that occurs. Storms can induce strong wave activities, elevate water levels and cause severe erosion of tidal flats due to enhanced bed shear stress and carpet-like erosion of surface biofilm. The associated high suspended sediment concentrations and long inundation period increase the turbidity of the water column and inhibit photosynthesis. The model presented here can be easily adapted to account for the seasonal variability of storms by incorporating the combined hydrodynamic effects of occasional storms and periodic tidal forces.

Simulations presented here demonstrate that the biostabilising effect due to the presence of biofilm decreases the time needed for the bed to reach equilibrium compared to clean, abiotic sediment (Fig. 5a). The deposits are finer than the initial bed condition at the landward boundary, which is particularly relevant as physically cohesive sediment like mud facilitates salt marsh survival and MPB growth by supplying nutrients to the bed (Smith and Underwood, 2000; Underwood, 2002). The temperature model implemented in this study in turn promotes further sedimentation and can limit mud erosion (Brückner et al., 2020), which is fundamental for the stabilisation of wider estuarine morphological features, bank accretion and stability, predicting estuarine and deltaic development, and coastal protection. Consequently, an increasing extent and thickness of mud cover might lead to a stabilisation of large-scale estuarine morphology. Although not directly modelled in this study, our findings suggest that the sediment bed would become “muddier” as biostabilisation is increased, and these changes may influence wider estuarine morphology as channels are stabilised, attract more mud and influence the evolution of channel morphology. Even when trends are observed between the amount of benthic biofilm and the grain size distribution at the bed, the relationship between these two parameters is not straightforward. These relationships are strongly modulated by the role played by a complex interaction of other factors, such as the light reaching the bottom; nutrient fluxes; human activity; and the community composition of the primary producers present such as diatoms, cyanobacteria, and green algae (Cahoon et al., 1999; Schmidt et al., 2018). Furthermore, at energetic and sandy sites, the frequent reworking of the substrate results in removal of the biofilm and more mobile bed, while under less dynamic conditions, even small increases in the fine and muddy sediment fraction can promote sediment stability (Hope et al., 2020).

In aquatic environments, spatial variability in water temperature can be natural (e.g. geothermal activity, source of water) or it can result from direct changes in local land use and activities (e.g. deforestation, industrial activities) or indirect and global changes (e.g. climate change) (Caissie, 2006; van Vliet et al., 2011). Alteration of thermal regimes can be a major determinant of changes in the diversity and resilience of aquatic biota from primary producers to consumers. The temperature model implemented in this study used surface sediment temperature as a key parameter for the growth of biofilm (Pivato et al., 2019, 2018) in temperate areas (orange line in Fig. 7a) and tested for ±5 ^∘C variations in temperature from the reference case. Favourable temperature conditions would result in changes in biofilm biomass production and affect fine sediment dynamics via local stabilisation and accretion, reduce the turbidity in the water column, and change the hydrodynamic conditions (reduce the bed roughness). In summary, the goal of these simulations was to investigate bed morphology in colder and warmer climates as reported in the literature. In the first case the activity of MPB is restricted to the warmer periods (blue line in Fig. 7a), while in the second case biological biomass at the bed can develop more during cooler months (green line in Fig. 7a, Hope et al., 2019). Temperature-regulated biofilm development at the bed strongly influences the final morphology of the channel. Dissolved oxygen levels are directly linked to water temperature, with low values of saturated dissolved oxygen for higher water temperature (Pivato et al., 2019). Projected future temperature increases could thus lead to a decrease in ocean oxygen solubility and have a direct effect on organismal physiology and on biofilm development, especially in shallow-water basins located in temperate regions (Kent et al., 2018). Nutrient cycling and carbon flows through benthic communities are influenced by chemical and biological processes, which are regulated by sediment temperature and light availability. Therefore, by considering the effect of biostabilisation here, this model indirectly accounts for the effect of water and sediment temperature on the morphodynamic evolution of coastal shallow bays (Marani et al., 2007, 2010; Mariotti and Fagherazzi, 2012).

Biofilm growth rate and seasonality are key parameters when modelling biostabilisation (Figs. 5 and 6). Large variation in biostabilisation between seasons is reported in the literature, with the highest values in spring and the lowest in late autumn (Underwood and Paterson, 1993; Marcarelli et al., 2008; Thom et al., 2015; Schmidt et al., 2016; Waqas et al., 2020). This is due to the differences in biofilm growth and composition resulting in mechanically diverse responses to the increased bed shear stress. Experiments conducted by Thom et al. (2015) reported a 10-fold increase in sediment stability, depending on boundary conditions and investigated season, and the hydrodynamic erosional process can also be influenced by seasonality, highlighting the heterogeneity of the process. Biostabilisation is considerably higher in spring than in summer, supported by the fact that EPS protein and carbohydrate contents increase (Amos et al., 2003; Dickhudt et al., 2009; Thom et al., 2015; Schmidt et al., 2016, 2018). Seasonality also affects bed morphology from early spring until the onset of summer, with 80 % of the surface of the intertidal flats covered in biofilm, which can enhance the formation of a hummock–hollow pattern (Weerman et al., 2011a). This trend is observed in temperate humid climates (cold winters and mild summers, Figure 7a blue line, Fig. 7b and e), while in warm temperate climates (mild winters and hot summers, Fig. 7a green line, Fig. 7d and g) the seasonal MPB biomass maximum is most likely to occur in late autumn (Haro et al., 2022). Friend et al. (2003) also observed strong seasonally dependent relationships between critical shear stress for erosion (τ_c), habitat type, Chl-a and bed elevation; in fact the seasonal activity of the species contributes significantly to increasing or decreasing the sediment stability (Thom et al., 2015). This aspect has been parameterised in this study in the maximum growth rate parameter (P_max), which accounts for the effect of seasonality according to a sediment temperature model (data available in the literature assume that this parameter ranges between 0.0078 and 1.10 d⁻¹) (Labiod et al., 2007; Mariotti and Fagherazzi, 2012; Uehlinger et al., 1996).

Experimental and field studies have attempted to identify the roles of biological and physical processes in sediment stability using regression analyses to relate the erosion threshold to biological and physical parameters (Defew et al., 2003; Amos et al., 2004; Droppo et al., 2007; Grabowski et al., 2011). The presence of biofilm can increase bed stability up to 500 % compared with non-colonised sediment (Le Hir et al., 2007; Zhu et al., 2019), and the effect of EPS is much greater than physical cohesion (Malarkey et al., 2015; Parsons et al., 2016). As Chl-a has a strong functional relationship with stabilising EPS (Friend et al., 2003; Paterson et al., 1994; Underwood et al., 1995), Chl-a is often regressed against erodibility (Le Hir et al., 2007). However, there remains no universal relationship available in the literature for the sediment bio-cohesivity parameter (α) that correlates the critical shear stress for erosion and Chl-a, besides the observation that the critical shear stress for erosion increases as the Chl-a content increases (Le Hir et al., 2007). Sample techniques, timing and in situ versus laboratory measurement limitations play an important role in the variable relationships observed between stability and biological cohesion (Tolhurst et al., 2000a, b; Le Hir et al., 2007; Perkins et al., 2003; Widdows et al., 2007).

While remote sensing techniques can help capture spatial and temporal variability in MPB dynamics (Méléder et al., 2020; Haro et al., 2022; Hakvoort et al., 1998, Méléder et al., 2020, Paterson et al., 1998), which is often missing from in situ investigations, modelling studies, such as the current study, can help to elucidate long-term changes to the bed morphology. Fully numerically reproducing the influence of biological cohesion in different habitats is still a challenge due to the complexity of intertidal systems, and a more detailed parameterisation of MPB effects is required to properly describe these complex environments (Hope et al., 2020). Significant knowledge gaps remain relating to how small-scale biological activity can impact large-scale cohesive sediment dynamics and overall landscape evolution. The results herein demonstrate that biophysical scale-dependent feedbacks are crucial in regulating the substrate and the spatial self-organisation of intertidal ecosystems. This process is fundamental not only for understanding the development of present channels but also for dating ancient biogenic deposits (Van de Vijsel et al., 2020) and ancient biostabilisation processes (e.g. microbialites; Burne and Moore, 1987; stromatolites; Hohl and Viehmann, 2021; Noffke et al., 2013). Further development of the model is required to account for the long-term effect of sticky microbial biofilms on the substrate and their effect on the landscape development.

While this study provides a sensitivity analysis of the biofilm model parameters, several assumptions and simplifications relating to the complexity of the biogeomorphology of these environments have been made. Resuspension of MPB in the water column in highly productive ecosystems will promote the settlement of MPB and the establishment of surface biofilms. In intertidal environments, MPB, macrophytes and fauna are heterogeneously distributed, i.e. found in patches, though typically on rather small scales. While biofilm patchiness can influence the grain size of the bed, this study assumes a uniformly distributed biofilm. Similarly, macrofaunal activity can increase the bottom roughness and surface heterogeneity (Borsje et al., 2009; Coco et al., 2006; Brückner et al., 2021), enhancing the complexity of the interactions between stabilising microbes and macrofauna. For example, infauna excrete essential nutrients that stimulate the growth of MPB and therefore bed stability (Murray et al., 2014). Further, their burrows, mounds and tube mats increase the surface area of sediment, creating a patchy distribution of nutrients on the substrate and enhancing the spatial complexity of biofilm distribution. These positive effects on MPB can negate the destabilising effects of burrowing or grazing (Hope et al., 2019). The presence of grazers and the abundance of nutrients can work differently at different spatial and temporal scales, and this often creates complex interactions that are difficult to quantify (Posey et al., 1999). On intertidal flats, spatial self-organisation of microbes observed during the early spring months can be destroyed as the season progresses. This shift towards a more homogenous surface is attributed to the presence of herbivores, bioturbation activity and the increase in grazing activity as the season progresses (Weerman et al., 2010, 2011a, b).

5 Conclusion

The study presented here has provided a novel insight into the biomorphodynamic evolution of intertidal channels. Biofilm effects and the influence of seasonality and temperature changes on biostabilisation potential were included. The 1D biostabilisation shallow-water model was implemented under different hydrodynamic conditions to investigate different climate scenarios and identify biofilm development parameters that influence the final channel morphology.

The model can be utilised to investigate the bed and deposit evolution in tidal-dominated channels, starting from a horizontal bed until it reaches equilibrium. The output suggests that high hydrodynamic disturbances play a fundamental role in shaping the channel equilibrium profile by creating carpet-like erosion of the biofilm layer, which exposes the clean sediment underneath. Low hydrodynamic forces (e.g. supratidal area) allow the steady development of biofilm, and the consequent biostabilisation can inhibit sediment mobility. The frequency and intensity of the hydrodynamic disturbances therefore regulates the growth and stability of the biofilm.

Changes in the annual sediment temperature profile (e.g. due to climate change) or in the biofilm maximum growth factor (regulated by, e.g. nutrient availability) strongly influence the amount of surface biofilm and as a consequence the bed profile and stratigraphy. Increasing and decreasing the sediment temperature from the optimum for photosynthesis both result in a less stable and less developed biofilm, and as a consequence the bed is more mobile.

It is concluded that hydrodynamic forces also play a decisive role in shaping the geometry of the channel in the uniform presence of surface biofilm, but the stratigraphy of the deposit is significantly affected by the biofilm conditions.

Appendix A: Model validation

The one-dimensional shallow-water equations modified for partially dry areas are solved simultaneously using the explicit, second-order-accurate in space and time predictor–corrector MacCormack scheme (Chaudhry, 2008; Viparelli et al., 2019). The numerical model is implemented on tide-dominated horizontal channels subject to tidal fluctuations at the ocean boundary, which results in erosion in the ocean part and a landward-migrating shoal, depositing and forming a beach until it reaches equilibrium conditions. The domain is divided into N-cells of width Δx, set equal to 0.5 m to have enough spatial resolution. The bed and water surface elevation with respect to the data are denoted by η(i) and ξ(i), respectively.

An impermeable wall is assumed at the landward boundary ($Q{|}_{x=\mathrm{0}}=$ 0, $Q_{\mathrm{b}}{|}_{x=\mathrm{0}}=$ 0). An open-ocean or tidal basin is assumed at the ocean boundary (x= L) with amplitude α_t and periodicity ω_t, from where tides propagate into the domain:

$$\begin{array}{}\text{(A1)}& {\mathit{\xi }}_{\mathrm{d}}={\mathit{\xi }}_{\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{t}}\cos (\mathrm{2}\mathit{\pi }t/{\mathit{\omega }}_{\mathrm{t}}).\end{array}$$

Extra points are added at the land and ocean boundaries of the domain to compute the predictor and corrector terms, respectively; zero gradient for discharge and water surface elevation is assumed at the land boundary (x= 0), while at the ocean boundary the flow rate and the water surface elevation are set equal to the value at (x= L) (Viparelli et al., 2019).

The final numerically modelled bed profile after 2000 tidal cycles shows good agreement with the temporal evolution of the cross-sectional averaged bed profile Tambroni et al. (2005) obtained from laboratory investigation of the process, whereby an equilibrium morphology is established in a tidal system consisting of an erodible channel connected through an inlet to a tidal sea (Fig. A1). The bed profile generated from the numerical model shows weaker concavity, resulting in better match with the theoretical predictions suggested by Seminara et al. (2010). Seminara et al. (2010) proposed two theoretical predictions for tidal-dominated channels, assuming the Chézy coefficient to be constant (C_constant= 12) or a function of the outer bottom profile at equilibrium (D₀; $C_{\mathrm{variable}}=C_{\mathrm{0}}{D}_{\mathrm{0}}^{\mathrm{1}/\mathrm{6}}$) (dashed lines in Fig. A1). The numerically simulated channel slightly underestimates the bed elevation at the entrance at the landward boundary (Fig. A1).

Figure A1 The experimental bed profile (grey symbols) observed by Tambroni et al. (2005) after 2000 tidal cycles in a straight, tidal channel with constant width and the theoretical predictions (two dashed lines) resulting from equations suggested by Seminara et al. (2010), computed with a constant and variable Chézy flow conductance, are compared with the modelled bed profile (red line).

Figure A2 Logarithmic root-mean-square error (RMSE) from the comparison between the model run with different numbers of computational cells (x axes) and the analytical solution computed with a variable Chézy flow conductance (Seminara et al., 2010).

Grid sensitivity analysis has been performed by investigating different range of computational grid points in the streamwise direction. Increasing the grid resolution did not show any significant effect on the results (Fig. A2).

A1 Stratigraphy of the final deposit

The sediment grain size distribution of the deposit after 30 000 simulated tidal cycles is presented here.

For the reference case of clean sediment, the grain size distribution of the landward deposit associated with shoal reflection coarsened in the upward direction and from the ocean to the land (Fig. A3a). Coarse sediment is transported upstream of the shoal and is deposited in the landward part of the channel forming the coarse basal part of the deposit. As the shoal approached the landward boundary, fine sediment is deposited on the basal layer. Sediment deposited after the shoal reflection presented a fining upward profile for decreasing velocities associated with beach formation.

First, the case of spatially uniform and stable biofilm on the bed surface for the entire duration of the simulation (Fig. A3b) is modelled. The sediment mean diameter at the bed surface at the ocean boundary is coarser than the initial condition (Fig. A3c). Figure A3d shows the mean diameter of the deposit emplaced with the water depth as a constraint for the development of biofilm (H < 0.05 m, see Sect. 3.2). The mean sediment diameter at the bed surface results are coarser than the initial condition (0.35 mm).

Considering the effect of seasonality (Fig. A3e) in the case of uniformly distributed surface biofilm, coarse sediment is found on the seabed (Fig. A3f). In the case of water depth constraining the development of biofilm, coarse sediment is found at the ocean boundary, while landward the bed is characterised by finer deposit (Fig. A3g).

Figure A3i and j show the evolution of the bed with surface biofilm periodically removed by the tidal-induced stresses on the bed (carpet-like erosion, Fig. A3h). In both examples, spatial distributions constrain conditions for surface biofilm. The stratigraphy of the deposit is analogous to the clear sediment case (Fig. A3a), with initial coarse sediment deposited landward until the shore reflects creating a lens of fine material, and after that point more coarse sediment is deposited. Assuming that biofilm uniformly develops spatially, the deposit emplaced after shoal reflection at the landward boundary is coarse (Fig. A3i), and the sediment at the bed is finer overall compared with the case of biofilm that only develops in shallow-water areas (Fig. A3j). Analogous observations can be made for the pattern of the final grain size distribution of the deposit in the case of a combined effect of seasonality and carpet-like erosion.

Figure A3 Final stratigraphy of the deposit after 30 000 simulated tidal cycles under different biofilm conditions. Panel (a) shows the stratigraphy of the deposit in the reference case with clean sediment. The panels in the left column represent the time evolution of biofilm during the simulation, showing stable biofilm (b) and the effect of seasonality (e) and carpet-like erosion (h). The rows represent, respectively, the stratigraphy of the deposit under these biofilm conditions. The panels in the central column show the stratigraphy of the deposit in the case of biofilm that is uniformly distributed spatially. The panels in the right column show the stratigraphy of the deposit in the case of biofilm that has developed only in locations where the water depth is smaller than 0.05 m. The initial geometric mean size of the sediment at the bed is 0.30 mm.

Appendix B: List of symbols

A	Water surface albedo
ar	Characteristic length scale of the bed irregularities
Ai	Cross-sectional area averaged over bed irregularities
	(Wξ)
B	Biofilm biomass
Bmin	Background biofilm
c0i	Near-bed concentration of suspended sediment in the
	generic grain size range averaged over turbulence
Cf	Friction coefficient
ci	Volumetric sediment concentration
D50	Median diameter of the bed material
Dg	Geometric mean sediment grain size
Di	Characteristic diameter
E	Catastrophic erosion
Ei	Grain-size-specific entrainment rate under equilibrium
	of suspension
Ek	Light saturation parameter
ET	Entrainment rate per unit bed summed over all
	grain sizes
Ei	Entrainment rate per unit bed for each grain size i
FH	Wet fraction of the channel bed
g	Acceleration of gravity
H	Water depth
Hres	Light availability
KB	Half-saturation constant for biofilm growth
L	Channel length
N	Number of computational nodes
pi	Fraction of sediment in each grain size range
PB	Effective maximum growth rate for biofilm
Pmax	Maximum growth rate for biofilm
$P_{\max }^{\mathrm{B}}$	Biofilm growth rate under light saturation conditions
Q	Flow discharge (AcU)
Qb	Total material load as the sum of the contribution of
	bed load and suspended load summed over all
	grain sizes
Qb,bi	Total volumetric bed material load as the contribution of
	bed load for the generic grain size i
Qb,si	Total volumetric bed material load as the contribution of
	suspended load for the generic grain size i
R	Submerged specific gravity of the bed material
RH	Hydraulic radius ($A_{\mathrm{c}}/\mathit{\chi }$)
Rsun	Solar irradiance reaching the water surface
Sf	Friction slope
t	Temporal coordinate
ti	Time detachment due to high hydrodynamic forces
τb	Average bed shear stress
Ts0	Surface sediment temperature
Tmax	Maximum temperature for photosynthesis
Topt	Optimal temperature for photosynthesis
U	Flow velocity
u*c	Critical shear velocity
u*s	Shear velocity due to skin friction
vsi	Fall velocity in each grain size range
W	Channel width
X	Longitudinal coordinate
Y	Effective flow depth

α	Bio-cohesivity parameter
αt	Tidal amplitude
β	Shape parameter
ε	Global decay
δ	Dirac function
η	Bed profile
ηo	Initial bed elevation
λ	Extinction coefficient
ξ	Mean water surface elevation
ρs	Density of the sediment
τb	Bed shear stress
τbc	Critical shear stress for erosion
τbc,0	Clean sediment critical shear stress
τbs	Bed shear stress due to skin friction
${\mathit{\tau }}_i^{*}$	Grain-size-specific Shields number
${\mathit{\tau }}_{\mathrm{bc},i}^{*}$	Grain-size-specific reference Shields
	number for significant bed load transport
χ	Wetted perimeter
ωt	Tidal period
$\mathit{\zeta }=z/b$	Dimensionless upward normal coordinate