Explaining multiple patches of aquatic vegetation through linear stability analysis

Aquatic vegetation is ubiquitous in lowland rivers, and it is typically present in the shape of spatial self-organized patches of biomass. In this work, we mathematically define the threshold conditions for the incipient formation of self-organized vegetated patterns in the shape of central or multiple row patches. The analysis is carried out through a linear stability analysis whereby the 2D eco-hydrodynamic model is linearized and the growth rate of small-scale perturbations is evaluated considering a basic state represented by an initially uniformly vegetated and straight channel having a certain aspect ratio and Froude number. Results illustrate that, for given vegetation properties, instability arises when both the Froude number and the aspect ratio are higher than a given threshold; in this case, self-organization occurs, and spatial patterns of patches will develop according to the wavelength associated to the maximum growth rate. Moreover, instability and self-organization take place when the undisturbed vegetation density is lower than upper bound; this suggests that densely vegetated channels, as in the case of rivers populated by invasive species, will not experience the formation of any spatial patterns.


Introduction
Aquatic vegetation is ubiquitous in water bodies providing many beneficial ecological services such as creating favorable habitat for fauna, improving water quality, reducing sediment turbidity, stabilizing sediment bed [23].
The so-called hydro-bio-geomorphic feedbacks play a fundamental role in the formation and plano-altimetric evolution of many different types of landscapes including alluvial rivers [28], intertidal areas [30], eolian dunes [2] and drylands [5].
In the case of the aquatic vegetation in lowland rivers, surprisingly, very little is still known regarding the mutual interaction between flow hydrodynamics and biological processes which might produce self-organized vegetated patterns [16]. In particular, various studies suggest that vegetation can grow over an initially bare bed in the form of isolated patches [17,27]. These patches can further expand either in the shape of self-organized patterns characterized by densely vegetated bed areas separated by unvegetated channels, or in a uniformly homogeneous distributed vegetation covering the entire channels. The growth and evolution of patches is dictated by interactions with the flow field: the divergence of flow around isolated patches can produce a laterally accelerated velocity field eroding the channel bed and further preventing the expansion of neighboring patches. However, when the distance between patches is relatively small, the side flow acceleration becomes negligible thus promoting sediment deposition, vegetation growth and, eventually, patch merging [30]. Additionally, the interactions between the wakes of vegetated patches might produce positive feedback for their streamwise growth [32].
Cornacchia et al. [15], while studying patch dynamics of Ranunculus penicilatus macrophythes in a lowland stream, discovered that patches organize themselves in V-like shapes to reduce drag forces. Importantly, self-organization of vegetation according to quasi-periodic spatial discontinuous patterns, represents an adaptive capacity of rivers to withstand a wider range of discharges than if it were homogeneously uniformly distributed over the entire river bed [16].
Recently, Calvani et al. [9], for a given vegetation species, formulated an analytical approach to define the threshold conditions between vegetated and sediment-covered river beds. In the case of Froude numbers below the critical value for removal of vegetation, Calvani et al. [8] extended the approach by Bärenbold et al. [3] and Crouzy et al. [18] to mathematically disclose the mutual interactions between flow field and vegetation evolution through a linear stability analysis. In analogy with the well-known case of free sediment bars developing in cohesionless straight channels [14], the formation of free vegetated patterns, consisting of migrating repetitive alternate sequences of isolated patches, was mathematically predicted as a function of the hydrodynamic parameters (e.g., Froude number) and on the specific plant-species properties, in terms of growth, spread and decay rates and carrying capacity, as well as the initial vegetation density.
In this work, we further expand the analysis of Calvani et al. [8] by mathematically investigating self-organized patchiness associated to the formation of more complex spatial distributions. In particular, we here focus on vegetated patches associated to higher transverse structures, such as central or multiple row patches, that can be typically found in lowland rivers; examples are the Zwarte Nete in Belgium [27] and various rivers in the UK (Fig. 1) .

