Vertical roots are modeled as circular obstacles within a two-dimensional rectangular domain with an optimized irregular mesh. The mesh was created with the mesh generation tool in the iRIC solver suite. Mesh convergence tests were performed by comparing the sizes of hydrological features in the wake region to those previously published in order to assure model accuracy while minimizing computational cost. The mesh generation tool of the iRIC allows mesh refinement regions to be included allowing much finer mesh resolution around each vertical root than regions far from the root array. After convergence tests to within a 7 percent difference with similar models, the maximum area for mesh elements around each vertical root was set to 0.0006m ² for 10 centimeter and 14 centimeter diameter roots and 0.001m ² for roots 20 centimeters in diameter or larger (see Figure 1 panels d and e for a visual example of the two-dimensional model). The maximum area of mesh elements is then gradually increased to 1.3m ² far from the root array. Minimum angle between cell vertices is set to 20 degrees. This results in the total number of elements in each model to range between 11,000 and 18,000 depending on the number of root refinement regions being used.

For most simulations, root diameter was kept constant at 20cm (root diameters are varied between 10cm and 80cm) and are evenly spaced at the vertices of an equilateral hexagonal array. Initial bed elevation was set to ten meters across the entire model area, and initial water depth is set to one meter. Boundary conditions were selected to simulate a typical hydrologic regime that would be experienced by a riparian tree in a floodplain forest at peak flood levels as depicted in Figure 1E ). The blue arrows represent the boundary conditions on each of the four sides of the model. Flood discharge from the left to the right side of the model area was varied between 40m ³/sec to 10m ³/sec. Boundary conditions on all other three sides were set to a constant water height of 1 meter, allowing both sediment and water movement across the boundary. To achieve results applicable to a variety of soil types, we used a mixed grain diameter sediment model with a ‘well-mixed’ sediment profile consisting of particles in the size range from 0.001mm to 34mm and a D ₅₀of 3.2mm. The morphological factor was kept at 1, relative weight of bedload material was set to 1.65, void ratio of bed material was 30 percent, and the exchange layer thickness was 0.5m. Flow was allowed to stabilize for 500 seconds to avoid any effects from initial flood wave motion. Computational time steps ranged from 0.01 seconds to 0.08 seconds depending on the density of the root array and flood time was kept constant at 25000 seconds for a total of 24500 seconds of morphological change. Flood time was determined by allowing the discharge event to continue until the bed elevation reached a steady state.

In MFlow02, water continuity is described by

where t is time, V is the 2-D flow velocity vector, and h is the water depth. The equations describing the momentum of the flow are,

Here, D t ≡ ∂ ∂ t + V · ∇ is the material derivative, D x y = ∂ x ( ν ∂ x ) + ∂ y ( ν ∂ y ) , u and v are the flow velocity components in the x and y direction, ν is the kinematic eddy viscosity, g is the gravitational acceleration, H is the water surface level (depth + ground elevation), f is a Coriolis parameter, ρ is the water density, and τ b x , τ b y are the bottom shear stress component in the x and y direction given by τ b x / ρ = C f u u 2 + v 2 and τ b y / ρ = C f v u 2 + v 2 . C f is the riverbed friction coefficient given by C f = g n 2 / h 1 / 3 . The kinematic eddy viscosity ν , is calculated from the standard k-e model as, ν = C μ κ / ϵ where ϵ = C e κ 3 / 2 / l is the energy dissipation rate, l is the length scale of the turbulence, C μ is a constant equal to 0.09 , κ is the turbulent energy given by, κ = 2.07 u * 2 , where u * is the bottom friction velocity given by, u * = n g ( u 2 + v 2 ) / h 1 / 6 , where, n is the Manning roughness coefficient. and The dimensionless shear stress, τ * used to calculate sediment discharge is calculated as, τ * = u * 2 / s g d , where s is the submerged specific gravity of the suspended sediment particle and d is the diameter of the particle. Bed load sediment discharge is calculated from the Meyer-Peter-Muller formula,

