Shear-induced instabilities of flows through submerged vegetation

TL;DR

This paper develops a fluid-structure interaction model for submerged vegetation flows, revealing how shear at the canopy surface influences flow instabilities like monami, and identifies conditions for stability based on plant flexibility and density.

Contribution

It introduces a reduced framework coupling elastic beam models with fluid equations to analyze flow instabilities in submerged vegetation, providing new insights into the onset of monami.

Findings

Shear at the canopy surface is key to instability onset.

Sufficient plant flexibility or sparsity stabilizes the flow.

The model predicts conditions for stable and unstable vegetative flows.

Abstract

We consider the instabilities of flows through a submerged canopy, and show how the full governing equations of the fluid-structure interactions can be reduced to a compact framework that captures many key features of vegetative flow. By modelling the canopy as a collection of homogeneous elastic beams, we predict the steady configuration of the plants in response to a unidirectional flow. This treatment couples the beam equations in the canopy to the fluid momentum equations. Our linear stability analysis suggests new insights into the development of instabilities at the surface of the vegetative region. In particular, we show that shear at the top of the canopy is a dominant factor in determining the onset of instabilities known as monami. Based on numerical and asymptotic analysis of the generalised eigenvalue problem, the system is shown to be stable if the canopy is sufficiently sparse or if the plants are sufficiently flexible.

Tables

Table 1. Table 1: Modelling approaches of a selection of previous work on flow through a fully-submerged canopy or a single plant. The ‘Coupling’ category states whether the stability analysis takes into account of both flow perturbations and plant reconfiguration.

	Flow	Plant model	Stability analysis	Coupling
Alben et al. (2002)	solved	elastic strip	no
Ghisalberti & Nepf (2004)	solved	rigid cylinder	no
Poggi et al. (2004)	solved	rigid cylinder	no
Dupont et al. (2010)	solved	mechanical oscillator	yes	yes
Luhar & Nepf (2011)	imposed	elastic strip	no
Luminari et al. (2016)	solved	rigid cylinder	yes	no
Singh et al. (2016)	solved	rigid cylinder	yes	no
Zampogna et al. (2016)	solved	rigid cylinders	yes	no
	solved	rigid porous medium	yes	no
Sharma et al. (2017)	solved	dynamic cluster	yes	yes
	solved	rigid porous medium	yes	no
Leclercq & de Langre (2018)	imposed	elastic beam	yes	no
This work	solved	elastic beam	yes	yes

Table 2. Table 2: A summary of the dimensionless parameters in the governing equations of flow through a homogenised canopy ( 18 ).

	Symbol	Expression
Reynolds number	$\text{Rey}$	$\rho UH/\mu$
Froude number	$Fr$	$U/\sqrt{gH}$
Submergence ratio	$h$	$l/H$
Canopy density	$\lambda$	$C_D\overline{N}bH$
Cauchy number	$C_Y$	$\rho b{C}_D{H}^3{U}^2/\left(EI\right)$

Equations

$$\nabla \cdot \mathbf{u}$$

$$\rho \left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}\right)$$

$$\mathbf{F}\equiv \sum _{k=1}^N{\mathbf{F}}^k(\mathbf{x},t),$$

$${\mathbf{F}}^k=-\frac{1}{2}\rho b{C}_D{\int }_0^l\delta (\mathbf{x}-{\mathbf{r}}^k)\left(\mathbf{u}-\frac{\partial {\mathbf{r}}^k}{\partial \tau }\right)\cdot {\mathbf{n}}^k\left|\left(\mathbf{u}-\frac{\partial {\mathbf{r}}^k}{\partial \tau }\right)\cdot {\mathbf{n}}^k\right|{\mathbf{n}}^k\mathrm{d}s.$$

$$\frac{\partial {\mathbf{T}}^k}{\partial s}+{\mathbf{q}}^k$$

$${\mathbf{M}}^k=EI\frac{\partial {\mathbf{r}}^k}{\partial s}\times \frac{{\partial }^2{\mathbf{r}}^k}{\partial s^2},$$

$$EI\frac{{\partial }^2{\theta }_k}{\partial s^2}+T_1^k\cos {\theta }_k-T_3^k\sin {\theta }_k$$

$${\mathbf{q}}^k$$

$${\overline{\mathbf{F}}}_R(\mathbf{x},t)=\frac{1}{\text{upi}{R}^2}{\iint }_{C_R(x,y;z)}\mathbf{F}({\mathbf{x}}^{\prime },t)\mathrm{d}{x}^{\prime }\mathrm{d}{y}^{\prime },$$

$$C_R(x,y;z)=\{(x^{\prime },y^{\prime },z):(x-x^{\prime }{)}^2+(y-y^{\prime }{)}^2\le R\}$$

$${\overline{\mathbf{F}}}_R(\mathbf{x},t)=\frac{1}{\text{upi}{R}^2}{\iint }_{C_R}\left[\sum _{k=1}^N{\int }_0^l{Q}_R^k({\mathbf{x}}^{\prime },t)\delta ({\mathbf{x}}^{\prime }-{\mathbf{r}}^k){\mathbf{n}}^k\mathrm{d}s\right]\mathrm{d}{x}^{\prime }\mathrm{d}{y}^{\prime },$$

$$Q_R^k({\mathbf{x}}^{\prime },t)=-\frac{1}{2}\rho b{C}_D\left(\mathbf{u}({\mathbf{x}}^{\prime },t)-\frac{\partial {\mathbf{r}}^k}{\partial \tau }\right)\cdot {\mathbf{n}}^k\left|\left(\mathbf{u}({\mathbf{x}}^{\prime },t)-\frac{\partial {\mathbf{r}}^k}{\partial \tau }\right)\cdot {\mathbf{n}}^k\right|.$$

$$\delta ({\mathbf{x}}^{\prime }-{\mathbf{r}}^k)=\delta (x^{\prime }-r_1^k)\delta (y^{\prime }-r_2^k)\delta (z-r_3^k).$$

$$\frac{\mathrm{d}}{\mathrm{d}s}=\frac{\mathrm{d}z}{\mathrm{d}s}\frac{\mathrm{d}}{\mathrm{d}z}=\cos {\theta }_k\frac{\mathrm{d}}{\mathrm{d}z}.$$

$${\overline{\mathbf{F}}}_R(\mathbf{x},t)$$

$$\overline{\mathbf{F}}=\underset{R\to 0}{\lim }{\overline{\mathbf{F}}}_R$$

$$[x]=[z]=[s]=H,[\mathbf{u}]=U,[t]=\frac{H}{U},[p]=\rho U^2,[\mathbf{T}]=\rho b{C}_{D}H{U}^2.$$

$$\nabla \cdot \mathbf{u}$$

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}$$

$$\frac{\partial \mathbf{T}}{\partial s}$$

$$\frac{{\partial }^2\theta }{\partial s^2}$$

