Numerical Investigation into the effects of rainfall and long stem plants spacing on Root Water Uptake (RWU)

. Abstract. Woody plants on earthen slopes are a bioengineering solution for the prevention of shallow landslides that occur mostly during a wet season. From a soil-hydrological point of view, slope stability is influenced by plant roots reducing soil water content through transpiration. Despite this, conventional engineering practice tends to ignore the effects of Root Water Uptake (RWU), in part due to the complexity of soil-vegetation-atmosphere interactions. This paper investigates the hydrological effects of plants, which involved seepage simulations performed on two different soil types. Each soil was exposed to different rainfall intensities, and the influence of plants over time was seen from the RWU over time for different configurations of plant spacing and canopy densities. This information with in-situ rainfall data, is useful to assess the effectiveness of plants for slope stability. Further, the relative importance of different mechanisms acting in soil-plant-atmosphere interactions was seen in the RWU data. Although the conducted simulations refer to a horizontal soil profile, the results are useful in more complex geometries such as earthen slopes and may help the design of bioengineering solutions (woody plants) and slope stability assessment. Future research is aimed to investigate additional soil-vegetation-atmosphere mechanisms and additional model geometries and plant species.


Introduction
Soil bioengineering is primarily used in the context of slope stability for water infiltration and erosion reduction, usually by introducing living organisms (e.g., vegetative cover plants and woody plants) or a combination of living and inorganic material (e.g., vegetated riprap). Woody plants affect (positively and negatively) shallow soil mass movement through a variety of mechanisms. The water absorbed through plant roots (transpiration) increases soil suction and strength. However, the Root Water Uptake (RWU) required for transpiration could potentially increase the hydraulic gradient between the root bulb and the surrounding soil, causing an increase in the rate of water transport with possible adverse effects on slope stability. The rate of water transport near plants can also be reduced, which is beneficial for slope stability, due to plant roots reducing permeability by blocking soil pores. Conversely, plant roots can also create preferential water flow paths that increase surface infiltration [1].
For each mechanism affecting water transport rate discussed, the individual qualitative effects are clear. However, when considering all mechanisms simultaneously the overall effects on suction distribution in soil (and consequently slope stability) are dubious due to the interactions between them. Additionally, it is worth mentioning that plants RWU is highly dependent on a variety of factors ( [2], [3], [4] among others), including soil water retention * Corresponding author: andrew.vidler@uon.edu.au characteristics, atmospheric conditions and root bulb architecture, with most of these factors changing over time, further increasing the complexity of the problem to estimate the overall effects.
This study aims to add to the body of knowledge quantifying the overall effect of plants on slope stability and presents an approach to assess the usefulness of plants under rainfall events of different intensity. This is done by investigating the influence of RWU from woody plants on the suction distribution in different soil types, within a numerical model of a horizontal 2D soil profile.
Such a soil profile was modelled so that the expected behaviour can be clarified in a simple situation before the complexity of a slope is added. However, changes in the suction distribution for horizontal soil profiles can still indicate consequences on slope stability.
The study quantitatively assesses to what extent plant RWU can reduce the effects of rainfall events on a shallow soil layer and investigates the interacting effects of plants on moisture movement (velocity of the saturation front and moisture distribution). The consequences of various plant spacing and plant canopy coverage (which affect prevalence of RWU over surface evaporation) were also investigated. To perform the desired study, a native Australian evergreen plant, Melaleuca Styphelioides, was selected for numerical modelling of water transport in two different soils (a silty clay and a silty sand), subjected to a variety of different rainfall intensities.

Methodology
The numerical modelling was performed using the Finite Element commercial software Hydrus 2D by Pc-Progress, which can simulate transient seepage in variably saturated porous media, taking into consideration plant water uptake.
Hydrus also provides parameters for a variety of soil types, of which a relatively fine grained and coarse grained soil were used in this study, the soils used are named silty clay and silty sand within the software. Subsection 2.1 details model parameters that are common to all models presented herein, and subsection 2.2 presents the modelling program along with any model specific parameters.