where d is the sediment grain diameter, σ , is the gravel density Meyer-Peter and Muller (1948). τ * c is the critical tractive force calculated from the Iwagaki formula Iwagaki (1956). τ * ′ is calculated from the Kishi and Itakura formula Itakura and Kishi (1980). The total sediment discharge set by the formula above is divided into a sediment discharge in both the normal (n) and tangential (s) direction of the river flow streamline as, q i = q b ( v b i / V b − κ ∂ z ∂ i ) where i stands for either the n or s direction, κ = τ * c / μ s μ k τ * , μ s is the static friction factor, μ k is the kinetic friction factor, z is the height of the river bed, and v b i is v b s and u b n and are the (s) and (n) directional components of the flow velocity near the river bed and are calculated as, v b = 8.5 u * and u b = − 7 v b h / r . V b is the absolute value of the flow velocity near the river bed. Here r is the radius of curvature of the river flow streamline calculated as,

The velocity of buoyancy of suspended sediment is calculated from the Itakura-Kishi formula Itakura and Kishi (1980), E s = K [ α * ρ u * Ω / ( τ * σ ω 0 ) − 1 ] where, K is a constant equal to 0.008, α * is a constant equal to 0.14, and w 0 is the sediementation speed calculated from the Rubey formula. The concentration of suspended sediment at the referential level can be calculated using the Lane-Kalinske formula Lane and Kalinske (1941),

where, C is the concentration of suspended sediment, z α is a referential height equal to 0.05 h . The conservation of mass of the depth-averaged concentration of suspended sediment is described by,

where, D x and D y are the diffusion coefficients of the suspended sediment. The change in elevation, z is then,

Where, q x = q n c o s ( θ ) − q s s i n ( θ ) and q y = q n c o s ( θ ) + q s s i n ( θ ) and θ = tan − 1 ( v / u ) .

A mixed grain diameter model was run, where the accumulation curve of riverbed grains is divided into n hierarchies. A representative grain diameter d_k indicates the existence possibility p_k of a particular representative grain. Central grain diameter is defined as, d m = ∑ k = 1 n p k d k . Additionally, the sheltering effect when calculating the non-dimensional critical tractive force of each grain diameter is calculated as,

The conservation of sediment volume is then expressed as,

First we compare the sediment height ( Figures 2A, C ) and fluid velocity profile ( Figures 2B, D ) around a tree with vertical roots with one without. When a tree is surrounded by vertical roots, a sediment bar of approximately one meter above the surrounding topography forms between one and eight patch lengths downstream of the tree. When the tree has no vertical roots, no downstream bar is formed. The elevation versus time plot in Figure 2 shows the elevation at 1, 4, and 8 patch lengths downstream for the tree with knees. Only a single curve is given for the tree without knees as no spatial variation in sediment height was observed.

The difference in downstream elevation can be explained by the increased drag on the fluid causing a delay in patch-scale turbulence formation due to bleed flow through the cylinder array. A more detailed description of this wake hydrology can be found in Yagci et al. (2016); Chang et al. (2017); Kitsikoudis et al. (2020). This consequence is particularly advantageous for plants whose seeds require deposition onto an exposed surface. As such, vertical roots can be viewed as organs that generate downstream protected microhabitats where the seeds can successfully germinate.

To determine how knee density (number per area, ρ=N/A) affects downstream bar production, we position equally sized knees (20cm in diameter) at equal distances from one another at the vertices of a equilateral hexagonal grid. We vary one of the three variables, ρ , A , N , and observing another, while keeping the third variable constant. In Figure 3 , panels a and d are the total downstream area above three different threshold depths, and the percent velocity change respectively, for trials holding the patch area constant and varying the density. In panels b and e, we plot the same two quantities, but this time, holding the number of roots constant and varying their separation. Similarly, c and f vary the number of roots, while holding root spacing constant. In Figure 3A ), a maximum is observed where an increasing vertical root density stops producing more viable downstream area. This point will be referred to as the optimal root density.

