## Abstract

The presence and continued existence of tidal morphologies, and in particular of salt marshes, is intimately connected with biological activity, especially with the presence of halophytic vegetation. Here, we review recent contributions to tidal biogeomorphology and identify the presence of multiple competing stable states arising from a two-way feedback between biomass productivity and topographic elevation. Hence, through the analysis of previous and new results on spatially extended biogeomorphological systems, we show that multiple stable states constitute a unifying framework explaining emerging patterns in tidal environments from the local to the system scale. Furthermore, in contrast with traditional views we propose that biota in tidal environments is not just passively adapting to morphological features prescribed by sediment transport, but rather it is ‘The Secret Gardener’, fundamentally constructing the tidal landscape. The proposed framework allows to identify the observable signature of the biogeomorphic feedbacks underlying tidal landscapes and to explore the response and resilience of tidal biogeomorphic patterns to variations in the forcings, such as the rate of relative sea-level rise.

## 1. Introduction

Tidal landforms in lagoons and estuaries exist in a constant pursuit of sea-level rise, offset by an interplay of inorganic and organic sediment deposition (e.g. [1–4]). An outcome of such pursuit is the generation of some striking biological and morphological patterns at different scales. At the large scale, marshes, tidal flats and subtidal areas form the ‘universal’ texture of the tidal zone [3–7]. At a smaller scale, zonation patterns, a mosaic of sharply bounded and nearly homogeneous vegetation patches, are widespread in marshes worldwide [8–10]. Indeed, many aspects of tidal biogeomorphodynamics and of the associated patterns have been seen to be tightly linked to interactions between biotic and abiotic processes [11–18].

We point out here that intertidal vegetation and microalgae act as ecosystem engineers [19] by exerting a control on their substrate elevation, by stabilizing the sediment against erosion and by directly producing organic sediment. We also suggest that ecosystem engineering plays a fundamental role in determining tidal biogeomorphic patterns, in contrast with their traditional interpretation as the result of a passive adaptation of biota to controlling physical processes and to the associated edaphic conditions [1]. In particular, we show that depositional intertidal patterns emerge from, and are the signature of, two-way feedbacks between accretion rates, inorganic sediment transport and biological processes controlling organic sediment productivity. To this end, we first briefly review contributions on biogeomorphic equilibria in a conceptual lumped setting [20,21]. Subsequently, building on new observational and modelling results [22], we introduce original analyses on the spatially extended equilibrium configuration of salt marshes, showing how multiple competing stable states, largely regulated by biota, provide a unifying framework to explain emerging patterns in tidal environments. Through this framework, we finally explore how tidal biogeomorphic patterns change with vegetation characteristics and competition dynamics, as well as their responses and resilience to variations in the forcings, such as rates of relative sea-level rise.

## 2. Depositional tidal biogeomorphology

The time evolution of topographic elevation *z*(**x**,*t*) (computed with respect to mean sea level, hereinafter MSL) at position **x** (in a planar cartesian reference frame) is jointly determined by the rate (e.g. in (mm yr^{−1})) of inorganic sediment deposition, *Q*_{s}(*B*,**x**,*t*) (the sum of sediment settling owing to gravity and trapping by standing biomass), the rate of organic soil production *Q*_{o}(*B*,**x**,*t*) (associated chiefly with below-ground biomass production), the erosion rate *E*(*B*,MPB,**x**,*t*) (owing to wind waves and, possibly, to tidal flows) and the rate of relative sea-level rise, *R* (herafter, RSLR, the algebraic sum of sea-level rise and local subsidence, possibly a function of time and space):
2.1MPB, accounts for the stabilizing effect of MicroPhytoBenthos, the benthic microalgae whose polymeric excretion forms a biofilm with high mechanical strength [3,23]. *B*(**x**,*t*) is the biomass produced by halophytic vegetation, i.e. vegetation adapted to living in the hypersaline and hypoxic marsh soils [10,24], which affects inorganic deposition (by increasing flow energy dissipation, reducing the turbulent kinetic energy and by direct trapping, e.g. [25–27]), produces organic soil (e.g. [4,28,29]) and reduces erosion (by damping wind waves and stabilizing the sediment, e.g. [30–32]). In turn, vegetation development is strongly linked to marsh elevation (summarizing the effects of sediment salinity and aeration, the stress factors that are more directly acting on marsh vegetation), such that a strong feedback exists between geomorphic characteristics and biomass production. Equation (2.1) inherently refers to a ‘vertical’ accounting of deposition and erosion contributions, and hence does not include edge erosion processes. We note, however, that neglecting such processes does not affect the validity of the results presented here, as ‘vertical’ dynamics and horizontal erosion of tidal landforms are largely independent [33,34].