$$w{|}_{z=0,1}=0,$$

$${\left[\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right]}_{z=0,1}=0,{\left[\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right]}_{z=0,1}=0,$$

$$\displaystyle\quad\mathbf{T}\rvert_{s=h}=\mathbf{0},\ \theta\rvert_{s=0}=0,\ \frac{\partial\theta}{\partial s}\biggr{\rvert}_{s=h}=0,\$$

$$z_0$$

$$\mathbf{u}=u(z){\hat{\mathbf{e}}}_x.$$

$${\text{Rey}}^{\ast }=\frac{UH}{{\nu }^{\ast }},$$

$$\frac{{\mathrm{d}}^{2}u}{\mathrm{d}{z}^2}$$

$$\frac{\mathrm{d}{T}_1}{\mathrm{d}s}$$

$$\frac{\mathrm{d}{T}_3}{\mathrm{d}s}$$

$$\frac{{\mathrm{d}}^2\theta }{\mathrm{d}{s}^2}$$

$$\frac{\mathrm{d}x}{\mathrm{d}z}=\tan \theta .$$

$$\frac{{\mathrm{d}}^2\tilde{u}}{\mathrm{d}{z}^2}\sim \{\begin{array}{ll}-\tilde{P}+\frac{1}{2},& \mathrm{if}z\le h,\\ -\tilde{P},& \mathrm{if}z>h.\end{array}$$

$$u=1+{\text{Rey}}^{\ast }\lambda \tilde{u}\sim \{\begin{array}{ll}1-\frac{{\text{Rey}}^{\ast }\lambda }{4}(1-h)(h-z^2),& \mathrm{if}z\le h,\\ 1-\frac{{\text{Rey}}^{\ast }\lambda }{4}h{\left(z-1\right)}^2,& \mathrm{if}z>h,\end{array}$$

$$\frac{{\mathrm{d}}^{2}u}{\mathrm{d}{\hat{z}}^2}=-{\text{Rey}}^{\ast }P+\frac{{\text{Rey}}^{\ast }\lambda }{2}u|u|,$$

$$u_{\text{below}}(z)\sim {\left(\frac{2P}{\lambda }\right)}^{1/2}\text{for 0≤z≤h and h−z=O(1)}.$$

Full text

1 Introduction

The study of fluid-structure interactions with vegetation has a wide range of industrial and environmental applications, including flood control, environmental conservation, and energy production. However, there are a number of challenges in modelling such flows, in particular due to the fact that the vegetation both affects and is affected by the flow. The focus of this work is to develop compact mathematical models that describe flows through submerged vegetated regions and their resultant instabilities. In particular, we predict the critical parameters in different regimes for instabilities which resemble monami—the synchronous waving of vegetation.

1.1 Flow through aquatic vegetation

Climate change is increasing the frequency and severity of hydrological disasters. With the rise of more challenging flow management problems, there is an increasing demand for new solutions. While there are many infrastructural solutions on flow management, there is also an emerging interest in utilising aquatic vegetation, as it is part of the natural habitat helping to sustain our ecosystems. Above all, compared to artificial measures, aquatic vegetation has the promising abilities of adapting to the local environment and self-repairing after destructive events (Morris et al., 2018). On the other hand, flows through vegetation are challenging to model due to interactions between the two. As such, the efficiency of aquatic vegetation in protecting coastal regions and the physical mechanisms involved are yet to be fully understood mathematically (Marion et al., 2014).

Flows through aquatic vegetation are an interesting class of flows to study due to the geometric and mechanical properties of aquatic vegetation. For plants to avoid mechanical failure or being uprooted, they have evolved to be typically flexible and streamlined so that they can passively reconfigure to reduce the fluid load (Vogel, 1994). In addition to their deformable nature, they can also have complex geometries. Their components can have length scales that differ by multiple orders of magnitude. Finally, distinct from terrestrial flows, submerged vegetation has a typical height comparable in magnitude to the water depth in order to photosynthesise (Marion et al., 2014). As a result, a significant proportion of the flow is obstructed by the canopy—a community of vegetation.

There are numerous modelling challenges in capturing the macro-scale and micro-scale properties of the system, and this primarily relates to the feedback mechanism between flow and vegetation. In a complete dynamic model, the fluid will apply a load on each vegetative structure, which causes a resultant deformation and this, in turn, must affect the flow. Thus in general, it appears that the fluid must be solved simultaneously with the configurations of each structure, and it is not clear what reductions can be applied. These challenges associated with modelling flow through vegetation have demanded sophisticated studies; see e.g. experimental works by Dunn et al. (1996); Ghisalberti & Nepf (2004); Hu et al. (2014); Mandel et al. (2019) and numerical works by Mattis (2013); Sundin & Bagheri (2019); Sharma & García-Mayoral (2019).

The question that we seek to answer in this work is whether there exist simpler mathematical models that are able to capture a number of the key physical features obtained in more complicated formulations or experiments. In regards to the development of a simplified model, the vast majority of previous work has fallen in two categories: either models of flow over a specified set of rigid obstacles [see e.g. Ghisalberti & Nepf (2004)]; or models where deformation can occur, but only under a known imposed flow [see e.g. Luhar & Nepf (2011)]. There have been fewer models that emphasise the coupling between deformation and flow. We highlight, in particular, the work of Alben et al. (2002) on flow past a single elastic strip in 2D and Dupont et al. (2010) on flow past an array of rigid and straight elements that are free to tilt. A summary of previous work, and their key features is presented in table 1.

1.2 Instabilities in aquatic vegetation

Part of the emerging interest in instabilities of flow through a canopy is sparked by a fascinating phenomenon known as monami (see the schematic in figure 1)—the progressive, synchronous oscillation of aquatic vegetation (Nepf, 2012). Honami, its counterpart in terrestrial flows is readily observable in daily life when wind blows across a patch of grass or a crop field (de Langre, 2008). Only recently has this phenomenon been explained by Raupach et al. (1996), pointing out that such instabilities arise due to a shearing mechanism that resembles the Kelvin-Helmholtz instability on free mixing layers at the top of the canopy. Once the shear exceeds a threshold, waves develop in the flow which evolve into vortices over time (cf. insets of figure 1). The argument that such instabilities are distinct from boundary layer instabilities is supported by comparing statistics of turbulent kinetic energy in experiments (Finnigan, 2000; Poggi et al., 2004).

Since this explanation has been accepted, many studies have attempted to explain and predict monami mathematically by analysing the stability of steady states in their respective model (cf. Table 1). In particular, Singh et al. (2016) established the dependence between viscous effects and flow instabilities by analysing flows through an array of rigid beams.

1.3 Objectives of this paper

This article is dedicated to analysing the mechanical aspects of flow through submerged vegetation and its resultant instabilities. We develop a coupled model for the fluid flow and the mechanical deformation of the canopy, where the plants are allowed to have large angles of deflection. Using this model, which also accounts for viscous effects, we assess criteria and mechanisms for the onset of instability. Furthermore, determine for which flow regimes de-coupling between flow and vegetation is a reasonable approximation, which is currently unclear.

The structure of this article is as follows. In §2, we derive a mathematical model which couples the dynamics of the flow with the reconfiguration of the canopy, modelled as an array of elastic beams. Using this coupled model, we analyse flows which are steady and unidirectional in §3. In §§4–6, we assess the temporal stability of such steady configurations. The analysis attempts to predict the critical parameters for instabilities which resemble monami and highlight how these parameters differ when the beams are rigid and vertical. We summarise our findings in §7 and discuss limitations and future work in §8.

2 Mathematical model

Consider a three-dimensional domain with fluid contained between $0\leq z\leq H(x,y,t)$. The bed, at $z=0$, is covered by a fully-submerged vegetative canopy that consists of $N$ identical plants of length $l$, which shall later be specified as linear elastic beams. A schematic of the setup and the coordinate system is given in figure 2.

2.1 Equations for the fluid

In a more complete model, the location of each physical plant surface must be calculated as part of the problem. However, provided that the cross-sectional length scale of the plant is much smaller than the separation of neighbouring plants, we can consider the far-field approximation of the drag due to individual plants using a distribution of point forces in the Navier-Stokes equations. We thus consider an incompressible fluid in three dimensions with velocity $\mathbf{u}=(u,v,w)$ at time $t$ satisfying

[TABLE]

where $\rho$ is the density of water, $p$ is pressure, $g$ is acceleration due to gravity and $\mu$ is the dynamic viscosity of water. The additional sink term, $\mathbf{F}$, that appears in the momentum equation (2) incorporates the contribution of the $N$ plants; the precise form of $\mathbf{F}$ is discussed below.

We model each of the $N$ plants as an inextensible linearly elastic beam with diameter $b$, which undergoes pure bending in the $xz$-plane. We parametrise the centreline of the $k$-th beam by $\mathbf{x}=\mathbf{r}^{k}(s,\tau)$ where $s$ is arc length and $\tau=t$ is time. We use $\tau$ rather than $t$ to emphassise the implicit change of variables from $(\mathbf{x},t)$ to $(k,s,\tau)$; this is a Lagrangian description of the canopy. We shall return to this later in §4.1 where we discuss the relationship between Eulerian and Lagrangian frames.

We assume that the collective drag term can be written as

[TABLE]

with each individual term, $\mathbf{F}^{k}$, accounting for the drag on the $k$-th plant. We suppose $\mathbf{F}^{k}$ depends on the difference between the fluid velocity and the beam velocity $\partial\mathbf{r}^{k}/\partial\tau$. In particular, we take

[TABLE]

Here, $\delta$ denotes the Dirac delta function, $C_{D}$ is the drag coefficient and $\mathbf{n}^{k}$ is the upstream normal of the $k$-th beam’s centreline in the $xz$-plane (cf. figure 2); if $\theta$ is the local angle of deflection at $s$, as measured from the upward vertical, then $\mathbf{n}^{k}=(-\cos\theta_{k},0,\sin\theta_{k})$.

Our drag law (4) can be interpreted as follows: for each element of arc length $\mathrm{d}s$, the drag $\mathrm{d}\mathbf{F}^{k}$ is that a tiled cylinder of length $\mathrm{d}s$ would experience (Sumer & Fredsøe, 2006). It has been proven both experimentally and numerically that this drag law is accurate until the beam approaches a configuration parallel to the flow (Ramberg, 1983; Vakil & Green, 2009; Zhao et al., 2009). Analogous expressions for drag have been used by Luhar & Nepf (2016) for steady-state flow and by Leclercq & de Langre (2018) for unsteady flow. In the case of rigid and vertical beams ($\theta\equiv 0$), (4) reduces to the formulation used in Singh et al. (2016). We refer readers to Zhou et al. (2010) for a more extensive review of drag on tilted cylinders.

2.2 Equations for the vegetation

Each individual plant is modelled as a linearly elastic beam with one end clamped perpendicularly to the substrate and the other end left free. For each beam, let $\mathbf{T}^{k}=(T_{1}^{k},T_{2}^{k},T_{3}^{k})$ be the internal stress and $\mathbf{M}^{k}$ be the moment on a cross-section given by position vector $\mathbf{r}^{k}$. By considering the momentum balance on this cross-section which has infinitesimal thickness, $\mathbf{T}^{k}$ satisfies (Landau et al., 1960)

[TABLE]

where $\rho_{B}$ is the density per unit volume of the beam and $\mathbf{q}^{k}$ is the external force per unit length on the beam.

We also consider the angular momentum balance on this cross-section. Noting that $\partial\mathbf{r}^{k}/{\partial s}=(\sin\theta_{k},0,\cos\theta_{k})$ is the local tangent vector of the plant, if we assume that the standard constitutive relation between $\mathbf{M}^{k}$ and $\mathbf{r}^{k}$ holds under dynamic conditions, namely (Landau et al., 1960)

[TABLE]

then $\theta$ satisfies (Howell et al., 2008)

[TABLE]

In this equation, $E$ is the Young’s modulus and $I=\upi b^{4}/64$ is the moment of inertia of the cross-section about the $y$–axis.

Finally, we must determine the load on the $k$-th plant, $\mathbf{q}^{k}$, in our fluid-structure interaction problem. By neglecting the effect of buoyancy in this analysis (by assuming that $\rho_{B}=\rho$), the total load is solely due to the drag in (4), so that

[TABLE]

where $\Omega$ is the fluid domain.

2.3 Homogenisation of the canopy

In real-world scenarios, the canopy is a medium with a complex micro-structure. Even with the simplifications we have made so far, it is impractical to monitor the positions and effects of individual plants in the flow. Instead, we use a simpler averaged model in which the canopy is an effective medium which contributes a bulk volumetric drag term (see e.g. Nepf (2012) and the references therein for previous models in this direction). Here, we briefly present the derivation of such a model by volume averaging.

Consider a fixed $\mathbf{x}=(x,y,z)$. We define $\overline{\mathbf{F}}_{R}$ as the local average of the collective sink over a disk of radius $R$, namely

[TABLE]

where $C_{R}(x,y;z)$ is the two-dimensional disk of radius $R$ centred at the point $\mathbf{x}$ i.e.

[TABLE]

(see figure 3). By the definition of $\mathbf{F}$ in (3),

[TABLE]

where

[TABLE]

Now, by definition,

[TABLE]

In order to sample plants that pass through $C_{R}$, we rewrite the $s$–integral in (11) with respect to $z$. Assuming that the beam configurations do not overturn,

[TABLE]

Therefore,

[TABLE]

where $\mathbbm{1}_{k}$ is an indicator function with $\mathbbm{1}_{k}(\mathbf{x},t)=1$ if the $k$-th plant passes through $C_{R}$ (and zero otherwise).

We now consider the limit in which $R\rightarrow 0$ (with the neighbouring plant separation also tending to zero so that there are still many plants crossing $C_{R}$). Since we have assumed that the plants are identical in the canopy, we approximate $\mathbf{u}$ to be uniform in $C_{R}$ and similarly for the configuration and motion of the plants. With reference to the local averaging argument presented by Chapman (1995), we deduce the continuum approximation of $\mathbf{F}$ as

[TABLE]

where $\bar{N}(x,y)$ is the number of plants planted per unit area (along the bed) and $\mathrm{H}$ is the Heaviside step function. The dependence of $\overline{\mathbf{F}}$ on $\sec\theta$ physically corresponds to the perpendicular distance between neighbouring plants being reduced when they tilt, so that the effective density of the canopy increases. For the remainder of this work, we will consider uniform canopies where $\overline{N}$ is constant.

2.4 Simplification and non-dimensionalisation

Before continuing, we highlight two key simplifying assumptions that we will make. Following the approach by Dupont et al. (2010) and Singh et al. (2016), we neglect the free surface that exists in reality and instead, treat the water depth $H$ to be constant and impose a shear-free condition along the top of the domain. Secondly, we impose a shear-free condition along the substrate by focusing on regimes with the Reynolds number based on a single element (plant) being much larger than unity (Dunn et al., 1996; Ghisalberti, 2002). We shall revisit some of these aspects in the final discussion in §8.

We non-dimensionalise using the variables in our problem with the following scales:

[TABLE]

Foreseeing the calculations ahead in our stability analysis, we take $U$ to be the velocity at the top of the domain for steady unidirectional flows, which we will specify in §3. For the beam equations, we will neglect the inertial terms in the beam equations, (5) and (7), by focusing on the regimes in which $\pi b/(4C_{D}H)\ll 1$ and $\rho U^{2}/E\ll 1$. The former inequality is a good approximation for slender vegetation with $b/H\ll 1$. The latter inequality holds whenever the velocity of the fluid is small compared to the speed of sound through the beams, which is also typical for aquatic vegetation (Lei & Nepf, 2016).

With all the variables being henceforth dimensionless, we have following system of dimensionless equations:

[TABLE]

where the dimensionless parameters $\Rey$, $Fr$, $h$, $\lambda$, and $C_{Y}$ are defined in table 2. In particular, $\lambda$ is the product of the dimensionless planting density, $\bar{N}H^{2}$, the aspect ratio, $b/H$, and the geometry factor, $C_{D}$. Finally, $C_{Y}$ characterises the balance between static deflection and loading due to drag. We have chosen the convention where we incorporate the geometry of the plant by scaling $C_{Y}$ with $C_{D}$ and $I$ instead of introducing a slenderness parameter (de Langre, 2008).

We have seen in (18) that it is natural to write the Navier-Stokes equations in Eulerian coordinates but it is more natural to write the beam equations in body-fitted coordinates (i.e. arc length $s$). These systems are related via $\partial\mathbf{r}/\partial s=(\sin\theta,0,\cos\theta)$. We examine the translation between the two coordinate systems more fully in §4.1.

We presented a complicated system of coupled non-linear partial differential equations for the model of the fluid-vegetation system. The purpose of this work is to analyse the instability mechanism of monami. Therefore, in the next section, we will develop tractable unidirectional reductions of the governing equations so that steady configurations are readily solvable.

3 Steady unidirectional flow

In this section, we seek solutions of the governing system (18) where the flow is steady and unidirectional along the $x$–axis (the streamwise direction), with

[TABLE]

The velocity only depends on the distance from the substrate. Historically, field studies, controlled experiments, and more recently, numerical simulations typically measure $u(z)$ [cf. Nepf (2012)]. To account for eddies in the flow, we follow the approach by Singh et al. (2016) and assign a constant eddy viscosity, $\nu^{*}$, to replace $\mu/\rho$ in (2). Thus, for the remainder of this work, we use the effective Reynolds number,

[TABLE]

to characterise unidirectional flows. In particular, the characteristic velocity, $U$, is chosen such that $u(1)=1$ (see discussion in §2.4).

3.1 Theoretical formulation

Since the flow is unidirectional, we can deduce from the momentum equation of the fluid (2) that the flow is driven by a constant pressure gradient, $P>0$. The governing equations (18) reduce to the following system of differential equations:

[TABLE]

[TABLE]

with $P$ chosen such that $u(1)=1$. To solve for the steady configuration, we rewrite the equations (21) as a system of ordinary differential equations in $z$ and solve the equations numerically [with the corresponding boundary conditions in (18)] in MATLAB using the boundary value problem solver bvp4c. Once we obtain the solution, the centreline of the homogenised plant configuration, $\mathbf{r}$, can be determined in Cartesian coordinates from the relation

[TABLE]

As an aside, we note that our model predicts a parabolic flow profile above the canopy rather than logarithmic as in classic boundary layer flows (Nikora, 2010). Since the canopy enhances flow mixing above the canopy, it has been experimentally shown that the logarithmic scaling is only recovered when $z\gtrapprox 2h$ [see Sharma & García-Mayoral (2018) and references therein]. Provided that $h=O(1)$, which is typical for aquatic vegetation (Nepf, 2012), we can assume that we are not in the logarithmic regime. For flows with $h\ll 1$, the transition zone where $h\leq z\lessapprox 2h$ is known as the roughness sublayer (Finnigan, 2000).

3.2 Unidirectional flows through rigid canopies ($C_{Y}=0$)

In the case of flow through a rigid canopy, the momentum equation of the fluid naturally decouples from the beam equations at leading-order. Below, we present the asymptotic analysis in the limits of sparse and dense canopies.

3.2.1 Asymptotic approximation of $u$ for sparse canopies ($\lambda\ll 1$)

To begin, observe the velocity profiles of figure 4, shown at different values of the canopy density parameter. In the limit the canopy density tends to zero, the velocity profile approaches that corresponding to a uniform and unobstructed flow. Noting that a more convenient perturbation parameter is $\Rey^{*}\lambda$, we let $u=1+\Rey^{*}\lambda\tilde{u}$ and $P=\lambda\tilde{P}$. Then by (21a), we have

[TABLE]

By solving for $\tilde{u}$ (and $\tilde{P}$), we deduce that

[TABLE]

and this confirms the behaviour illustrated in figure 4. In particular, we remark that as the canopy density increases, the flow velocity reduces everywhere due to drag; however, as expected, the reduction is greatest within the canopy itself.

3.2.2 Asymptotic approximation of $u$ for dense canopies $(\lambda\gg 1$)

We also observe from figure 4 that in the dense-canopy limit, $\lambda\gg 1$, the velocity, $u(z)$, is apparently divided into two outer regions ($z<h$ and $z>h$), as well as an inner region near $z=h$. In the limit $\lambda\rightarrow\infty$, we observe that the flow becomes approximately uniform in the canopy, where $0\leq z\leq h$ and $h-z=O(1)$. From (21a), we have that in this region,

[TABLE]

and thus to leading order, the pressure gradient balances the drag, and the velocity below the canopy satisfies

[TABLE]

This matches the result of Poggi et al. (2004, §5) and Singh et al. (2016, §3). Before deriving the solution in the boundary layer, we note that the solution for $z\geq h$ can be found exactly: integrating (21a) and applying the surface boundary conditions, $u(1)=1$ and $u^{\prime}(1)=0$, we find that for the solution above the canopy,

[TABLE]

Note that the constant quantity, $P$, itself must be expanded in $\lambda$ and will be determined through boundary conditions. In the inner region, we substitute $z=h-\lambda^{-1/3}\eta$ and $u=\lambda^{-1/3}\tilde{u}$ into (21a) with $\eta\geq 0$. The inner solution satisfies

[TABLE]

where the integration constant $C=0$ in order to match (26) to this order of approximation. As a result, taking the relevant branch of (28) and requiring that the solution is continuous with (27), we have that for the solution within the boundary layer,

[TABLE]

which is valid for $z\leq h$ and $h-z=O(\lambda^{-1/3})$. Finally, in order to determine the leading-order behaviour of the constant pressure gradient, $P$, we require that the gradient of (27) matches that given by the boundary layer of (28) at $z=h$. This yields

[TABLE]

Note that in the above expression, we have retained the first two orders in $P$ so as to ensure higher accuracy in the above-canopy solution for a wider range of $\lambda$ values. In figure 4, we observe good agreement between the exact numerical flow profiles and our matched-asymptotic approximations given in (26), (27), and (29).

3.3 Unidirectional flows through flexible canopies ($C_{Y}\neq 0$)

Having gained some intuition on flows through rigid canopies, we now explore the differences in flows through flexible canopies—the main motivation of this work. Before we continue, we first define how we will vary the flexibility. We recall from §2 that the Cauchy number, $C_{Y}$, characterises the amount of deflection due to drag. Hence, it varies with the velocity scale (cf. table 2). In order to vary flexibility independently, we vary the ratio

[TABLE]

in our analysis for the remainder of this work. Flow profiles for varying $C_{Y}/\Rey^{*2}$ and fixed canopy density are shown in figure 5.

As for flows through rigid canopies, every flow profile in figure 5 increases monotonically in $z$ and inflects at the top of the canopy. As we increase the flexibility of the vegetation (by increasing $C_{Y}/\Rey^{*2}$), less of the domain is obstructed and we get faster flows at any given $z$. For applications such as flood control, if we use the (dimensional) maximum velocity as a simple measure for damage, the results suggest that upright obstacles attenuate a steady flow most effectively. We will revisit this conclusion in the next section when we consider the unsteady problem.

As an aside, we note that for regimes where the canopy is dense ($\lambda\gg 1$), we expect the approximations in §3.2.2 to hold for plants that are sufficiently stiff. In this limit, since $u\sim\lambda^{-1/3}$, the load and hence the deflection of each plant is negligible.

4 Stability analysis of the steady configuration

Our main task in this section is to extend the stability analysis by Singh et al. (2016) on rigid canopies to flexible canopies. By modelling the canopy as an array of elastic beams, we are interested in the temporal evolution of both the perturbed flow and the canopy configuration.

By considering the steady configurations of §3 as base states, we first derive the system of equations the perturbations satisfy at leading order. We then impose a spectral decomposition on the perturbations and solve the corresponding generalised eigenvalue problem numerically. In this work, we consider two-dimensional disturbances in the $xz$–plane, which is sufficient if we are primarily interested in the critical conditions for instability (Drazin & Reid, 1982).

4.1 Relating Eulerian and Lagrangian frame

Before we consider perturbations of the base flow, we again note that the spatial variables of the fluid are in the Eulerian frame but the variables of the homogenised plant are parametrised by the arc length, $s$. Since it is natural for us to perform stability analysis in Cartesian coordinates, we first express the variables in the system (18) in terms of $(\mathbf{x},t)$.

Firstly, although the transformation between $s$ and $z$ (14) still holds in the dynamic problem, namely

[TABLE]

this transformation is perturbed when $\theta$ is perturbed and solved as part of the problem. Expanding all variables as $f=f^{*}+\epsilon\hat{f}+\dots$, where $f^{*}$ is the steady state and $\epsilon\hat{f}$ is the perturbation, and linearising in $\epsilon$, we find that

[TABLE]

Secondly, the Lagrangian time derivative, $\partial/\partial\tau$, is distinct from the Eulerian time derivative, $\partial/\partial t$. As a result, by (32),

[TABLE]

With the transformations above, we can now rewrite every dependent variable with respect to $(\mathbf{x},t)$. Before we write down the linearised equations, we need to determine the expansions of two key quantities: the velocity of the homogenised plant, $\partial\mathbf{r}/\partial\tau$, and the height of the canopy, $z_{0}$.

4.1.1 Plant velocity

We would like to express the plant velocity $\partial\mathbf{r}/\partial\tau=(\partial r_{1}/\partial\tau,0,\partial r_{3}/\partial\tau)$, at a fixed point on the plant in terms of $(\mathbf{x},t)$. In the Lagrangian frame, for a given $s$,

[TABLE]

It remains to express $\partial\theta/\partial\tau$ as a function of $(\mathbf{x},t)$. By (34),

[TABLE]

Therefore,

[TABLE]

Define $V_{1}=\partial r_{1}/\partial\tau$ and $V_{3}=-\partial r_{3}/\partial\tau$. We have that $V_{1,3}^{*}=0$ for the base state, with

[TABLE]

and

[TABLE]

4.1.2 Height of the canopy

Recall that the height of the canopy, $z_{0}$, also varies as we perturb the base state. This has to be taken into account when we impose boundary conditions on the differential equations of the perturbations at the original height.

Using the expansion of $\theta$ in the integral constraint on $z_{0}$ (18h), gives

[TABLE]

Having derived for expressions of $\hat{V}_{1}$, $\hat{V}_{3}$ and $\hat{z}_{0}$, we now proceed to derive the system of equations the perturbations satisfy.

4.2 Linear stability analysis

By substituting in the perturbed base state into the original system (18) and collecting the linear terms, we find that the perturbations satisfy

[TABLE]

We will give the boundary conditions explicitly in the next section when we consider travelling wave perturbations.

4.3 Travelling wave perturbations

Anticipating the calculations ahead, it is convenient to first define the stream function of the flow, $\psi=\psi^{*}+\epsilon\hat{\psi}+\dots$, such that

[TABLE]

By rewriting $\hat{u}$ and $\hat{w}$ in (41a)–(41c) as derivatives of $\hat{\psi}$, the incompressibility condition (41a) is then naturally satisfied. Furthermore, by equating the derivatives of $\hat{p}$, we can combine (41b)–(41c) into a single differential equation for $\hat{\psi}$. With $\hat{u}$, $\hat{w}$, and $\hat{p}$ being eliminated, we consider a Fourier decomposition in the streamwise direction by letting

[TABLE]

real part understood, with $k$ being the wavenumber of the perturbation along the domain and $\sigma$ being the eigenvalue of the problem. By substituting the ansatz (43) into (41) and use primes (′) to denote derivatives in $z$, we find

[TABLE]

where

[TABLE]

are functions of $z$, $\theta$, and $\phi$. The stream function satisfies a modified Orr-Sommerfeld equation (Drazin & Reid, 1982), modified by additional momentum sinks, $\mathcal{S}_{\mathrm{c}}^{\prime}$ and $\mathcal{S}_{\mathrm{s}}$, in (44a) that are only active in the obstructed part of the domain (at steady state). The expressions for these sinks are coupled to the perturbed beam equations to account for canopy deformation. The corresponding boundary conditions are

[TABLE]

Finally, although $\phi$ and $\phi^{\prime}$ are continuous at $z=z_{0}^{*}$, due to the discontinuous momentum sink in (18b), we have

[TABLE]

Further discussion on how the jump conditions are derived can be found in Appendix A. For a given base state and a given $k$, we seek for values of $\sigma(k)$ such that there are non-trivial eigenmodes of the system of integro-differential equations (44) for $\phi$, $\vartheta$, $\mathcal{T}_{1}$, and $\mathcal{T}_{3}$. In particular, we are interested in the most unstable (or the least stable) mode in each spectrum. We solve this eigenvalue problem numerically.

5 Numerical results

In this section, we visualise some typical solutions for the unstable modes and compare our numerically determined eigenvalues to those of simplified problems with the hope of gaining some physical insights. We present a more extensive exploration of the solution space and discuss the critical parameters for the onset of stability in §6.

5.1 Numerical method for solving the eigenvalue problem

We first rewrite the eigenvalue problem (44) as a system of ordinary differential equations in $z$. In this setting, $\phi$ is defined in $[0,1]$ while other variables are only defined in $[0,z_{0}^{*}]$. Furthermore, $\phi$ has discontinuous derivatives at $z=z_{0}^{*}$. Therefore, we partition the domain into $[0,z_{0}^{*}]$ and $[z_{0}^{*},1]$ and split $\phi$ into two separate functions, namely

[TABLE]

Following the practice of previous work in solving Orr-Sommerfeld problems, we solve this system of equations numerically using a spectral method with Chebyshev polynomials of the second kind (Orszag, 1971). To implement this method, in each interval, we discretise the dependent variables in (44) into their function values at $N$ Chebyshev nodes. Note that $z=z_{0}^{*}$ is an edge node in each interval, which represents $z_{0}^{*-}$ and $z_{0}^{*+}$ when we impose the jump conditions (44m)–(44n). Once we have defined the nodes, we then construct the discrete version of the differential operators in (44) with Chebfun, an open-source package (Driscoll et al., 2014). Finally, we can rearrange the discrete equations so that $\sigma$ satisfies a generalised eigenvalue problem of the form

[TABLE]

where $\mathsfbi{A}$ and $\mathsfbi{B}$ are known matrices and $\mathbf{X}$ is the eigenvector of function values. We solve for all possible pairs of $\mathbf{X}$ and $\sigma$ with eig, the in-built eigenvalue solver in MATLAB.

The classic Orr-Sommerfeld problem is known for its non-normality—its eigenvalues are highly sensitive to perturbations (Reddy et al., 1993). To ensure that the eigenvalues converge, we eliminate rows in the matrix problem (46) that are independent of $\sigma$ (Weideman & Reddy, 2000)—this removes spurious modes which disrupts the convergence of the physical spectrum (Goussis & Pearlstein, 1989). Furthermore, we precondition the problem by rescaling each row with $\|\mathbf{a}_{j}\|_{1}$, where $\mathbf{a}_{j}$ is the $j$-th row of $\mathsfbi{A}$ (Wathen, 2015). We note that without pre-conditioning, the numerical solutions of the eigenvalue problem can exhibit convergence issues between consecutive values of $N$ if $N$ is sufficiently large.

We solve the eigenvalue problem with a starting value of $N=80$ Chebyshev nodes. We increase the number of nodes in intervals of $20$ and recalculate the spectrum until the most unstable eigenvalue is within $0.1\%$ of the previous estimate in the $L_{2}$–norm.

5.2 Typical unstable modes of the eigenvalue problem

Typical results for $\phi$ and $\vartheta$ are shown in figure 6, which correspond to flow perturbation and plant deflection respectively.

We first note from the behaviour of $|\phi|$ that for all cases, the energy of the flow perturbation is localised near the top of the canopy. For slower flows, $|\vartheta|$ increases monotonically in $z$ so that the perturbations to the deflection angle are also largest at the top of the canopy. Moreover, for faster flows, the largest perturbation to the deflection may be in the middle of the canopy. The temporal evolution of such a mode is illustrated in figure 7.

5.3 Temporal evolution of a perturbed steady configuration

We can visualise how the physical configuration evolves in time by perturbing the base flow with a single unstable mode. We plot the streamlines of the flow given by

[TABLE]

where $\gamma$ is the (arbitrarily chosen) initial amplitude.

In figure 7, we present the perturbed configuration at four different instances. The amplitude of the travelling wave grows in time as it convects downstream. If we plot the streamlines in figure 7 without the base flow, we will find closed contours along the top of the canopy corresponding to rolling vortices. Regarding the canopy configuration, plants oscillate synchronously, resembling monami. In particular, we see that when the deflection near the base of the canopy is increased, the canopy becomes more aligned with the flow, reducing the drag higher up the canopy, so that the deflection angle higher up is reduced. We highlight that our flexible canopy model is able to capture the streamwise variation of the canopy height and the local angle of deflection for each plant.

In the remainder of the paper, we are interested in the following aspects: (i) in the limits of small and large canopy density, and small and large flexibility, what is the resultant behaviour of the eigenvalue problem and consequently the flow instabilities; (ii) once the numerical procedure in §5.1 is applied at numerous values of $\Rey^{*}$ and $k$, we can form the neutral stability curve and thus examine the critical conditions in the parameter space for the onset of instability.

We shall begin in §§5.4 and 5.5 with discussion of the asymptotic reductions in the limits of small and large canopy density.

5.4 Stability of flows for sparse canopies

In the limit where the canopy is sparse, in the sense that $\Rey^{*}\lambda,\,\Rey^{*2}\lambda\ll 1$, we know from §3.2 that $u_{0}\sim 1$ and the momentum sink in (44a) becomes negligible. Since the perturbations of the beams are then decoupled from the fluid, we can approximate the eigenvalues of the full problem by solving

[TABLE]

with

[TABLE]

We can solve this reduced problem analytically and deduce that

[TABLE]

for $n\in\mathbb{N}$. The uniform flow is globally stable. We compare the leading eigenvalue with that of the rigid-canopy problem in figure 8a.

In the distinguished limit where $\Rey^{*}\lambda\ll 1$ with $\Rey^{*2}\lambda$ being $O(1)$, we can apply the approximations for the velocity profile in §3.2 for rigid canopies and rewrite the eigenvalue problem as

[TABLE]

with $\phi$ satisfying the boundary conditions (49), $\tilde{u}$ given by (24), and $\sigma_{\mathrm{new}}=\Rey^{*}(\sigma+k\mathrm{i})$. This scaling determines the critical Reynolds number for instability when the canopy is sparse—we will verify this numerically in §6.

5.5 Stability of flows for dense canopies

Recall from our dense-canopy analysis in §3.2.2 that in the limit $\lambda\to\infty$, the flow inside the canopy satisfies $u^{*}=O(\lambda^{-1/2})$ from (26). Moreover, the flow above the canopy is parabolic and given by (27) and (30). Thus with this approximation, we should recover the eigenvalues of the classic Orr-Sommerfeld problem for a plane Poiseuille flow between $[h,2-h]$ with $u_{0}=0$ at $z=h$ (the channel bottom) and $u_{0}=1$ at $z=1$ (the flow centreline). This gives [compare with Drazin & Reid 1982]

[TABLE]

Recall from Chapman (2002) that the above problem has two types of modes: symmetric and antisymmetric. Since we have imposed a no-shear condition (44j) at $z=1$, eigenmodes of the full problem (44) in this limit will only correspond to the antisymmetric modes of the classic problem (52), as seen in figure 8b. These modes are always stable (Orszag, 1971). Therefore, our problem is also linearly stable in the dense-canopy limit.

6 Critical conditions for the onset of instability

We have thus seen that the system is linearly stable when the canopy is absent or infinitely dense. However, we found that unstable modes exist for intermediate canopy densities. Therefore, it is natural for us to investigate the critical conditions for the onset of instability.

6.1 Evolution of the neutral curve

We show in figure 9 contour plots of $\Real(\sigma)$ for various points in the $(C_{Y}/Re^{*2}$, $\lambda)$-parameter space. Firstly, for a given canopy density and flexibility, the base states first become unstable for a critical wavenumber $k$ of $O(1)$. Secondly, although not plotted here, all of the unstable modes have $\Imag(\sigma)<0$, with $|\Imag(\sigma)|$ increasing non-linearly with $k$ for any given $\Rey^{*}$, corresponding to the instabilities being convective downstream and dispersive. This is in agreement with the predictions by (Singh et al., 2016; Luminari et al., 2016) on flows through rigid canopies.

We comment now on a key feature of figure 9 that will be further discussed in §§6.2 and 6.4. Notice that for a fixed flexibility, $C_{Y}/\Rey^{*2}$, and decreasing density, $\lambda$, i.e. as we transition from (b) to (c) to (d), the topology of the neutral curve may change: the area enclosed by the neutral curve tends to zero and the system becomes globally stable. This trend also occurs for fixed $\lambda$ and for increasing $C_{Y}/\Rey^{*2}$ i.e. as we transition from (e) to (f) to (g). Thus, there exists a critical curve,

[TABLE]

in the $(C_{Y}/\Rey^{*2},\lambda)$-plane that separates the two possibilities—globally stable for all $\Rey^{*}$, or unstable for some $\Rey^{*}$. This bifurcation is plotted as a line in figure 9. In particular, our numerical results suggest that $\lambda_{c}\sim C_{Y}/\Rey^{*2}$.

6.2 Behaviour of the critical Reynolds number

We notice from the contour plots in figure 9 that if there are any unstable modes in the parameter space, then for every neutral curve, there exists a critical value for $\Rey^{*}$, $\Rey_{\mathrm{crit}}^{*}>0$, such that flows will be first unstable if $\Rey^{*}>\Rey_{\mathrm{crit}}^{*}$. Therefore, it is natural for us to analyse the behaviour of $\Rey_{\mathrm{crit}}^{*}$ as a function of $\lambda$ and $C_{Y}$, shown in figure 10.

For fixed flexiblity, $C_{Y}/\Rey^{*2}$, we see that the critical Reynolds number, $\Rey_{\mathrm{crit}}^{*}$, exhibits a minimum as the density, $\lambda$, changes. This can be understood intuitively through the asymptotic analysis of §5.4 and §5.5. In particular, as the canopy density decreases to zero or increases to infinity, the base flow becomes globally stable. We will revisit this statement when we analyse the critical shear in §6.3.

The critical values $\lambda=\lambda_{c}(C_{Y}/\Rey^{*2})$ below which the flow is globally stable are indicated by the vertical dotted lines in figure 10.

6.3 Behaviour of the critical shear

In the previous discussion, we have sought to understand the effects of plant flexibity and density on the critical value of $\Rey^{*}$ separating unstable and stable flows. However, there are interesting insights once this is viewed in regards to the shear.

We examine the critical shear in our model for the onset of instability. Suppose we let $\chi$ be the (dimensionless) Navier slip length at steady state at $z=z_{0}^{*}$, namely

[TABLE]

For every flow that is marginally stable, we can then quantify the critical shear with $\chi_{\mathrm{crit}}^{-1}$, the reciprocal of the Navier slip length.

We observe from figure 10b that for any given canopy density:

1. (i)

The critical shear, $\chi_{\mathrm{crit}}^{-1}$, increases monotonically with $\lambda$. 2. (ii)

The critical shear, $\chi_{\mathrm{crit}}^{-1}$, is nearly identical across different plant flexibilities (cf. the inset of figure 10b) even when $\Rey_{\mathrm{crit}}^{*}$ and the configuration of the canopies are distinct. This is highlighted by the inset of figure 10a.

The collapsed data suggests that shear at the top of the canopy is the relevant criterion in determining the stability of steady unidirectional flows. This conclusion is consistent with previous statistical analysis on turbulent flows through vegetation, which suggested that canopy-scale instabilities are shear induced (Raupach et al., 1996; Finnigan, 2000).

For sparse rigid canopies, we found that stability depended on $\Rey^{*2}\lambda$, implying $\Rey_{\mathrm{crit}}^{*}\sim\lambda^{-1/2}$ as $\lambda\rightarrow 0$. Since $\chi\sim(\Rey^{*}\lambda)^{-1}$, this implies $\chi_{\mathrm{crit}}\sim\lambda^{-1/2}$ as $\lambda\rightarrow 0$. The scaling laws that we predicted in §5.4 for the limit $\Rey^{*}\lambda\ll 1$ with $\Rey^{*2}\lambda$ being $O(1)$ are indicated by the dashed lines in figure 10.

Finally, we interpret the U-shaped curve in figure 10a physically. Denoting by $\Lambda$ the canopy density corresponding to the minimum $\Rey_{\mathrm{crit}}^{*}$, we find that for canopies with $\lambda>\Lambda$, the shear threshold is greater and therefore a faster flow is required for instability. For $\lambda<\Lambda$, the shear threshold is lower but because of the reduced drag, a faster flow is required to generate this shear at the top of the canopy.

6.4 What happens if the perturbation of the canopy is ignored?

Despite the complexity of the solutions presented above, it appears that the overall qualitative behaviour of the neutral stability curves shown in figure 9 can be understood through largely hydrodynamical effects. In particular, we may attempt to decouple the canopy and the flow by considering only perturbations in the flow. This is equivalent to setting the perturbations on the canopy, $\vartheta$, and stresses, $\mathcal{T}_{1,3}$, to zero in (43) and solving the single equation (44a) for $\phi$. Once this is done, the neutral stability curves for this decoupled problem are derived, and these are shown in figure 9 as dashed lines.

When the canopy is dense or the plants are stiff, the difference between the two problems are negligible, as discussed in §3.3. However, as we decrease the canopy density or increase the plant flexibility, the modes of the fluid problem are always more stable than the coupled problem, indicating that the movement of the canopy contributes towards instability.

7 Conclusions

We have shown that the mechanics and instabilities of flows over a fully submerged vegetative region can be studied using relatively simple models. In particular, we have presented a continuum approximation that effectively homogenises the vegetation, and identifies five key dimensionless parameters: the Reynolds number, Froude number, submergence ratio, canopy density, and the Cauchy number. The study of two-dimensional steady unidirectional flows allows for the derivation of a number of interesting results, and asymptotic analysis can be used to study particular cases of rigid or flexible canopies, and dense or sparse vegetation.

Temporal stability of steady unidirectional flows can be studied by numerically solving an eigenvalue problem for perturbations in the form of travelling waves. In both limits where the canopy density, $\lambda\rightarrow 0$ and $\lambda\rightarrow\infty$, we have demonstrated that the base configurations are stable. However, for intermediate canopy densities, there exists a critical Reynolds number and critical shear (determined by the flow at the top of the canopy) that separates unstable and stable flows. These critical thresholds vary as functions of the canopy density and flexibility. In particular, we have shown that plant deflection encourages the temporal growth of perturbations. Our flexible canopy model can also capture the temporal evolution of monami by studying the growth of a single unstable eigenmode. In contrast to typical terrestrial flows, the local angle of deflection may not increase monotonically in height.

8 Discussion

One important conclusion from the present study is that the shear at the top of the canopy is a dominant factor in determining the stability of a flow; moreover, we have shown how this result can be ascertained using relatively compact models of vegetative flows. Thus one interesting question that follows is whether there may be particular flow regimes where more involved models are required, and for which the central mechanism for instability may be different. As a particular example, we highlight the review by Nepf (2012), who notes that as the canopy becomes increasingly sparse, the flow transitions from a mixing-layer-like flow to a boundary-layer flow (see e.g. figure 10 as $\lambda\to 0$). In our current model we have neglected the effects of bed shear, and this hints at the need for a model for which the sparse-canopy limit can be more accurately captured.

Another important avenue for progress is the consideration of the particular turbulence model used, and again, there are signs that a more complete model is required. Note that we analysed flow stability using a single mixing length at the top of the canopy (54) and our constant eddy viscosity (20) closure model can be interpreted as effectively averaging the various length scales involved. However, as we have highlighted in the asymptotic analysis of §3, in limits of small or large elasticity/density, it is unclear whether it may be necessary to consider a variable eddy viscosity in order to capture disparate length scales in the flow.

Which length scales might be involved in a more complete model? Firstly, the bed can be rough in reality and the boundary layer along the bed may contribute one length scale. Secondly, there is an element length scale inside the canopy, which is determined by the wakes of individual plants. Finally, for deeply submerged canopies ($h\ll 1$), turbulence that is generated far below is expected to be negligible far above ($z/h\gtrapprox 2$); thus there is a length scale determined by the decay of the canopy-scale vortices. A more detailed description of each scale can be found in Marion et al. (2014), and we envision that a more refined model will impose a closure model that captures a number of these length scales. However, in order to implement such a scheme, we require a better understanding of the importance of these length scales, and moreover, we would hope to determine these scales as part of the model, rather than verified a posteriori from given parameters [see discussion in Poggi et al. (2004)].

In addition to the typical vegetative fluid-structure questions we have highlighted above, there is a great deal of interest as well in the study of intrinsically time-dependent scenarios which include the case of oscillatory flows through vegetation. Compared to unidirectional flows of the same magnitude, oscillatory flows can have higher in-canopy flow velocities (Lowe et al., 2005) and they can potentially create more shear than flows over a bare bed (Zhang et al., 2018). For these situations, vegetation can move out of phase with the waves (Gijón Mancheño, 2016). Moreover, drag can be affected by the wave frequency (Mattis et al., 2019) and the presence of a current (Zeller et al., 2014; Luhar & Nepf, 2016; Lei & Nepf, 2019). A natural extension of our work is to utilise our coupled model in §2 to analyse the interaction between oscillatory flows and vegetation.

Acknowledgement

This publication is based on work partially supported by the EPSRC Centre For Doctoral Training in Industrially Focused Mathematical Modelling (EP/L015803/1) in collaboration with HR Wallingford and US Army Engineer Research and Development Center. We would like to thank A. Dimakopoulos from HR Wallingford and C. Kees from US Army ERDC for their important contributions. We are also grateful to I. Hewitt and D. Vella for insightful discussions.

Appendix A Jump conditions for the Orr-Sommerfeld equation

In this appendix, we discuss the derivation of the jump conditions of $\phi$ at $z=z_{0}^{*}$ (44m)–(44n) presented in §4.3. By equating the derivatives of $\hat{p}$ in (41b)–(41c), and consider travelling wave perturbations (43), $\phi$ satisfies

[TABLE]

with

[TABLE]

The expressions for $\mathcal{A}$, $\mathcal{B}$, and $\mathcal{C}$ are given by (44e)–(44g) in the main text. The modified Orr-Sommerfeld equation (55) here can be rewritten into equation (44a) by splitting the equation into cases above and below $z=z_{0}^{*}$. The corresponding jump conditions can be derived from integrating (55) across the interval $[z_{0}^{*}-d,z_{0}^{*}+d]$, with $d\rightarrow 0$.