In previous studies, clump-type vegetation has been modeled as regularly spaced cylinders described by the solid volume fraction, ϕ = N ( d / D ) 2 where N is the number of cylinders, d is the diameter of a cylinder and D is the diameter of the circular patch Tanaka and Yagisawa (2010); Nicolle and Eames (2011); Chang and Constantinescu (2015); Chang et al. (2017). For comparison purposes, we also provide the ϕ values of all the arrangements.

As we see in Figure 3 , downstream elevation is best preserved at densities between three and four roots per m 2 ( ϕ between 0.12 - 0.14). Percent velocity change was calculated from the ratio of average fluid velocity directly downstream of the root array to the average fluid velocity far upstream of the root cluster. In Figure 4 we show the contour plots of the downstream elevation and fluid velocity around each of the three root clusters denoted by roman numerals in panel d) of Figure 3 . We show in Figure 3 that as patch density increases, the percent velocity change for the fluid entering and leaving the root array also increases. For densities lower than the optimal density, fluid leaving the root array is slowed, acting as a wake barrier and preventing the re-connection of the faster outer streamlines. This decreased velocity leads to less erosion and increased sedimentation, resulting in an elevated downstream bar. When the optimal density is reached, the root array begins to generate downstream patch-scale turbulence causing a truncation of the sedimentation region and bringing the elevated region nearer to the root array. As seen in panels c and f, the solid volume fraction ϕ is not always a good predictor of downstream area preservation. The diameter of the patch must also be large enough to separate the flow long enough to create a bar.

Another interesting feature of our simulation results depicted in Figure 4 is that, further increasing the root density beyond the optima decreases the distance between the root array and the elevated peak. This suggests that trees with less-dispersive propagules would be selected for higher root density as to draw the bars closer, rather than simply maximizing the downstream area for their seedlings.

Species whose vertical roots grow continuously (such as the baldcypress), would start out with a lower root density, encouraging seedling growth farther away from themselves; however, as the tree ages and its knees grow larger and wider, the downstream bar will move closer to itself, ultimately leading to a seedling replacing its parent.

In Figure 4 , we observe a transition in the velocity contour plots from a slowed steady wake into patch scale turbulence, which occurs when root density is increased. This trend is consistent with the wake structures experimentally observed in Tanaka and Yagisawa (2010); Chen et al. (2012); Zong and Nepf (2012). Furthermore, the resulting increase in downstream elevation observed in Figure 4 agree with the trend empirically shown in Tanaka and Yagisawa (2010), except, in our simulations, we are able to observe that the total viable soil area stops increasing once the vertical root density hits an optimum value.

Nicolle and Eames (2011) observed that patch-scale turbulence occurred at ϕ values between 0.0884 and 0.1451. The optimal root density that we find here ( Figure 3 ) lies in the same regime. Also, our observed pattern in wake characteristics with increasing ϕ agrees qualitatively with the fully three-dimensional large eddy simulations of Chang and Constantinescu (2015). However, our finer sampling of root density as well as our more detailed description of the downstream elevation profile reveals an optimal peak in soil area as well as the dependence of bar position on root density.

The different colored curves in Figures 3A–C , depict the results for when a different threshold is used during analysis. The red curve depicts the amount of area downstream of the root array that is above a threshold of 9.97m, the orange curve for a threshold of 9.84m and the yellow curve for a threshold of 9.7m. As we see, remarkably, the threshold value does not change the location of the optimal point on the curve.