Because temporal variations in the bottom profile occur over time scales typical of morphological changes, i.e. time scales of one to several years, equation (2.1) can be regarded as the annually averaged mass balance equation. *Q*_{s}(*B*,**x**,*t*), *Q*_{o}(*B*,**x**,*t*) and *E*(*B*,MPB,**x**,*t*), however, must be obtained by first resolving the relevant biotic and abiotic processes at the time scale of a single tidal cycle and by then integrating the corresponding sediment fluxes to the yearly scale. Steady-state solutions to equation (2.1) will be described below for a point model, in which the landform being analysed is ‘collapsed’ into a single point and *z* can be interpreted as its mean elevation, and for a one-dimensional representation, which allows the description of the smaller scale morphological properties.

## 3. A point model: large-scale tidal patterns

### (a) Modelling framework

Marani *et al.* [3,20] and D’Alpaos *et al.* [21] have examined the numerical solution of equation (2.1) when this is integrated in space to describe the time evolution of the mean elevation of a tidal platform (scales indicatively larger than 1 km). Hydrodynamic flows and sediment transport in this lumped approach are parametrized to describe the exchange of water and sediment between the channel (where the tidal wave propagates and suspended sediments are transported) and the tidal platform. This is done according to the following conservation equation [3,35]:
3.1where is the instantaneous tidal elevation with respect to the local MSL (*H* being the tidal amplitude and *T*=12 h the main tidal period). *D*(*t*)=*h*(*t*)−*z* is the instantaneous water depth. *C* is the instantaneous suspended sediment concentration (SSC), when d*h*/d*t*>0 (*C*_{0} being the prescribed SSC in the channel), whereas when d*h*/d*t*<0.

Note that *w*_{s} *C*(*t*)/*ρ*_{b} and *αB*^{β}*C*(*t*)/*ρ*_{b} are the instantaneous settling and trapping fluxes, respectively, which, by integration over all tidal cycles in 1 year yield the yearly mean total inorganic flux, *Q*_{s}(*t*), from the water column to the surface [20]. We assume *w*_{s}=0.2 mm s^{−1} for the settling velocity (typical of a sediment size of 20 μm [36]), whereas *α* and *β* synthetically account for suspended sediment trapping by the canopy (see [20] for details). *ρ*_{b}=*ρ*_{s}⋅(1−*p*) is the bulk density of the sediment after compaction has taken place, with *ρ*_{s}=2650 kg m^{−3} and porosity *p*=0.5.

Erosion is assumed to be zero when vegetation is present [30,31], otherwise it is *E*(*B*,MPB,*t*)∝(*τ*(*z*,*t*)−*τ*_{c}(MPB)), where *τ*_{c}(MPB) is the threshold shear stress above which erosion starts (*E*(*B*,MPB,*t*)=0 if *τ*(*z*,*t*)<*τ*_{c}(MPB)), and *τ*(*z*,*t*) is the wave-induced bottom shear stress (a function of the local wind ‘climate’, and of time-varying water depths). *τ*_{c}(MPB) is a function of the possible presence of microphytobenthos.

Vegetation may be assumed to respond over time scales (1 to a few years) which are shorter than the characteristic response times of geomorphological processes (usually of the order of several to tens of years [20]). Hence, rather than introducing a full dynamical equation for the time evolution of biomass, we assume it to be in equilibrium with the current landform elevation, i.e. we simply prescribe it as *B*(*z*). From observations in Mediterranean marshes, where several vegetation species compete for marsh space, one can assume that there a biomass production which is an increasing function of elevation (for marsh platform elevations between MSL (*z*=0) and *z*=*H*). We call this the *multi-species* case: *B*(*z*)=*B*_{0}⋅2⋅*z*/(*H*+*z*) (considering *B*_{0}≅1 kg m^{−2}, *B*(*z*)=*B*_{0} if *z*>*H* and *B*(*z*)=0 if *z*<0 [3]). Field observations in *Spartina-dominated* areas [4,37] give, on the other hand, a biomass productivity which decreases with *z* in the range 0<*z*<*H*: *B*(*z*)=*B*_{0}⋅2⋅(*H*−*z*)/(2*H*−*z*); *B*(*z*)=0 if *z*>*H* or *z*<0. In either case the organic contribution to accretion is expressed as *Q*_{o}(*z*)=*γ*⋅*B*(*z*), where *γ*=2.5×10^{−6} m^{−3} yr^{−1} g^{−1} incorporates typical vegetation characteristics and the density (after compaction and partial decomposition) of the organic soil produced [20,29].