Formulation of the problem
The starting point of our analysis is the two-dimensional form of the governing equations of eco-hydrodynamics in a straight channel with a non-erodible bed. By analogy with the naming adopted by Bärenbold et al. [3], we label eco-hydrodynamic the minimal model composed of vegetation dynamics coupled with the classical two-dimensional Saint-Venant and continuity equations for hydrodynamics. We consider a flow in a straight channel with a rectangular cross-section and fixed lateral banks; the bed is non-erodible, actually sediment transport may take place but no net aggradation and degradation is accounted for. Indeed, we decide not to include a sediment transport equation in the model in order to focus on the vegetation dynamics, with the aim of addressing the eco-morphodynamic problem in a further step. Finally, the non-erodible bed presents favourable condition for plant colonization [9], and at equilibrium it is covered by submerged aquatic vegetation with a uniform density. The set composed by the fluid mass density , the uniform depthaveraged flow velocity U * 0 and depth Y * 0 , the river width W * =2B * (with B * the half-width of the river), and the carrying capacity * m is used for the non-dimensionalization of the equations. We here outline the physical meaning of carrying capacity that is actually very well known in ecology, but much less classical in river hydraulics: * m is the maximum density or biomass of a population that can live in a particular habitat given the limited resources of the habitat itself (for population modeling in ecology the reader is referred to Hartvigsen [20] , whereas we refer to Camporeale and Ridolfi [10] and Bärenbold et al. [3] for vegetation modeling in river environments using the concept of carrying capacity). As to notation, here and in the following, an asterisk indicates dimensional variables and the subscript 0 denotes variables evaluated at the base state consisting of normal flow. Further details on dimensionless quantities can be found in Bärenbold et al. [3] and Calvani et al. [8].

The flow model
The governing equations for the water flow are the steady depth-averaged two-dimensional Saint-Venant and continuity equations: where �� ⃗ V = {U, V} is the flow velocity, Y is the flow depth, s and n are the longitudinal and transverse coordinates, respectively, Fr is the Froude number, is the bed elevation, = B * ∕Y * 0 is the half-width to depth ratio (also labelled aspect ratio), and the total shear stress. A closure relation for the total shear stress is needed and it is built using the Chézy nondimensional conductance coefficient, {U, V} . We adopt for the Chézy coefficient the relationship proposed by Luhar and Nepf [21] which accounts for the presence of submerged vegetation: where C f = (5.75 log 10 (2Y * 0 ∕d * )) −2 is the bed friction coefficient [31] being d * the sediment grain size, plant height, * m is the carrying capacity, and is the vegetation density (the dimensional vegetation density, * , is evaluated as the number of plants per unit area, and is made dimensionless by means of the dimensional carrying capacity, * m ). Values of vegetation characteristics, together with hydraulic conditions, are representative of dense canopy (i.e., 1 , according to Nepf [23]). The corresponding value of solid volume frac- is in the range analyzed by Zong and Nepf [33] and Ciraolo et al. [11] in flume experiments. Additionally, the vegetation height, h * v , represents the effective height of plants after reconfiguration, which is dependent both on flow conditions and plant properties such as stiffness and life stage [21]. For the sake of simplicity, we consider h * v to be constant, without losing the general applicability of the proposed framework. Values of the parameters involved in Eq. (4 ), and in the problem in general, are reported in Table 1.
We anticipate that vegetation dynamics is much slower than hydrodynamics, therefore, equations 1-3 are in the steady form (see Calvani et al. [8]). Finally, we point out the boundary conditions for the governing equations: by considering an infinitely long channel, we can assume periodic boundary conditions in the streamwise domain; whereas in the transverse direction we consider a rectangular cross-section, we impose a vanishing flux of water between the channel and the lateral banks, and we assume the vegetation to be confined in the wetted section.