The root arrays in the previous simulations could just as well be thought of as individual vegetation stems, indicating that an optimal vegetation density for maximal downstream propagation success exists. Vegetation with canopy size comparable to its stem diameter or that have shade tolerance might be able to reach such optimal densities. However, vegetation with large canopies must grow at lower densities and thus would need to produce vertical roots to obtain maximal downstream reproductive success. To determine the cutoff of when it is beneficial for vegetation to produce vertical roots, a series of additional simulations were run where smaller diameter knee roots are added into larger arrays of tree trunks. Tree trunks are modeled as either 40 or 80cm diameter obstacles in the computational mesh and the knee roots are modeled as 10 cm diameter obstacles. For each given trunk configuration, knee roots are added and the area above a 9.71 m threshold is calculated. These data can be seen in Figure 5 as the open symbols and dashed lines. For configurations that are already of significant solid volume fraction, any additional knees results in a decrease in downstream area above the threshold. However, configurations of low density see a large benefit in producing knees. This shows that vegetation growing below the optimal solid volume fraction for their diameter will see a benefit in growing knees, and vegetation able to grow at or above this density would see a decreasing benefit.

Our next set of simulations quantify the cooperativity between neighboring trees, as mediated by their vertical roots. In these simulations, two root patches consisting of 37 twenty centimeter diameter roots at their optimal separation were spaced at varying distances from one another.

The root patches were offset so that the edges of the patches aligned. Figure 6 (left) depicts the downstream area above a 9.71 m threshold for these two adjacent root patches as a function of the distance between the two patch centers. Contour plots of the downstream elevation profile for the labeled data points as well as that of a single patch are given in Figure 6 . An increasing cooperative advantage is seen as the patch distance decreases, indicating that trees should seek to close the gap between themselves and their neighbors to reap the maximal cooperative effect. However, a gap up to 3 m will still provide additional cooperative effects. Our theoretical results agree with the experimental flow tank observations, in that the wake interactions between neighboring patches enhance sand bar growth Yamasaki et al. (2021a), and that enhanced downstream deposition is observed for two closely neighboring patches Meire et al. (2014).

In Figure 6 (right), we show how clusters of two and three trees with knees work together to produce considerably more viable area than a single tree. In addition to cooperating with neighboring trees, the vertical roots of an individual tree may also benefit from this cooperativity effect. The inset image in the lower right of Figure 6 shows a knee grouping, which may be a self-cooperating strategy to reduce biomass expenditure while maintaining large sedimentation effects.

Another pattern seen in Figure 6 is that the root patch furthest downstream reaps a greater benefit from the cooperative effect due to sandbar location. Future simulations with more vegetation patches in different arrangement angles could reveal further interesting competitive and cooperative dynamics between multiple root patches.

Figure 4 suggests that a tree can modify its geomorphological influence to suit its life history by altering the density of its vertical roots. By increasing the root density beyond 3.7 roots / m 2 , a tree can draw its downstream bar towards itself, which would be more advantageous for trees whose propagules are adapted for close germination, such as bruguiera gymnorhiza which produces cigar shaped propagules that are often viviparous and germinate after sticking into the mud. In contrast, trees whose propagules more readily float and dispersed hydrochorically such as Xylocarpus moluccensis and Sonneratia alba would perhaps benefit from maximizing the overall area suitable for germination. Long lived species like Taxodium distichium may reap both rewards, where early in their life when their knees are of smaller diameter and the vertical root patch is at a lower solid volume fraction, offspring will germinate farther away from the parent. As the tree ages and the knees become larger and the solid volume fraction increases, the downstream bar will move closer towards itself, where one of its offspring can eventually overtake the parent’s spot in the forest.

Our results on the cooperativity between two root patches emphasizes the advantage of having neighbors who also produce vertical roots. In cypress forests, the fluctuating recruitment rate of seedlings often leads to patches of similar aged trees. After an initial competition for access to sunlight, the surviving neighboring trees turn into allies as their combined footprint of knees enhances the amount of fertile germination grounds for their seeds.