### (b) Stable equilibria and corresponding biogeomorphic patterns

Use of the above mathematical representations of the controlling physical processes makes the deposition and erosion fluxes on the right-hand side of equation (2.1) dependent just on *z*, such that one obtains an equation describing a one-dimensional dynamical system: d*z*/d*t*=*F*(*z*)−*R*. *R* plays the role of an external forcing and *F*(*z*) is a nonlinear function as a result of trade-offs between biotic and abiotic geomorphic processes. The nonlinearity of *F*(*z*) in itself implies that the system will display a set of stable and unstable fixed points, corresponding to an equal number of competing equilibrium states [38]. Hence, from these considerations alone, we can infer the presence of a path to pattern formation: only specific and recurring landforms may exist within a tidal environment, corresponding to the available equilibrium states as mediated by the type of vegetation (*Spartina-dominated* or *multi-species*) and by space-varying forcings (e.g. RSLR and available SSC).

If the rate of RSLR is constant in space and time, it is possible to easily identify the accessible equilibria by geometrically identifying the values of *z* corresponding to the condition d*z*/d*t*=0. Stable fixed points are defined by d/d*z*(d*z*/d*t*)<0, such that, depending on *R* and *C*_{0}, up to three stable equilibria may occur on the three descending branches of *F*(*z*) illustrated in figure 1 (for the Spartina-dominated case). Such stable equilibria correspond to the basic tidal landforms: (i) subtidal platforms, with elevation lower than the minimum tidal level ( m in this case); (ii) tidal flats, with elevation between and MSL; and (iii) vegetated marshes (colonized by *Spartina* in the case of figure 1), located above MSL. Note that changes in *R*, determined by spatial (e.g. owing to heterogeneous subsidence rates) or long-term (induced by climatic changes) variations of RSLR rates, vertically shift *F*(*z*) and may vary the number and types of accessible equilibria, possibly triggering abrupt biogeomorphic regime shifts [3,20]. In particular, it is seen that thresholds exist above which the marsh and the tidal flat equilibria no longer exist (*R*>9.7 mm yr^{−1} and *R*>15.5 mm yr^{−1}, respectively, for the values of *H* and *C*_{0} assumed in figure 1). These thresholds signal transitions to regimes characterized by a reduced number of possible biogeomorphic patterns and mark limits to the resilience of tidal systems as we currently know them [2,20,39]. In the *multi-species* case (not reproduced here for brevity), because organic soil production increases monotonically with elevation, no stable marsh equilibrium is possible for *C*_{0}=20 g m^{3} and *R*=8.5 mm yr^{−1}, and only subtidal platform and tidal flat equilibria exist (a marsh equilibrium can be recovered, however, either by lowering the value of *R* or by increasing the value of *C*_{0}). This further illustrates the strong control exerted by vegetation on the tidal landscape (the interested reader will find a more detailed discussion of this point in [20–22]).

Competing stable states thus constitute the pattern formation mechanism for the large-scale biogeomorphic features characteristic of tidal environments worldwide. It is now interesting to explore the mechanisms that lead to the formation of well-known, smaller scale, patterns associated with marsh vegetation species distributions.

## 4. An one-dimensional model: marsh-scale biogeomorphic patterns

Vegetation *zonation* is possibly the most recognizable characteristic feature of tidal marshes over scales between a few centimetres and about 100 m. It is traditionally described as a spatial vegetation pattern in which patches composed of single species or of mixtures of species with approximately fixed composition display sharp transitions in the seeming absence of appreciable topographic gradients [8–10]. We show here that zonation is, in fact, a biogeomorphological pattern (as illustrated in figure 2), rather than simply a biological one, and that it is the visible symptom of the underlying feedbacks between biomass productivity and soil accretion.

### (a) Modelling framework

We summarize here the approach introduced in Marani *et al.* [22] to study the time evolution of a one-dimensional marsh transect perpendicular to the channel feeding the marsh with water and suspended sediment (figure 2)). This configuration may represent both a fringing marsh, delimited landward by a natural/artificial levee, or half of a marsh delimited by channels/creeks on both sides. The system is governed by the one-dimensional version of equation (2.1) (where the erosion of the surface is neglected as waves are effectively dampened by marsh vegetation), complemented by a hydrodynamic/sediment transport formulation and by a description of the biomass contribution to soil accretion.