Common model parameters
A 2D soil profile 16m long and 3m deep, with root bulbs spaced at different intervals, has been modelled. Figure  1 summarizes the boundary conditions adopted in the numerical model: • A no flux boundary condition (B.C.) on the left and right edges of the domain as it is assumed that the model is a small segment of a larger soil profile. • A free drainage B.C. for the lower boundary (i.e., zero pressure head gradient) to represent the situation where the water table is located far below the spatial domain. • An atmospheric boundary is assumed at the soil surface that experiences inflow and outflow by precipitation and evapotranspiration. The initial soil conditions are defined by setting the pressure head in the soil profile to vary linearly with depth, from -2 m at the surface to -4 m at a depth of 3 m. The assumed initial conditions are representative of a relatively wet state of the soil during a rainy season. The atmospheric conditions of the model were defined by autumn conditions of the Australian temperate zone, using climate statistics of Kangaroo Valley (NSW), and were kept constant during the simulation period of 22 days. The conditions adopted include: a relative humidity of 75%, a wind speed of 3.25 ms -1 , a mean daily solar exposure of 15.9 MJm 2 and an average, minimum and maximum temperature of 18.9; 15.6 and 25.10 °C, respectively. The Potential Daily Reference Evapo-Transpiration contribution (ETo) has been computed according to the Penman-Monteith method in the last revision by Allen et al. (1994) [5]. Beer's law was used to partition ETo (=0.15 mmh -1 ) into potential evaporation (Ep) and potential transpiration (Tp), using a radiation extinction coefficient of 0.8 and a Leaf Area Index (LAI) value of either 3 or 8. Hydrus 2D solves the Richards equation for a variably saturated flow, incorporating a sink term S [T -1 ] that accounts for RWU [6]. The 2D formulation of the Richards' equation is a function of the volumetric water content θ [L 3 /L 3 ], the hydraulic conductivity function K [LT -1 ] and the sink term S. θ and K are in turn dependent on water pressure head h [L] and defined in this work, according to the van Genuchten-Mualem model. The hydraulic and retention parameters adopted in the model are shown in Table 1 for each soil. Table 1. Parameters of the van Genuchten-Mualem model adopted to describe the soil water retention curve of the silty clay and the silty sand soils. Hydrus2D estimates the water uptake of roots by relating the potential transpiration with soil restriction to water flow, using a water stress response function α(h). The formulation of Feddes et al (1978) [7] has been adopted to determine α(h) using the following parameters: • h1 (= 0.1m) an arbitrary anaerobiosis point, i.e., when h < h1 the water uptake is assumed to be zero. • h2 (= 0.25 m) and h3 (= 4 m) define a pressure range (h2 ≤ h ≤ h3) where plant roots can absorb water at the maximum rate, i.e., actual uptake is equal to the potential uptake. • h4 (= 160 m) corresponds to the wilting point pressure head, i.e., when h > h4 the water uptake is assumed to be zero. The potential uptake is distributed over the root zone according to a water uptake distribution function b(x,z) [L -2 ] which depends on root architecture and root length density. The parameters of the b(x,z) function adopted in the numerical model were calibrated using experimental measurements of the sandy clay soil surrounding a 2year-old Melaleuca Styphelioides plant with LAI of 1.3. The plant was placed in a large-scale apparatus subjected to different boundary conditions, and the calibration was based on suction records in proximity of the root bulb. For more details refer to [8]. The parameters of b(x,z) were kept constant over time, as the simulation period is short and significant changes in the root bulb are unlikely. The size of the root bulbs used in the models (i.e, ~1.4 and ~1m in the vertical and horizontal dimension, respectively) are three times the size of the root bulb calibrated in the large-scale experiment (Figure 2), which was done to more closely represent the root bulb of an adult plant.

Series 1
The purpose of this series of modelling is to observe the decay in RWU due to soil wetting from rainfall events of different intensities (2.5 to 250 mmh -1 ). For each intensity, the rainfall duration that causes a reduction of the RWU to zero or to a steady value is identified. To this end, the following model parameters and configuration were applied to each soil type: plants spaced at 4 m intervals (as in Figure 1) and a LAI value of 3, corresponding to a potential evaporation and transpiration equal to 1.43×10 -2 and 1.44×10 -1 mmh -1 , respectively.

Series 2
The purpose of this series of modelling is to observe how the spacing interval of plants affects the decay in RWU and suction distribution for each considered soil. The model parameters adopted for each soil are a plant spacing of either 2 m (7 plants), 4 m (4 plants) or 6 m (3 plants), with the LAI value kept constant and equal to 3. A single rainfall intensity of 10 mmh -1 has been imposed on the models.

Series 3
The purpose of this series of modelling is to observe how the influence of the hydraulic gradient between the root bulb and the soil surface affects the RWU in each soil. The models have the plant spacing interval kept constant at 2 m, with the LAI value adopted as either 3 or 8 (in this latter case, potential evaporation and transpiration equal to 2.63×10 -4 and 1.58×10 -1 mmh -1 , respectively). A single rainfall intensity of 10 mmh -1 was imposed on the models.

Series 4
The purpose of this series of modelling is similar to Series 2, however this time a change in the LAI value with plant spacing is considered, which is more realistic. The model parameters used are a plant spacing of 6 m with a LAI of 3 and a plant spacing of 2 m with a LAI of 8. A single rainfall intensity of 10 mmh -1 was imposed on the models.