The vegetation dynamics model
Vegetation dynamics on the river bed is modeled through an equation describing changes in time and space of the plant density, (Eq. 5). The equation of vegetation dynamics reads: where, in addition to the variables already introduced, we find the param- Fr g 1∕2 which are the dimensionless parameters for growth, diffusion and decay, respectively, being * g , * P , * D , and * d their relative dimensional coefficients (reported in Table 2 after Calvani et al. [8]). Such dimensional coefficients are rates that characterize the behaviour of the processes of growth, decay, diffusion; for instance high values of the decay rate, * d , model plants prone to uprooting, whereas low values of * d model plants more resistant to uprooting [7] . These rates mainly depend on the vegetation species and on the hydraulic conditions. For a detailed analysis of these vegetation parameters and coefficients the reader is referred to Calvani et al. [9] and Calvani et al. [8].
On the right-hand side, the equation contains three terms. The first one models the growth of vegetation density, usually represented in condition of limited resources by a logistic function [10]. The second term represents vegetation diffusion, namely the spatial spread of plants on the river bed described, as usual, through the laplacian operator [18, 19, 26, among others]. While the first term models the growth of vegetation density according to the density itself in a particular cell, the diffusion increases or decreases such growth according to the density in neighbouring cells. The diffusive term models the ensemble of processes that ecologists consider among positive feedbacks [27]. For instance, stress reduction, which may promote sediment deposition and eventually colonization, represents a short-ranged positive-feedback effect exerted by a patch on its neighbouring plants [30,32]. The last term on the right-hand side of equation 5 models the plant decay by accounting for vegetation mortality due to flow uprooting [24]. Equation 5 presents a rather complex growth term, which is able to describe how the growth of vegetation density depends on the deposition and resprouting of vegetative propagules generated by uprooting [8]. The underlying mechanism is that not all the uprooted plants die, some of them can resettle and resprout [4]; we performed a mass balance of these plants, i. e. the propagules, and it was included in Eq. 5 since it influences the growth term. If resettlement dynamics of propagules is neglected, the growth term accounts for the presence of limited resources only, and the growth term simplifies to g (1 − ) , and the whole equation of vegetation dynamics takes the form of the standard equation [3,18]. However, Calvani et al. [8] demonstrated that the resettlement of uprooted propagules is essential for the formation of vegetation patches in the stability analysis framework.

Linear stability analysis
The equations of the water flow and the vegetation dynamics are coupled and together form the eco-hydrodynamic model; the whole set of equations is linearized by expanding all the variables as the superposition of a base state plus some small periodic perturbation where the parameter is chosen to be small in accordance with the assumption of small perturbations of the base state; precisely, the subscript 0 denotes the base state ( O( 0 ) )  ) ). In Eqs. 6 and 7, Ω = Ω R + i Ω I represents the complex wave speed of the perturbation, k s and k n are the streamwise and transverse wavenumber, respectively, and c.c. stands for complex conjugate. In particular, the wavenumbers k s and k n are non-dimensional quantities expressed by the relationships k s = 2 B * ∕L * s and k n = 2 B * ∕L * n , where L * s and L * n are the physical wavelengths in the streamwise and transverse directions, respectively.
At the base state the variables {U 0 , V 0 , Y 0 , 0 } assume the values {1, 0, 1, 1 − h v d ∕ g } describing normal flow conditions with uniform vegetation density. At the linear level the following algebraic eigenvalue problem is recovered [8]: with and where the complex quantity Ω is the eigenvalue of the vegetation dynamics. In the system, C 0 is the Chézy coefficient (Eq. (4)) and C 1,Y = C Y | 0 and C 1, = C | 0 are its partial derivatives in Y and , respectively. The subscript 0 denotes that the three terms are evaluated at the base state i.e., Y = Y 0 , = 0 . The solution of the problem (9) is the so called dispersion relationship, a linear algebraic expression of unknown Ω . The real part of the eigenvalue, Ω R , represents the growth rate of the perturbation and the imaginary part of the eigenvalue, Ω I , the migration rate of the perturbation. Positive values of Ω R indicate growing perturbations, and thus unstable base state and stable vegetation patches; on the other hand, negative values of Ω R indicate stable equilibrium conditions and therefore no spatial patterns are expected to form. Positive (negative) values of Ω I denote upstream (downstream) migration of the perturbation.
Calvani et al. [8] provided the analysis of the growth rate and the migration rate for the pattern mode m equal to 1 and therefore lateral boundary conditions imposed at the riverbanks posed by k n = m ⋅ ∕2 = ∕2 . This condition corresponds to alternate patterns, as in classical stability analysis of river sediment bars [14]. In this work we extend such analysis by exploring patterns of higher modes, precisely m = 2 and m = 3 corresponding to k n = and k n = 3 ∕2 , respectively. These conditions match field observations of multiple vegetation patches in the river cross-section, as showed in Fig. 1. In this case too, the analogy with river sediment bars is clear (for instance, see the recent work by Ali et al. [1] that addresses instability modes from 1 to 4). For the aim of comparison, we also report the solution for m = 1 . A schematic representation of the the three pattern modes is reported in Fig. 2 ; the figure also illustrates the condition of patches' absence, namely the uniform T vegetation cover (top plot, m = 0 ). Finally, in this study we do not report information on the patches' migration rate because no relevant differences with respect to Calvani et al. [8] were found.