The local inorganic sediment deposition rate is computed by integrating the instantaneous settling flux, *w*_{s}⋅*C*(*x*,*t*)/*ρ*_{b} (*w*_{s} and *ρ*_{b} having the same value assumed in the point model case) over a whole year, to obtain the mean annual settling rate *Q*_{s}(*x*,*t*). The trapping flux is here neglected, following Marani *et al.* [3] and Mudd *et al.* [27], who showed that particle settling largely dominates inorganic sediment deposition for flow velocities commonly observed in tidal marshes, i.e. less than 0.05 m s^{−1} in microtidal environments (i.e. areas where the tidal range is less than 2 m). The instantaneous SSC, *C*(*x*,*t*), needed for the evaluation of *Q*_{s}(*x*,*t*), is obtained by assuming an advective–dispersive sediment transport process (solved at time scales shorter than a tidal cycle) of the form
4.1where *D*(*x*,*t*)=*h*(*t*)−*z*(*x*,*t*) is the local instantaneous water depth; *h*(*t*) is the spatially uniform instantaneous water elevation obtained by assuming a semidiurnal sinusoidal fluctuation of the tidal level as described earlier; *k*_{d} is the dispersion coefficient (assumed here to be 1.5 m^{2} s^{−1} [40]); and *u*(*x*,*t*) is the local instantaneous fluid advective velocity, computed by assuming a quasi-static propagation of tidal levels from the continuity equation ∂*D*/∂*t*+∂(*Du*)/∂*x*=0. The separation of time scales between equations (4.1) and (2.1) allows us to decouple the two solutions: sediment transport is solved through a finite volume numerical method over a tidal cycle, by setting the following boundary conditions *D*(0,*t*)=*h*(*t*)−*z*(0,*t*) and *u*=0|_{x=L} in the continuity equation; and ∂*C*/∂*x*|_{x=L}=0 and *C*(0,*t*)=*C*_{0} in equation (4.1).