Series 1
To analyse the results of this modelling series, the RWU of each model is presented in what is herein called a RWU decay curve. This curve represents the evolution of the actual RWU as a ratio of its initial value (at time 0 of the simulation). The decay curve is useful as it indicates when a null RWU is reached, marking when the plant is no longer beneficial in terms of soil strength increment for suction increase. Fig. 3. The actual RWU decay curves for the silty clay (A) and the silty sand (B) models subjected to rainfall events with intensities from 2.5 mmh -1 to 250 mmh -1 .
The RWU decay curves are shown in Figure 3, and all have an S-shape in which three main parts can be identified: • The initial portion of the curve (e.g., in Figure 3A between 0 and 900 min) in which potential and actual RWU coincides. • The central portion of the curve (e.g., in Figure 3A between 900 min and 10000 min) in which an abrupt decrease in actual uptake is observed. • The final portion of the curve (e.g., in Figure 3A for a time greater than 10000 min) in which a flattening in RWU to either a residual or null value is reached. A decay of the actual RWU to a residual value is observed where an equilibrium is reached between the rainwater downward flux and the outward flux by evaporation and transpiration. In this case, the soil volume is unsaturated and anaerobic conditions are not reached (h > h1). On the contrary, a null value of the actual RWU is reached when soil is saturated (or almost saturated) and the restriction of oxygen in the soil prevents transpiration (h ≤ h1).
A residual value of RWU is observed for all rainfall intensities for the fine-grained silty clay ( Figure 3A). whereas only for intensities less than 10 mmh -1 in the coarser silty sand ( Figure 3B). For rainfall intensities above 10 mmh -1 in the silty sand, actual RWU decays to zero.
It is interesting to notice that all the decay curves for the silty clay are coincident ( Figure 3A). This is due to the rainfall intensity being greater than the maximum infiltration rate of the low permeable clay, thus the infiltration rate (and decay of RWU) is equal for all the rainfall intensities. This is confirmed when observing the evolution of infiltration flux over time for each silty clay model, where, for all rainfall intensities, the maximum infiltration flux is reached in a very short time (e.g., for 2.5 mmh -1 rainfall after 13.5 min and for 250 mmh -1 after 0.32 min). For relatively high permeability soils, such as the silty sand, rainfall intensities lower than 150 mmh -1 , do not maximize the infiltration flux. Thus, the RWU decay curves for the silty sand vary with rainfall intensity with a progressive shift towards the right as intensity increases, until an intensity of 150 mmh -1 , after which all the RWU decay curves are coincident ( Figure 3B).
The determination of RWU decay curves for different rainfall events and for a given plant configuration and soil, can possibly show the duration over time in which plants may be useful at removing moisture from soil. The RWU decay curves in Figure 3B show that in a silty sand, plants may continually remove water from soil for rainfall intensities below 25 mmh -1 , and for more intense rainfall events (25 mmh -1 to 250 mmh -1 ), water removal will occur over a period of 3-16 hours. This information could be easily compared with in-situ intensity-duration-frequency curves, which can assess how beneficial plants may be to slope stability during different likelihood rainfall events. For example, if the plant transpiration reaches zero before the end of high likelihood rainfall events, it may indicate that a different plant spacing or other measures should be implemented to guarantee slope stability.

Series 2
In Figure 4, the RWU decay curves of the three different plant spacing intervals for each soil are presented. As the distance between plants decreases (from 6 to 2 m), the RWU decay curve progressively shifts toward the left for the silty clay ( Figure 4A). Similarly, the silty sand shifts to the left but only when plant spacing reduces from 4 m to 2 m. To consistently assess the change in RWU decay with plant spacing, the time to reach 50% of the initial RWU is used as a comparison point. Table 2 shows the time to reach 50% of the initial RWU for each plant spacing and soil, and highlights that the time required to reduce to half the RWU initial value is significantly increased for a smaller spacing, up to 100%. These results indicate that the root bulbs of the plants effectively reduce the velocity of the downward saturation profile, and to a greater extent than the counteracting effect of a higher hydraulic gradient between the soil surface and the root bulbs (from RWU), for the same LAI value.
In Figure 5, the progression of the saturation profile and the suction distribution are shown at two different instants of the simulation for the silty clay and the silty sand (at 8.7 days and 0.4 days, respectively) and for different plant configurations: no plants (images A -D), 6 m and 2 m plant spacing (images B -E and C -F, respectively). For the configurations with a vegetated topsoil, it is clearly observable that the saturation profile curves due to the root bulbs presence, creating an inhomogeneous downward progression. The curved profiles reduce the average depth of the saturation front at each given time for both soils and highlights the predominance of RWU over the increased hydraulic gradient. Table 2. Time (in days) required to reach 50% of the initial RWU for different plant spacings, in each considered soil. The values in brackets are the percentage changes with respect to the configuration with a 6 m plant spacing. A LAI of 3 was adopted in all the considered cases.