Vegetation patches for different pattern mode scenarios
The main outcome of stability analyses are the so called stability plots, graphics that illustrate the values of the growth rate of the perturbation in the parameter space, namely in ranges of the relevant parameters of the problem [3, 12-14, 22, among others]. Figure 3a reports the stability plots for the three pattern modes m = 1, 2, 3 in the − k s plane. Precisely, the marginal curves, namely the curve of null growth rate, are reported for m = 1, 2, 3 by means of thick solid lines in black, blue, and green respectively; the same legend of colors is used for the dashed lines representing the curves of maximum amplification rate. For each solid curve the plot must be read as follows: the solid curve is the locus of points where the perturbation neither grows nor decays ( Ω R = 0 ); the concave region outside the curve is the region where the perturbation decays ( Ω R < 0 ); the convex region inside the curve is the region where the perturbation grows ( Ω R > 0 ). This last zone of the plot is called unstable region since there the equilibrium state is unstable and thus the formation of spatial patterns is predicted. Inside the unstable region, the curve of maximum amplification rate ( Ω Rmax ) identifies, for a given set of values of , Fr, 0 , and m , the most unstable wavenumber k s , thus the analysis predicts the pattern wavelength that it is the most probable to form. The same line of reasoning holds true when interpreting the curves in Fig. 3b, where the real part of the vegetation eigenvalue is reported in the Froude number vs wavenumber plane.
Both Fig. 3a, b show that for all the three pattern modes (i) the unstable regions are bounded; (ii) there are critical lower thresholds both for and Fr and the mode 2 (3) thresholds have values that are twice (three times) the value of mode 1 threshold, like in sediment bar morphodynamics [1, among others]; (iii) a wavelength selection mechanism is provided by the analysis as highlighted by the maximum growth rate curves, each of them also showing that the range of most unstable wavenumber is narrow. For a more detailed analysis of the unstable region for the first mode m = 1 we address the reader to Calvani et al. [8], where we provided the stability plots for m = 1 only in the same parameter spaces of Figs. 3 and 4 of this work, in which there are no different results with respect to the m = 1 analysis of Calvani et al. [8] , indeed we recover here the same unstable regions; however in the present work we just report marginal curves ( Ω R = 0 ), whereas in Calvani et al. [8] we also represented unstable regions ( Ω R > 0 ) through colors and isolines. In this work, we focus our attention on the result that the three curves of maximum amplification rate are not coincident and the first mode marginal curve does not contain the marginal curves of higher modes. This is actually different from what is observed in classical stability plots of river bars morphodynamics where the three curves of maximum growth rate coincide and the unstable region of m = 1 contains the unstable regions of higher modes [3,6]. Instead, as regards the formation of higher-mode patterns of vegetation, Bärenbold et al's anlysis in conditions of fixed vegetated bed (just vegetation instability, no morphodynamic instability) does not predict the formation of multiple patches ( m > 1 ); however this outcome is not duly substantiated and the overall analysis is impaired by unrealistic, extremely high values of the vegetation parameters (growth, decay, and diffusion rates) that make the vegetation dynamics unrealistically rapid [8] .
The observations on the unstable regions as well as the maximum amplification rate curves of Fig. 3 can be extended to the stability plots reported in Fig. 4. In particular, Fig. 4a, b show the values of the real part of the vegetation eigenvalue in the 0 − k s space and in the m − k s space, respectively. A peculiarity that is present only in the stability plots in the 0 − k s plane is that each unstable region present an upper limit in terms of vegetation density 0 . Above this upper threshold no vegetation patches are expected to form, and the uniform vegetation cover at the base state is stable. In other words, it exists an upper threshold in terms of initial vegetation density ( 0 ≃ 0.42, 0.42, 0.34 for m = 1, 2, 3 respectively) for the onset of vegetation patches [8]. This may explain the fact that a very dense vegetation cover is stable to perturbation and permanently resists to a wide range of hydraulic conditions, as it is attested by field studies on the development of vegetation cover by invasive species [25,29] . Such field studies give evidence of the absence of vegetation patches, conversely uniform, resistant colonization is observed. Our mathematical finding, graphically represented by the upper thresholds of 0 in Fig. 4 a, agrees with such field observations; this outcome, already presented in Calvani et al. [8] for m = 1 , holds true also for higher pattern modes. We finally point out that in figure 4 as well as in Fig. 3 the curves of maximum growth rate are distinct and never superimposed.
The condition of distinct curves of maximum amplification rate, however, enables to always determine which is the most unstable pattern mode, as it is illustrated in Fig. 5. Figure 5 shows The three examples demonstrate that for every set of values of the relevant parameters it is possible to determine which mode of instability, 1, 2 or 3, has the larger growth rate. Figure 5 shows again the consequence of distinct curves of maximum amplification rate and first mode unstable region not containing unstable regions of higher mode (Figs. 3, 4). Indeed, there is a region of the parameter space where the instability mode 2 only is unstable and the instability mode 1 is stable (see Fig. 5b, 2.4 ≤ ks ≤ 3.5 ); similarly there is a region in which the pattern mode 3 is unstable and patterns mode 1 and 2 are stable (see Fig. 5c, 4.0 ≤ ks ≤ 4.6 ). This shows again the difference with respect to sediment bar modes: in the classical bar morphodynamics if mode 2 is unstable also mode 1 is unstable, as well as if mode 3 is unstable also modes 1 and 2 are unstable [Ali et al. 1, among others].

Conclusion
In this work we presented a linear stability analysis of 2D eco-hydrodynamic problem to identify the threshold conditions for the formation of self-organized patterns of patches of vegetation in straight channels. In particular, we investigated spatial patterns associated to higher order transverse structures, such as central or multiple row patches. The main conclusions can be summarized as follows: • Instability and self-organization occur when both the Froude number and the aspect ratio are higher than a given threshold. As in the classical linear-stability analyses, this threshold increases linearly with the transverse mode. As an example, the formation of central type of patterns (mode = 2) might be expected in channels having a Froude number higher than about 0.3 and aspect ratio greater than about 15, depending on vegetation properties. • The most unstable longitudinal dimensionless wavenumbers are in the order 1 indicating a longitudinal spacing of about 3 times the channel width; • Instability presents an upper bound for higher vegetation densities; thus, self-organization will likely not occur in those vegetated channels subject to a dense cover such as in the case of invasive species. In the case of Froude number of 0.3 and aspect ratio of 30, we found that the upper limit is in the order of 0.3, for given vegetation properties.
Future developments might regard field monitoring of aquatic vegetation in lowland rivers to provide further quantitative evaluations of the parameters (such as decay, diffusion and growth coefficients) that enters the equation describing the dynamics of vegetation density.