In order to evaluate the role of a set of competing vegetation species on marsh pattern formation, we express the local organic accretion rate as dependent on a species-specific *fitness function*: *Q*_{o}(*x*,*t*)=*γB*_{i}(*x*,*t*)=*γB*_{0}*f*_{i}(*z*(*x*)), ‘*i*’ being the species which happens to colonize site *x* at time *t*. 0<*f*_{i}(*z*)<1 summarizes the degree of adaptation of species *i* to elevation *z* and represents here the fraction of the maximum rate of biomass production, *B*_{0}, that species *i* is able to produce at elevation *z*. Marsh species are commonly observed to be maximally productive within specific ranges of optimal adaptation [6,7]. Following Marani *et al.* [22], we represent the decrease in productivity as elevation moves away from the optimal range by way of a simple analytical expression: , where *ζ*=*z*/*H*, corresponds to the maximum value of the fitness function for species *i*, and *λ* is a scale parameter, expressing the rate at which the fitness function tends to zero as elevation departs from the optimal range. A large value of *λ* represents specialized species, i.e. species which display a nearly maximum biomass productivity (and reproduction rate, as we shall see later) just within a narrow range of elevations. Conversely, a small value of *λ* represents species which are relatively well adapted to a broader range of marsh elevations.

Ecological competition for space is of course an important driver of vegetation distribution. Competition strategies in marshes, as in other biomes, involve investing energy and carbon, e.g. into increasing reproductive rates, or into an apparatus which increases the chance of plant survival (e.g. aerenchima, a system of conduits in the stem for ‘snorkelling’ air during tidal flooding, see [24]). It is reasonable to assume that such competitive abilities must operate most efficiently when a plant finds itself in the optimal environmental conditions for its species. Hence, we assume here that competition, i.e. the mechanism by which a generic site at elevation *z*_{k}=*z*(*x*_{k},*t*+*Δt*) may be colonized by species *j* at time *t*+*Δt*, while it was previously colonized by species *i* at time *t* (when elevation was *z*(*x*_{k},*t*)), should be driven by the value of the fitness function *f*_{j}(*z*_{k}) as compared to *f*_{i}(*z*_{k}). We consider two competition mechanisms. A *fittest-takes-all* mechanism, in which the species distribution throughout the system is updated at every (yearly) time step by placing at the generic site *x*_{k} the species *i* for which *f*_{i}(*z*_{k})>*f*_{j}(*z*_{k}) ∀*j*≠*i*. This mechanism allows the study of system equilibria and patterns in the ideal case in which stochastic heterogeneities (in the environment and in the organisms) affecting the outcome of competition are absent. We will see that this approach produces sharp biogeomorphic patterns, which serve to more clearly illustrate vegetation controls on marsh morphology. However, actual marshes are characterized by heterogeneous soil properties, individually varying plant responses and variable exposure to external forcings (wind, waves, insolation, nutrients, etc.), such that they cannot be expected to behave according to a perfectly regular the fittest-takes-all rule. In order to obtain a more realistic representation of actual marsh processes, we also consider the results of a *stochastic competition* mechanism, in which the ‘state’ of each site *x*_{k} is updated at each time step by randomly selecting species *i* with a probability .

The model simulations discussed below are started from a random distribution of individuals from a set of three species, whose fitness functions are preassigned in such a way as to span the elevation range between MSL and the maximum tidal level *H* in an approximately homogeneous manner. Several tests suggest that the salient properties of the final equilibrium configuration of the transect do not depend on the initial topographic profile, such that an initial topography which linearly decreases away from the channel is assumed in the following.

### (b) Stable equilibria and corresponding biogeomorphic patterns

We first consider a reference case (figure 3*a*) in which the platform is flooded by a tidal oscillation of period *T*=12 h and amplitude *H*=0.50 m and is forced with a prescribed SSC in the channel *C*_{0}=20 mg l^{−1} and a rate of RSLR *R*=3.5 mm yr^{−1} (e.g. consistent with the local rate of RSLR for the twentieth century in the Venice Lagoon [41]). We consider three vegetation species with the same degree of specialization (*λ*=5), which differ for different optimal elevation values ( m, m and m above MSL) for the ‘red’, ‘green’ and ‘blue’ species (dashed black line in figure 3*c,f,i*). Vegetation competition is modelled according to the fittest-takes-all mechanism. Theoretical analyses of the fixed points of equation (2.1) and simulations indicate that a stable steady-state configuration exists. Such configuration exhibits sharp transitions between neighbouring biogeomorphic terrace-like structures (figure 3*a*), which display very mild within-patch slopes and are colonized by a single vegetation species (emphasized by colour-coding in figure 3*a*). These zonation patterns are very reminiscent of the typical features observed in actual salt-marsh landscapes, and in particular of the sharp boundaries between vegetation patches observed in nature [8–10].

Following Marani *et al.* [22], the physical and biological processes responsible for the development of zonation patterns in salt-marsh landscapes can now be identified through the lens of the quantitative framework described. In particular, one finds that zonation structures emerge as a result of ecogeomorphic feedbacks involving inorganic deposition and organic accretion (via biomass production), which lead to pairs of stable and unstable equilibrium states.

Let us consider, for the purpose of illustration, the feedbacks that produce the geomorphic structure located at mid-elevations (species *i*=2, identified by m, in figure 3*a*, in green in the online version), and in particular the equilibrium of the lowest site within this patch (with coordinate, say, ). The equilibrium condition ∂*z*/∂*t*=*Q*_{s}(*z*)+*γ*⋅*B*_{0}*f*_{2}(*z*)−*R*=0 gives two solutions: a stable equilibrium above the maximum of *f*_{2}(*z*) ( m above MSL), with elevation m above MSL (solid circles in figure 3*b*,*c*, in green in the online version), and an unstable equilibrium below the maximum of *f*_{2}(*z*), located at m above MSL (open circles in figure 3*b*,*c*, in green in the online version).

Note that *Q*_{s}(*z*) is found (through numerical simulation) to be a monotonically decreasing function of *z* because an increase in *z*, by decreasing the amount of the ‘precipitable’ suspended sediment *C*⋅*D*=*C*⋅(*h*−*z*), decreases the amount of inorganic sediment deposited during a tidal cycle. Therefore, the existence of two solutions to ∂*z*/∂*t*=0 is not because of variations in *Q*_{s}(*z*), but is a consequence of the presence of a maximum in *f*_{2}(*z*), and hence of biological controls.

To determine the stability properties of the equilibrium state in (figure 3*a*), we numerically find the value of ∂*z*/∂*t* from equation (2.1) for a set of transect configurations obtained by perturbing the local elevation in within a neighbourhood of and (figure 3*b*), and we look at the sign of the changes in the sediment fluxes brought about by the perturbation.