Both the aeration and sedimentation properties of vertical roots could be the reason these appendages were selected for in wetland environments. In mangroves, the presence of lenticels and aerenchyma suggest a closer tie to gas exchange, while in baldcypress the lack of these structures suggests a closer tie to sedimentation. In their 150 million years of existence (Taylor and Taylor, 1993: Seward, 2010) baldcypress trees have endured intense flooding events, like those at the end of the quaternary ice age, possibly because of their vertical roots that allowed them to spread into flood scoured plains when other species could not.

It is interesting to note that their close relative, the pond cypress (Taxodium ascendens) appears to be losing their reliance on the production of knees and also have a tendency to grow to smaller dimensions. By reducing canopy diameter, they are able to grow in large density clusters called ‘cypress domes’. It could be the case that this recent speciation is a response to less severe flood patterns. By growing at greater densities they instead alter geomorphology to their advantage at a group (rather than individual) scale, consistent with our simulations and their interpretation, as discussed in the subsection “Vertical roots are only beneficial for low density vegetation”. It would be interesting to test this idea empirically: Trees with greater canopy spreads would need to grow at a lower densities and therefore would be in a greater need of vertical roots, so if our claims regarding knee function is correct, then cypress knee root production should correlate with the canopy spread of the tree.

The outcomes of this paper are consistent with a number of qualitative and quantitative earlier observations. For example, it was shown earlier that cypress knees lessen the storm runoff capabilities of streams Miroslaw-Swiatek and Amatya (2017) and since baldcypress seeds are most likely to germinate in regions of sediment exposed during flood drawdown Demaree (1932); Titus (1990); Rutherford (2015), an elevated downstream region would act as a prime location for its seeds to germinate. Additionally, baldcypress seedling mortality was tied to the elevation at which they germinate, with seedling death occurring after 60 days of submergence Conner et al. (1986); Souther and Shaffer (2000); Middleton (2000). Here, we showed that a tree with knees is much more effective at producing these downstream protected microsites than a tree without.

The role of vertical mangrove roots on sediment accumulation has also been noted in qualitative empirical observations, particularly in their value in coastline protection Feller et al. (2010); Lee et al. (2014). Interestingly, vertical roots (the pneumatophores) were better at conserving sediment in mangrove forests when compared to stilt roots Krauss et al. (2003). Additionally the amount of oxygen uptake was found to be less in their cone roots than in the stilt roots Kitaya et al. (2002), strengthening our hypothesis.

The bidirectional feedback loops between vegetation coverage and the geomorphodynamics of regions with high hydrodynamic influence is lately gaining increased attention Cornacchia et al. (2019); Larsen (2019); Yamasaki et al. (2021a); Huai et al. (2021) with vegetation found to self-organize itself in environments of hydrologic influence Schwarz et al. (2018); Cornacchia et al. (2020). Vegetation density has been determined to be an important factor in how these ecosystems develop and has led to a range of experimental and computational work on understanding how clump-type vegetation affects hydraulic processes Tanaka and Yagisawa (2010); Chen et al. (2012); Zong and Nepf (2012); Meire et al. (2014). Our findings are consistent with these previous studies and expand upon them by looking more closely at the parameter space which affects downstream sediment bar size and shape.

More recently, flume tank experiments have investigated the drag coefficients of various vegetation arrays Nair et al. (2022), and looked at quantifying the incipient motion of sediment, suggesting an optimal vegetation patch porosity to reduce erosion Kazemi et al. (2021). Because the focus of this study was on the incipient motion of downstream sediment, they gathered data only for approximately one patch length downstream, and used only three sizes of 9 cylinders which provides a porosity resolution too wide to adequately describe the downstream sedimentation characteristics in the flow regime where patch scale turbulence is induced.

Our results suggest that the primary adaptive advantage of vertical roots is to constitute a clump-type vegetative profile to maximize its impact on the fluid flow. This observation deepens our understanding of the vegetation and geomorphology feedback loop. It demonstrates that vegetation does not simply play a passive role in influencing flow patterns, but rather puts energy resources into harnessing geomorphodynamic processes to its own benefit.