Series 3
It is relevant to note that, the adoption of a smaller spacing between plants, generally increases the LAI value and the potential RWU for a given plant. Thus, as the plant spacing reduces, the hydraulic gradient between the root bulb and the surrounding soil increases along with the downward velocity of the saturation profile. To assess this effect, the RWU decay curves and the cumulative surface infiltration were plotted for each soil as shown in Figure 6. The silty clay ( Figure 6A and 6B) shows a RWU decay slightly faster for the model using a LAI value of 8. This is an indication of a greater hydraulic gradient driving a faster water recharge from the surrounding soil. However, the hydraulic gradient does not increase the infiltration rate at the surface, indeed the increment of the LAI appears to reduce the surface infiltration. This behaviour is merely caused by how Hydrus implements water evaporation. In fact, Hydrus can only remove water due to evaporation from the infiltrated water, however this is not accurate when the rainfall intensity exceeds the maximum soil infiltration rate (which is the case for the silty clay) and water ponds and runs off the soil surface. With this implementation of the evaporation, it is consistent (although not accurate with reality) that as the LAI reduces, the rate of evaporation increases and the calculated infiltration rate reduces for the silty clay. When observing the RWU decay curves and the cumulative surface infiltration flux for the silty sand ( Figure 6C and 6D), the results, related to a LAI of 3 and 8, are essentially coincident with one another. Given the results seen in Figure 6 and knowing the implementation of the evaporation contribution in Hydrus, the effects of a greater hydraulic gradient (driven by RWU) appear small on surface infiltration and plant spacing results the primary factor influencing the velocity of the downward moisture movement.    Table 3 shows the times necessary to reach a 50% reduction in RWU (with respect to its initial value) for all the cases considered in Series 4. The increment in time required to reach 50% RWU is ̴ 90% for both soils when plant spacing varies from 6 m to 2 m, which is similar to the 100% increase registered for the same plant configuration but a constant LAI of 3 (see Table  2). From the direct comparison between these values (in days) reported in the last columns of Table 2 and 3, an increase of ̴ 6% in RWU decay rate is observed for the fine-grained soil (5.0 vs 4.73 days) and of ̴ 3% for the coarse grained one (0.30 vs 0.29 days). These results confirm what in Section 3.3 was only supposed: the increase in hydraulic gradient (due to roots uptake) has a small effect compared to plant spacing.

Conclusions
When woody plants are introduced into soil as a bioengineering solution, there are several counteracting hydrological factors that affect soil water transport, which should be considered for a precise assessment of slope stability. The study addresses in particular two of them: changes in soil matric suction through RWU and changes in downward seepage gradient between (wet) topsoil and the root bulb due to root water absorption. The RWU of a plant generally reduces during sustained rainfall exposition and its evolution in time, as a fraction of its initial value, has been herein defined as the "RWU decay curve". When these curves are obtained for a variety of different rainfall intensities, the length of time in which plants, in a certain configuration, result useful for slope stability increment (via RWU) can be assessed.
Results show that the rate of RWU decay decreases as plant spacing reduces, with a 100% decrease moving from a 6 m to a 2 m plant spacing, for the same LAI value. Additionally, the average downward velocity of the saturation profile was clearly shown to reduce with decreasing plant spacing. A reduced plant spacing is associated with a greater plant canopy density and LAI value. When all these elements are taken into consideration, the numerical analyses show that the rate of RWU decay increases due to a greater hydraulic gradient between the root bulb and the surrounding soil; this partly counteracts the effects of plant spacing. However, this effect is relatively minor, as a 6 m plant spacing with a LAI of 8 causes an increase in RWU decay rate of only ̴ 3-6% compared to the condition with the same plant spacing and a lower LAI (=3).
Although the presented results relate to a horizontal soil profile, the general outcomes may be extended to more complex conditions such as earthen soil slopes. In particular, a twofold application is possible: in slope design to guide the choice of the best suited plant species and interspace in function of the in situ expected rainfall conditions; in slope stability problems to assess whether an accurate computation of plant transpiration has to be incorporated in the flow model as a function of the simulated rainfall duration and intensity.
The discussed results introduce several areas of future possible research, such as the incorporation of additional counteracting effects of plants on water transport (e.g., preferential flow paths and pores blockage by roots) and the investigation into the effects of slope geometry and of different plant species.