A positive perturbation above (, figure 3*b*) generates a decrease in biomass production (figure 3*c*) and in the related organic accretion rate (because d*f*_{2}(*z*)/d*z*<0), together with a decrease in the inorganic deposition rate, *Q*_{s} (see discussion above). Because both *Q*_{s} and *Q*_{o} decrease with respect to their equilibrium values, the accretion rate becomes negative, ∂*z*/∂*t*<0 (figure 3*b*), thus acting to dissipate the initial perturbation and bringing the system back to the original equilibrium. Similarly, a negative perturbation below (i.e. ) induces an increase in the total (organic + inorganic) accretion rate, such that ∂*z*/∂*t*>0 (figure 3*b*), again acting against the imposed perturbation and bringing the system back to the equilibrium state. is thus a stable equilibrium, a condition which can in summary be expressed as ∂(∂*z*/∂*t*)/∂*z*<0. A similar reasoning can be applied to all other sites belonging to the same vegetation patch (‘green’ in the online version), to find that these are all stable equilibrium states, whose elevation slowly increases from right to left because of a spatial increase in the inorganic deposition rate, *Q*_{s}, when approaching the channel, the source of inorganic sediment.

The second equilibrium solution, , is also important, as it is involved in the pattern-generating mechanism. Let us imagine site having elevation . A positive perturbation bringing the local elevation to enhances the organic accretion rate, *Q*_{o} (figure 3*c*), while reducing the inorganic deposition rate, *Q*_{s}. However, one finds numerically that the decrease in *Q*_{s} is small compared with the increased organic sediment production, such that ∂*z*/∂*t*>0 (figure 3*b*). This means that a positive perturbation is destined to grow, eventually leading the local elevation to . Similarly, a negative perturbation bringing the local elevation to generates a decrease in the organic accretion rate (because d*f*_{2}(*z*)/d*z*>0) which is faster than the increase induced in the inorganic deposition rate, *Q*_{s}, thus making ∂*z*/∂*t*<0 (figure 3*b*). Negative perturbations to are also destined to grow, thus pushing the elevation to , the next lowest stable equilibrium state, associated with species *i*=3 (solid circle in figure 3*b*, in blue in the online version).

The pattern-generating mechanism is thus clarified: step-like structures emerge because unstable states, controlled by the shape of the fitness functions, act as repellers, pushing topographic elevation either to the immediately higher stable equilibrium or to the next lowest one. This picture, highlighting the role of vegetation species in tuning soil elevation within ranges in which they are preferentially adapted, is valid when soil organic production is more sensitive to variations in *z* than the inorganic sediment deposition, such that |∂*Q*_{s}/∂*z*|≪|∂*Q*_{o}/∂*z*|. In this case, ∂*Q*_{o}/∂*z* controls the stability of the equilibria.

On the contrary, when |∂*Q*_{s}/∂*z*|≈|∂*Q*_{o}/∂*z*|, equilibrium stability is jointly determined by inorganic deposition and biomass production. This is for example the case for the part of the transect adjacent the channel in figure 3*a*. One may, in fact, note that the equilibrium state in the upper marsh zone (solid circle in figure 3*b*,*c*, in red in the online version) lies on the branch of the fitness function located below the maximum (, dashed black line in figure 3*c*), which would seem to imply an unstable equilibrium because d*f*_{1}(*z*)/d*z*>0. However, because |∂*Q*_{s}/∂*z*|>|∂*Q*_{o}/∂*z*|, then ∂/∂*z*(∂*z*/∂*t*)<0, implying a stable equilibrium controlled by inorganic sediment deposition.

### (c) The role of the rate of relative sea-level rise

We can now analyse the changes induced in biogeomorphic patterns when the marsh is subjected to increasing rates of RSLR. We explore here *R*=5 mm yr^{−1} and *R*=7 mm yr^{−1} for comparisons with the reference value *R*=3.5 mm yr^{−1}, representative of many tidal environments under current conditions. In the case of *R*=5 mm yr^{−1}, the biogeomorphic features of the equilibrium patterns are qualitatively similar to those of the reference case, although some differences emerge (figure 3*d*). The comparison of the two equilibrium profiles shows that as the rate of RSLR increases, the equilibrium marsh elevations tend to decrease, in agreement with the results of other modelling approaches [2–4,21,39]. It should be noted that, under increased RSLR conditions, the lowest sites within the two highest geomorphic structures, colonized by species *i*=1 and *i*=2 (figure 3*d*, in red and green respectively in the online version), lie on the branch of the fitness functions where d*f*_{i}/d*z*>0 and hence d*Q*_{o}/d*z*>0 (solid circles in figure 3*f*), very close to the maximum value (dashed black lines in figure 3*f*). Because ∂(∂*z*/∂*t*)/∂*z*<0 (solid circles in figure 3*e*), this suggests that the stability of the equilibria are here governed by the inorganic deposition flux, i.e. ∂*Q*_{s}/∂*z*>∂*Q*_{o}/∂*z*. On the contrary, the lowest equilibrium elevation in the lowest patch, colonized by species *i*=3, belongs to the branch of the fitness curve where d*Q*_{o}/d*z*<0 (solid circle in figure 3*f*, in blue in the online version), and is therefore governed by the organic flux, as in the reference case (i.e. ∂*Q*_{o}/∂*z*>∂*Q*_{s}/∂*z*). Interestingly, the general lowering of the equilibrium elevations induced by an increased rate of RSLR causes a sort of cascade effect, in which areas that, for a lower value of *R*, belonged to higher elevation structures now abrupt transition to lower equilibrium states. As a result, the area colonized by the species adapted to the lowest elevations is larger under increased *R*, with a corresponding reduction in the extent of the higher vegetation patches. A further increase in the rate of RSLR, however, may drastically change marsh patterns, owing to the lowest sites in the lowest patch becoming unstable and switching abruptly to a yet lower (unvegetated) equilibrium state (figure 3*g*,*h*,*i*). Part of the marsh has transitioned to a tidal flat equilibrium state, a transition suggested also by the point model. It is important to note, however, that the rates of RSLR suggested by the two models are rather different and refer to different phenomena. The point model, in fact, can only deal with equilibria of a tidal platform as a whole. The one-dimensional model in this case shows that for relatively low threshold rates of RSLR (approx. 6 mm yr^{−1} versus 7.2 mm yr^{−1} obtained from the point model, considering *H*=0.50 m and *C*_{0}=20 mg l^{−1}) part of a marsh can collapse to a tidal flat equilibrium with the consequent loss of the associated biodiversity and ecosystem services. Interestingly, this finding points to the existence of a maximum sustainable length of vegetated marshes dependent on the RSLR (and, similarly, on the local SSC in the feeding channel, not explored here).

### (d) The role of vegetation specialization

We have shown that vegetation acts as a landscape engineer by producing particular morphological–biological patterns in salt-marsh ecosystems. One might then wonder how vegetation characteristics influence the biogeomorphic features of salt-marsh landscapes. To address this question, we have considered vegetation species with different degrees of specialization, i.e. species capable of adapting to narrower or broader ranges of marsh elevations, as represented by different values of the scale parameter *λ*. In particular, we consider here vegetation species which are both more specialized (*λ*=10) and less specialized (*λ*=2) than the reference case (*λ*=5). The maximum productivity elevations are in all cases the same as in the reference case.

Model results show that when vegetation species are assumed to be more specialized (figure 4*a*–*c*) the emerging biogeomorphic structures are more ‘focused’ within narrow elevation ranges with smaller topographic gradients, the evidence of a stronger vegetation control on morphodynamic processes. Such a behaviour can be illustrated by comparing the slope of the fitness functions, d*f*_{i}/d*z*, for the reference case and the case of more specialized vegetation species in figure 4*c*. As the degree of specialization increases, the slope of the fitness function increases (compare dashed and solid lines in figure 4*c*). Therefore, smaller variations in the topographic elevation are required to determine the changes in *Q*_{o} necessary to match the spatial gradients of *Q*_{s} and thus balance the rate of RSLR. When less specialized (i.e. *λ*=2) vegetation species are considered, topographies displaying smoother, and possibly more realistic, transitions between patches are obtained (figure 4*d*–*f*). Interestingly, because of the reduced vegetation control on elevation, within-patch topographic gradients are greater for less specialized vegetation, suggesting the possibility of reading the morphological signature of vegetation physiological adaptations to environmental conditions. Our results, in fact, suggest that the slope of a vegetation patch may be indicative of the breadth of the range of elevations to which the colonizing vegetation species is adapted.

### (e) Stochastic competition

We finally address the possible role of environmental stochasticity, e.g. associated with heterogeneities in the soil characteristics or with the variability of individual plant properties. More realistic and less regular patterns emerge when the outcome of competition is modelled with a stochastic rule (figure 5*a*). Even though the marsh profile decreases as one moves away from the marsh border, as in the reference case, the transect is now characterized by a somewhat irregular trend and vegetation patches previously populated by single species are replaced by patches of mixed vegetation species. However, the underlying biogeomorphic feedbacks, responsible for forming the observed zonation patterns, may still be identified by studying the frequency distribution of topographic elevation along the transect (figure 5*b*): the frequency density of marsh elevations is multi-modal and each frequency maximum is robustly associated with just one vegetation species, as emphasized by colour-coding the bars of the graph according to the most abundant species colonizing each elevation interval. Hence, even though topographic and vegetation patterns are not visually evident in this case, the tight relationship between peaks in the frequency distribution and vegetation species is found to be a characteristic signature allowing the detection of the landscape-constructing role of biota [22].

## 5. Discussion and conclusion

We have analysed the spatial distribution of biogeomorphic patterns characterizing tidal environments over a wide range of scales. Our key result is that multiple competing stable states, governed by two-way biogeomorphic feedbacks, provide a unifying framework to explain the formation of patterns in tidal environments from the marsh to the system scale. Furthermore, we propose that throughout such a wide range of spatial scales vegetation is not just passively adapting to morphological features prescribed by sediment transport, but rather it is fundamentally tuning tidal patterns according to its physiological requirements. Vegetation is indeed the ‘Secret Gardener’ shaping the tidal landscape.

A conceptual point model, describing the time evolution of the mean elevation of a tidal platform, shows that subtidal platforms, tidal flats and salt marshes are possible equilibrium states resulting from the interplay of physical and biological processes. The existence of these competing stable states represents the large-scale pattern formation mechanism accounting for the emergence of subtidal platforms, flats and marshes, characteristic and ubiquitous biogeomorphic features of tidal landscapes worldwide.

A one-dimensional description of tidal biogeomorphology shows that biogeomorphic feedbacks are also responsible for the occurrence of smaller scale patterns associated with the spatial distribution of marsh vegetation. Zonation patterns are shown to emerge as a result of the feedbacks between inorganic deposition and organic accretion, via biomass production. We highlight the crucial role of vegetation in shaping the tidal landscape: marsh vegetation acts as a landscape engineer by tuning platform elevation, leading to the formation of zonation structures.

The tight feedback between primary productivity and topography suggests the possibility of decoding the signature of physiological plant adaptations in observed tidal marsh morphologies. In fact, we show here that specialized vegetation species, highly fit only within a narrow range of marsh elevations, exert a strong control on topography and constrain elevation within similarly narrow ranges. By contrast, vegetation species which are relatively well adapted to a broader range of marsh elevations more loosely tune the marsh landscape and produce more gradual transitions between adjacent vegetation patches. Vegetation controls on topography also hold in the presence of more realistic heterogeneous environmental forcings, such as soil properties and marsh microtopography. Biogeomorphic patterns are less sharply defined in this case, but the underlying presence of multiple equilibria as the pattern-generating mechanism can still be detected via the presence of multiple peaks in the frequency distribution of topographic elevation. The one-to-one correspondence between such peaks and species abundance is the unique signature of vegetation as a landscape constructor.

Possible biogeomorphic equilibria and the associated patterns are sensitive to changes in the forcings. Changes in the rate of RSLR produce the migration of marsh biogeomorphic patterns, eventually leading to their selective or complete disappearance when a critical rate of RSLR is exceeded. The results of the point model and the one-dimensional views of tidal biogeomorphology qualitatively agree in showing the possibility of threshold behaviour and abrupt shifts among alternative equilibria. However, the spatially extended description of the relevant biogeomorphological processes suggests that part of a marsh platform may transition to a tidal flat equilibrium well before the marsh as a whole ceases to exist. This result points to the possibility of widespread marsh degradation under the rates of sea-level rise projected for the coming century (up to 20 mm yr^{−1} [42] of eustatic SLR, not inclusive of local subsidence).

## Funding statement

We acknowledge support by Duke University, the PhD School of Civil and Environmental Engineering Sciences at the University of Padova and the Italian National Project ‘Eco-morfodinamica di ambienti a marea e cambiamenti climatici’. A.D’.A. thanks the CARIPARO Project titled ‘Reading signatures of the past to predict the future: 1000 years of stratigraphic record as a key for the future of the Venice Lagoon’.

## Acknowledgements

We thank Massimiliano Ignaccolo for discussions on the mathematical representation of vegetation fitness functions.

## Footnotes

One contribution of 13 to a Theme Issue ‘Pattern formation in the geosciences’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.