Model for water infiltration in vegetated soil with preferential flow oriented by plant roots

There is strong experimental evidence that root systems substantially change the hydraulic properties of soil. However, the mechanisms by which they do this remain largely unknown. In this work, we made the hypothesis that a preferential flow of soil moisture occurs in directions which follow the orientation and distribution of roots within the soil, and that this phenomenon alters soil moisture flow patterns. We modified Richards’ equation to incorporate root-oriented preferential flow of soil moisture. Using the finite element method and Bayesian optimisation, we developed a pipeline to calibrate our model with respect to a given root system. When applied to simulated root distributions, our model produced pore-water pressure profiles which agreed with those derived from experimental saturated hydraulic conductivity values of soils vegetated with willow and grass. Agreement improved for simulated root distributions where root segments were oriented in a more realistic way, suggesting that the hydraulic characteristics of vegetated soils are a consequence of root-oriented preferential flow. By incorporating root-oriented preferential flow, our model improves the ability to describe and analyse water infiltration through vegetated soil. This could help optimise irrigation, forecast flood events and plan landslide prevention strategies.


3
Vol:. (1234567890) to vary among plant species (Song et al. 2017). This indicates that the impact of a certain root system on infiltration may depend upon specific traits, for example root distribution and branching angle. However, the precise mechanism by which the soil's hydraulic properties are changed by plant roots is not well characterised.
Several studies have shown the presence of root systems to be associated with a preferential flow (PF) of moisture through soil (Luo et al. 2019;Holden 2005). When using dye to trace the patterns of water infiltration beneath a tropical rainforest, Noguchi et al. (1997) found that dye accumulated along the axes of plant roots, thus indicating a PF that follows the orientation of roots. Noguchi et al. (1999) concluded that this PF is induced by zones of more friable soil, and that the locations of these zones are determined by the position and orientation of roots within the soil. Figure 1 gives an illustration of this phenomenon. A number of experimental studies have found that root activity changes the porosity of the Rhizosphere and other regions of soil close to roots (Dexter 1991;Bruand et al. 1996;Feeney et al. 2006). There are several factors which could contribute to this. One explanation is that the penetration of growing roots fragments and loosens previously compact regions of soil (Angers and Caron 1998). Additionally, it has been shown that mucilage exuded by roots causes the wetting and drying patterns of surrounding soil to change. This is also said to cause increased soil fragmentation (Caron et al. 1992). Furthermore, roots, and associated decaying matter, can serve as a food source for underground fauna, whose activity may also cause fragmentation of the nearby soil (Angers and Caron 1998).
Unfortunately, experimental studies investigating root-oriented PF, and the resultant impact on a soil's hydraulic properties, are lacking (Ghestem et al. 2011). Traditional measures of a vegetated soil's hydraulic conductivity only consider vertical moisture transport, and neglect the contribution of PF induced by individual roots. Electrical Resistivity Tomography allows access to 3-dimensional moisture distributions, providing richer information about the impact of root systems on soil hydraulics (Beff et al. 2013). However, such techniques require specialised equipment and, at present, access is limited. Modelling and simulation based investigations of PF often involve dual-porosity models (Gerke and Van Genuchten 1993). These models comprise of an equation for PF occurring in a network of macropores, coupled with an equation for moisture transport in the bulk soil. They do not, however, explicitly account for the architecture of the macropore network, or the fact that this may be determined by the distribution and orientation of roots within the soil (Noguchi et al. 1999). With these points in mind, the 3 objectives of this study are: (i) To develop the PF model for soil moisture transport that explicitly incorporates root-oriented PF, under the assumption that representa- Vol.: (0123456789) tive volumes of the domain are composed of a volume of bulk soil and a volume of soil in which root-oriented PF occurs, see Fig. 1c. (ii) Define and demonstrate how the PF model is calibrated using experimental data. (iii) Investigate the capacity of the PF model to recreate experimental results, hence providing support to the hypothesis that PF is a contributing factor to water infiltration in vegetated soils.
We integrated the PF model into a software pipeline, which linked root system simulation, volumetric root density construction, finite-element approximation of partial differential equations, and model calibration, in order to simulate 3-dimensional porewater pressure profiles within vegetated soils. Using this pipeline, with the experimental data of Leung et al. (2018), we provided evidence that root-oriented PF is a contributing factor to the hydraulic properties of a vegetated soil. It should be noted that in this work we do not include root water uptake in the model we develop. This is because the experimental data we used to calibrate our model also did not consider uptake. We intend to incorporate root water uptake into future versions of the PF model.

Materials and methods
The PF model Richards' equation, Richards (1931), is the classic equation for moisture transport through fallow unsaturated soil. This equation defines the change with time of soil moisture content (L 3 L − 3 ), where the flux q (LT − 1 ) is induced by a negative gradient in hydraulic potential and given by the Darcy-Buckingham law, Darcy (1856) and Buckingham (1907), Here K s (LT − 1 ) is the saturated hydraulic conductivity of the soil and the variable h (L) is the pressure head, related to pore-water pressure p (P) by the equation p = ρgh, where ρ (ML − 3 ) is the density (1) t (h) + ⋅ = 0 in Ω × (0, T], (2) = −K s (h) (h + x 3 ).
of water and g (LT − 2 ) is gravitational acceleration. The function κ defines the relationship between the soil's hydraulic conductivity and the pressure head h. In Richards' equation (1), the soil domain Ω ⊂ℝ 3 is assumed to have Lipschitz continuous boundary ∂Ω, and T > 0 is a final time. Spatial coordinates are denoted as x = (x 1 ,x 2 ,x 3 ) ∈Ω and x 3 is positive in the upward direction. The constitutive relations in Richards' equation (1) and the Darcy-Buckingham law (2) are the soil moisture retention function and the hydraulic conductivity function κ. The formulas of Van Genuchten (1980) are popular choices in experimental studies of water flow through soil, and analysis of Richards' equation (Radu et al. 2008;Leung et al. 2018). In the moisture retention function (3) the residual and saturated moisture contents are given by r and s (L 3 L − 3 ) respectively, and α (L − 1 ), n > 1 (−) and m = 1 − 1/n are shape parameters. The values of all 5 of these parameters depend upon the soil type under consideration. In the hydraulic conductivity function (4), the parameter l is linked to pore connectivity (−).
We assume the soil domain Ω is such that all representative volumes are composed of a bulk soil volume, uninfluenced by root activity, and a volume of soil closer to roots through which root-oriented PF occurs. Therefore, to formulate a model for soil moisture transport which incorporates root-oriented PF, we propose a decomposition of the macroscopic soil moisture flux into a component for the unoriented flow occurring in bulk soil ( Fig. 1c blue arrows) and the oriented flow occuring in regions close to roots ( Fig. 1c red arrows): Here ψ : Ω → [0,1) is a function describing volumetric root density, which increases and decreases with the level of root abundance, therefore giving a continuous approximation of the distribution of root mass within the soil. The matrix function H : Ω →ℝ 3×3 is termed the "flow-anisotropy" matrix and causes soil moisture flux to be increased in directions parallel to nearby roots. Furthermore, the extent of this oriented preferential flow, induced by H, increases in regions of soil that are closer to roots, and also increases with increasing lateral surface area of nearby roots. This action of the flow-anisotropy matrix H is justified firstly by the fact that the low radial hydraulic conductivity of root tissue forms a barrier to water uptake (Frensch and Steudle 1989;Naseer et al. 2012), thus causing moisture flow to be diverted away from directions perpendicular to the root axis. Secondly, due to root activity altering pore distribution, increasing soil fragmentation and changing wetting and drying patterns, the hydraulic conductivity of soil regions that are close to roots is often enhanced (Feeney et al. 2006;Angers and Caron 1998;Carminati et al. 2010). These 2 factors combine to support the hypothesis, modelled by the action of the flow-anisotropy matrix H, that soil moisture flow is enhanced in directions parallel to the axes of nearby roots. The explicit formulation of both the volumetric density function ψ and the flow-anisotropy matrix H, for a given root system, shall be described in the remainder of this section. However, from this high level description alone, it can be seen from Eq. 5 that the macroscopic flux ̃ is formulated so that in regions of higher root abundance PF is increased in the direction of roots, and in regions of low root abundance soil moisture flux closely resembles the isotropic Darcy-Buckingham flux (2).
To formulate the volumetric root density ψ and the flow-anisotropy matrix H for a given root system ℜ ⊂ Ω , it is assumed that the root system can be discretised into a union of First consider a single root segment ℜ = ℜ 1 ⊂ Ω , as seen in Fig. 2a. The centres and radii of the circular ends of the segment are a 1 , b 1 , and r a 1 , r b 1 respectively. The length, midpoint, volume and lateral surface area of the segment are denoted by 1 , β 1 , V 1 and S 1 respectively. A normalised lateral surface area is defined as S 1 = S 1 /S max , where S max is the maximum possible lateral surface area of any segment for a given root system, and is determined by the criteria of the discretisation process. Using these dimensions, the volumetric root density ψ is formulated so that it is proportional to the probability density function f 1 of the multivariate normal distribution N( 1 ,C 1 ) . The covariance matrix C 1 is determined so that the principal axes of variation are oriented along the root segment, with variances proportional to the length and radius of the segment. This is achieved by setting C 1 = R 1 C 1 R −1 1 , where and R 1 is the rotation matrix such that R 1 (0,0,1) ⊤ produces the unit vector parallel to the segment ℜ 1 . Rotation matrices such as R 1 are computed as in Rodrigues (1840) (supplementary material). To also ensure that the volumetric density function integrates over the domain to give the volume occupied by the root system i.e. ∫ Ω dx = V 1 , the full formulation of the volumetric density function ψ is The volumetric density ψ is therefore a Gaussian function which has its peak at the midpoint β 1 of the root segment ℜ 1 , a spread in the radial directions of the segment that is proportional to the segment's average radius, and a spread in the axial direction of the root segment that is equal to the segment's length. The result of this formulation is illustrated in Fig. 2b and shows how the volumetric root density function ψ gives a continuous approximation of the root segment ℜ 1 in the soil domain Ω.
To define the flow-anisotropy matrix H for the root segment ℜ 1 , we must consider how the flux of soil moisture is affected by this single root segment. With this in mind, the matrix A 1 is defined such that for any vector = ( 1 , 2 , 3 ) ⊤ ∈ ℝ 3 , the operation A 1 ζ increases the magnitude of the vertical component: Here, c a > 1 is defined as the "facilitation" constant and it can be seen that the operation R 1 A 1 (c a )R −1 1 increases the magnitude of ζ in the axial direction of the segment ℜ 1 . The flow-anisotropy matrix associated with the root segment ℜ 1 is therefore computed as the product of the corresponding volumetric root density ψ (7), and the matrix A 1 rotated in the root segment's direction: Vol.: (0123456789) This means that the operation H(c a ;x)q has the effect of modifying the magnitude of the component of the Darcy-Buckingham flux q (2) that is parallel to the axial direction of root segment ℜ 1 . Furthermore, due to multiplication by the volumetric root density ψ, the extent of this modification is greater in regions closer to the root segment ℜ 1 , and also increases for larger values of the facilitation constant c a and segment lateral surface area S 1 . The zone of influence of the individual root segment ℜ 1 , on the orientation and magnitude of soil moisture flow, increases with the length and radii of the segment. This is because the spread of the associated Gaussian density function ψ is determined by the covariance matrix C 1 , the entries of which are proportional to the segment's dimensions. As an example, the macroscopic flux ̃ (c a ;x) = (1 − (x)) + H(c a ;x) for a segment at 45 ∘ from the downward x 3 axis, as seen in Fig. 2a, b, exhibits PF in a direction of 45 ∘ from the downward x 3 axis in regions close to the root segment. However, in regions away from the segment, less PF is exhibited and the macroscopic flux more closely resembles the isotropic Darcy-Buckingham flux. This principle can be extended to complex root systems that are discretised into multiple segments is the total number of segments. Each root segment ℜ i then defines its own volumetric density function ψ i and flow anisotropy matrix H i (c a ; ⋅), and these can be superposed to define global functions for the volumetric root density ψ, and the flow-anisotropy matrix H, corresponding to the entire root system: This method of constructing the volumetric density function ψ and flow-anisotropy matrix H is a type of kernel based method. It has been shown for these methods that increasing the number of kernels used, and decreasing their size, causes the density estimators to converge to the distributions they are designed to estimate (Silverman 2018). The approach described in this work uses kernel functions ψ i that are defined using the geometry of the root segments ℜ i to which they uniquely correspond. Furthermore, the overlapping supports of the kernel functions corresponding to connected root segments ensures that there are no noticeable transitions in density function value at the connection points between root segments. This construction method therefore yields a volumetric density function ψ that gives a physical representation of the whole root system, and integrates over the soil domain to give the volume occupied by the full root system. An illustration of ψ for a real root system is shown in Fig. 2c, d. It can be seen from the definition of the flow-anisotropy matrix for the full root system (10), that H(c a ;x)q aggregates the modifications that are made to the Darcy-Buckingham flux q by each component flow-anisotropy matrix H i (c a ;x). This means the macroscopic flux ̃ in Eq. 5 behaves like the isotropic Darcy-Buckingham flux in regions of low root abundance, but that PF is induced in the axial direction of nearby roots as root abundance increases, see Fig. 2d. As with the single segment root system, the zone of influence of a multi-segment root system depends upon the dimensions of the component segments. The larger the radii and lengths of the segments, the farther reaching their influence on the orientation and magnitude of soil moisture flow. As a result, if there are many large root segments, distributed relatively uniformly throughout the domain, then some degree of root-oriented PF will be induced in most regions of the soil. (10) The PF equation corresponding to a multiple-root system is then formulated as Despite the hypothesis that roots induce PF in directions parallel to their axes, the root tissue itself remains an obstacle to soil moisture transport (Frensch and Steudle 1989). This is incorporated into the PF equation (11) through the factor of η(x) = 1 − ψ(x). In numerical simulations and analysis of the PF equation To formulate the full PF model we equip the PF equation (11) with boundary and initial conditions where j denotes the flux at the boundaries (in the case of the PF model = (x)̃ ) and n denotes the outward unit normal vector. The soil surface is The lateral boundaries of the soil column are Γ 2 = A × L and the base of the soil column is Γ 3 = A × {D} . The Neumann boundary condition q in < 0 on Γ 1 reflects the rate of infiltration into the soil surface and we assume no-flux conditions on the lateral surfaces of the soil domain Ω. The Dirichlet condition on the lower boundary Γ 3 imposes that the base of the soil column is at the interface between the soil and the water table. These choices of conditions mean the model is well disposed to comparison with the experimental data of Leung et al. (2018). However, more general choices of boundary conditions would also be admissible, and could be easily integrated into our calibration pipeline.

Experimental data and the benchmark model
In the work of Leung et al. (2018), effective saturated hydraulic conductivities K * s were obtained for soils vegetated by willow and Festulolium grass root systems, see Table 1.
This allows the formulation of the equation of the benchmark model for the influence of root systems on water infiltration through soil with Here the influence of the root system on the infiltration through Ω is incorporated by way of an experimentally obtained, depth dependent, saturated hydraulic conductivity In this formulation (15), the value x * 3 < 0 is the depth that the root system reaches and K * s > 0 is the experimentally observed saturated hydraulic conductivity of the vegetated soil (Table 1). Equation 13 of the benchmark model does not directly incorporate the influence of root-induced PF on water infiltration. However, the K * s values for the saturated hydraulic conductivity of vegetated soil are experimental observations, and solution profiles can therefore be used as a benchmark to calibrate the value of the facilitation constant c a in the PF model (11), (12). Due to this, we consider (13) of the benchmark model over the same soil domain Ω and assign the same boundary and initial conditions as given in Eq. 12, with the flux given by =̂ .

Formulating the PF model for incomplete root system data
For the willow and grass root systems studied by Leung et al. (2018) The PF model", are unavailable. This means that the volumetric density function ψ and the flow-anisotropy matrix H, cannot be directly computed to formulate the corresponding parametrisation of the PF model (11), (12). However, the total root length, rooting depth, and distribution of root diameters are provided for each root system. As a result, we used a root system simulator (Algorithm 1), to construct distributions of willow and grass root segments which agree with these measurements. Algorithm 1 constructs a list of root segments whose start and end points are uniformly distributed throughout the section of the soil domain that lies above the recorded rooting depth (Leung et al. 2018). For each segment, the azimuthal angle i on the horizontal x 1 − x 2 plane is drawn uniformly from the interval [0,2π]. The polar angle φ i of each segment from the upward x 3 axis is drawn from the interval [ 2 , ] , with i = 2 giving a horizontal root segment and φ i = π giving a root segment pointing vertically downwards. In setting up simulations, we chose either to let all segment polar angles be drawn from a uniform distribution U( 2 , ) or to draw them all from a truncated normal distribution TN( , , 2 , ) . Drawing from a truncated normal distribution means that the segments are more likely to have polar angles φ i from the upward x 3 axis that are close to some mean value ∈ [ 2 , ] . The value of the standard deviation ∈ [0.01, 2 ] determines the extent to which segment polar angles vary from this mean. The third and fourth parameters 2 and π in TN( , , 2 , ) , state the interval from which segment polar angles can be drawn (Burkardt 2014). When the polar angles φ i , of root segments from the upward x 3 axis, are distributed uniformly from horizontal to vertical, i ∼ U( 2 , ) , the resulting root distribution is referred to as uniform. However, when polar angles φ i follow a truncated normal distribution, i ∼ TN( , , 2 , ) , the root distribution is referred to as truncated normal with mean and standard deviation σ.
Formally, Algorithm 1 constructs an array  the diameter of the distal end. Furthermore, if segment polar angles from the x 3 axis are drawn from the truncated normal distribution TN( , , 2 , ) , then the probability P(a ≤ φ ≤ b) that a polar angle φ lies between values a and b, is given by where Here f N and F N are the probability density function and cumulative density function of the parent normal distribution N( , ) respectively (Burkardt 2014). The while loop in lines 12 and 13 of Algorithm 1 states the procedure if the distal end point (M i,5 ,M i,6 ,M i,7 ) of a candidate root segment lies outwith the lateral boundary A of the soil domain Ω. Namely, the azimuthal i and polar φ i angles of the segment are redrawn from their respective distributions to generate a new end point (M i,5 ,M i,6 ,M i,7 ). This step is repeated until (M i,5 ,M i,6 ,M i,7 ) ∈ Ω. In a similar way, the while loop in line 3 of Algorithm 1 ensures that the distance between the experimentally observed rooting depth x * 3 and the deepest distal endpoint in the simulated root distribution is less than some preset acceptable tolerance δ.

Calibration pipeline
To calibrate the value of the facilitation constant c a in the PF model (11), (12) against the experimental data in Table 1, we first set the initial and boundary conditions (12), and values of the other parameters in the PF equation (11), to match the conditions reported in Leung et al. (2018). The soil column Ω = A × L has depth D = − 2m, and circular end surfaces We considered a simulation time of 2 days (d), where the water infiltration rate through the upper surface was g in = −0.01md −1 , with q in = η(x)g in in the PF model (11), (12), and q in = g in in the benchmark model (13), (12). The saturated and residual water contents were s = 0.383 (m 3 m − 3 ), and r = 0.17 (m 3 m − 3 ) respectively, with α = 1.47m − 1 and n = 1.43. The saturated hydraulic conductivity in the PF model (11), (12) and for fallow soil in the benchmark model (13), (12), was K s = 0.187md −1 , and the experimental values for the saturated hydraulic conductivity of vegetated soil K * s (md − 1 ) in the benchmark model (13), (12), were taken from Table 1. The initial pressure head profile in both models was set to hydrostatic equilibrium, with an initial pore-water pressure of − 20kPa at the soil surface Γ 1 , and 0 at the interface between the soil and the water table Γ 3 . The root distribution simulator (Algorithm 1) was then used to generate willow and grass root distributions that agreed with the total root length, rooting depth and diameter distribution recorded by Leung et al. (2018). Willow root distributions after 8 weeks growth were simulated using 25000 segments of the same length. So, with a total willow root length of 47.1m, this meant a segment length of 0.001884m. Each segment in simulated grass root distributions was also given this length, so the total 8 week grass root length of 63.1m meant that grass distributions were composed of 33492 root segments. For both 8 week old grass and willow plants, Leung et al. (2018) state the percentage of total root length which fell in the following root diameter ranges: 0 − 0.0001m, 0.0001 − 0.0002m, 0.0002 − 0.0003m, 0.0005 − 0.001m, 0.001 − 0.002m, 0.002 − 0.003m and 0.003 − 0.005m. This provided empirical cumulative distribution functions for the radii of root segments of both species. Using these distribution functions in Algorithm 1, the interval to which the radius of each segment belonged was randomly assigned, and each radius was drawn uniformly from its corresponding interval. We initially simulated willow and grass root distributions where the polar angles of root segments from the positive vertical axis x 3 are uniformly distributed between horizontal and vertically downwards. The 8 week old willow root system had a smaller total root length than the 8 week old grass (Leung et al. 2018), yet the willow root system reached greater depths within the soil. To reflect this, we also simulated root distributions where the polar angles of root segments from the positive vertical axis x 3 were likely to be closer to π (vertically downwards) for willow and 2 (horizontal) for grass. This was achieved by drawing segment polar angles φ i from truncated normal distributions: TN( w , w , 2 , ) for willow and TN( g , g , 2 , ) for grass. The mean values for the distributions of segment polar angles from the upward x 3 axis were set to w between 3 4 and π for willow, and g between 2 and 3 4 for grass. The values of σ w and σ g determined the standard deviation of segment polar angles from these means, for willow and grass plants respectively. As in the infiltration experiment conducted by Leung et al. (2018), we assumed that these root distributions were static over the 2 day infiltration scenario. The volumetric root density function ψ and the flow-anisotropy matrix H were then constructed for each simulated root distribution, and the PF model (11), (12) was parametrised, in accordance with the experimental conditions. For numerical solutions to the benchmark model (13), (12) and the PF model (11), (12), we used the first-order conformal finite element method with an implicit Euler discretisation in time. The finite element mesh was comprised of tetrahedra, where the minimum circumradius of a component tetrahedron was 0.0145m and the maximum was 0.0275m. The time-step size in the implicit Euler method was 0.01 days. For linearisation of the soil moisture retention function and hydraulic conductivity κ, we used an L-scheme (List and Radu 2016) with constant L = 0.067 (supplementary material). Calibration of the PF model (11), (12) for the optimal value of the facilitation constant c a then required the minimisation of a cost function u(c a ) which evaluates the difference between numerical solutions of the PF model (11), (12) and the benchmark model (13), (12), for a given value of facilitation constant c a , see the next subsection for details. These cost functions cannot be expressed analytically, their derivatives are hard to obtain and their convexity properties are unknown. A classical approach to minimising functions such as these is Powell's method which is a line search algorithm (Brent 2013). Evaluations of the cost functions involve time intensive steps, i.e. numerically solving nonlinear PDEs. Consequently, Powell's method requires long computation times to find minimisers, and a maximum number of function evaluations cannot be set a priori. This motivates the use of Bayesian optimisation which does not require derivatives of the cost function and efficiently explores the parameter space. Furthermore, the number of function evaluations performed can be fixed a priori (see supplementary material for details).
This framework of formulation and optimisation provides a pipeline to calibrate the PF model and generate simulations of moisture transport through soil which incorporate PF oriented by the roots of a root system (Fig. 3). The pipeline takes as input a discretisation of a root system into segments, and the method in the subsection titled "The PF model" is employed to construct the corresponding volumetric root density function ψ and flow-anisotropy matrix H, then parametrise the PF model (11), (12). The evolution of pore-water pressure over time, through the vegetated soil, is simulated by numerically solving the PF model (11), (12) given the hydraulic properties of the soil, an initial pressure distribution and specific infiltration conditions. By employing Bayesian optimisation to minimise the relevant cost function (see the next subsection), the pipeline calibrates the PF model (11), (12) with respect to the saturated hydraulic conductivity of the soil vegetated by the root system, see e.g. Table 1. The parameter values which minimise the cost function are given as the output of the pipeline, along with the numerical solutions to the corresponding calibrated PF model (11), (12).
The finite-element scheme was implemented using the FEniCS library (Alnaes et al. 2015) for Python 3. Applications of Algorithm 1 and the methods for constructing the volumetric root density function ψ and the flow-anisotropy matrix H were carried out using the NumPy and SciPy libraries (Harris et al. 2020;Virtanen et al. 2020). Bayesian optimisation was implemented with the scikit-optimise library for Python 3 (Head et al. 2018), and Powell's method was carried out using the optimize module of the SciPy library. Simulations were visualised using Paraview (Ahrens et al. 2005 We used the PF model (11), (12) to test the following questions: Q1 Is the calibration pipeline capable of identifying (i) the optimal facilitation constant c a when parametrising the PF model (11), (12) for a specific simulated root distribution? Green arrows indicate the pipeline of the PF model (ii) the optimal values of ( , ) with which to simulate a root distribution, and the optimal facilitation constant c a so that the corresponding parametrisation of the PF model (11), (12) reflects experimental results? Here is the mean and is the standard deviation of the truncated normal distribution from which segment polar angles from the upward x 3 axis in a simulated root distribution are drawn.

Q2 When considering simulated root distributions
where the polar angles of root segments from the upward x 3 axis were drawn from a uniform distribution between horizontal and vertically downwards, does the PF model (11), (12)  (ii) For simulated root distributions, where the polar angles i of root segments from the upward x 3 axis were drawn from a truncated normal distribution with a biologically informed mean value , can the PF model be parametrised for both species with a single facilitation constant c a , and still effectively replicate experimental results for water infiltration through vegetated soil?
Q1, Q2 and Q3 are addressed in the first, second and third subsection of "Results" respectively.
To calibrate the PF model (11), (12) with respect to experimental data, we first parametrised the benchmark model (13), (12) using the experimental saturated hydraulic conductivity values K * s in Table 1, and then solved numerically to obtain the pressure head profiles h e j , where the subscript j indicates whether a grass (j = g) or willow (j = w) root system was considered. Given a pressure head profile h e j from the benchmark model, the corresponding laterally averaged pore-water pressure profile p e j (kPa) is To address Q1 (i), we constructed the volumetric root density function ψ for a simulated grass root distribution. Given some value of facilitation constant c a j > 1 , we then constructed the flow-anisotropy matrix H(c a j ;x) and parametrised the PF model (11), (12). Solutions to the PF model (11), (12) were obtained numerically and, using the same formula as in Eq. 16, the laterally averaged pore-water pressure p j,c a j was computed. The optimal value of the facilitation constant c a j > 1 was then found by minimising the cost function where T > 0 is the final time and N D > 0 is a number of equally spaced points {x 1 3 , ..., x N D 3 } from the soil surface (x 3 = 0) to the water table (x 3 = D).
Testing Q1 (ii) involved determining the mean j and standard deviation σ j , of the truncated normal distribution from which the polar angles of root segments from the upward x 3 axis within a simulated root distribution of species j are drawn, and also the facilitation constant c a j with which to parametrise the corresponding PF model (11), (12) so that it most accurately reflected experimental results. For root species j, given a mean segment polar angle j from the upward x 3 axis and a standard deviation of σ j from this mean value, we used Algorithm 1 to construct a root distribution in which segment polar angles are drawn from a truncated normal distribution, and then generated the corresponding volumetric root density function ψ. With facilitation constant c a j the flow-anisotropy matrix H(c a j ;x) was constructed and the PF model (11), (12) was parametrised. The laterally averaged pore-water pressure from the PF model is denoted as p c a j , j , j and optimal parameters (c a j , j , j ) were found by minimising where p e j = g 1000 h e j .
We addressed Q2 by using Algorithm 1 to simulate grass and willow root distributions where the polar angles of root segments from the upward x 3 axis follow a uniform distribution. The cost function u j in Eq. 17 for each uniform root system was then minimised to find the optimal facilitation constant c a j so that simulations from the PF model (11), (12) reflected the experimental results in Table 1. In these root distributions, segments are equally likely to be oriented in any direction. Therefore, solutions do not provide information on the impact of root system structure on infiltration, they only indicate the influence of root abundance.
Question Q3 (i) was tested by using Algorithm 1 to simulate grass and willow root distributions where the orientations of segments are more biologicallyrealistic. The roots in grass root systems are deemed to be oriented more horizontally and those in willow root systems are deemed to be closer to vertically downward. Therefore, we simulated a grass root distribution where segment polar angles from the upward x 3 axis are drawn from a truncated normal distribution with mean g = 0.625 and standard deviation σ g = 1.0, and a willow root distribution where segment angles from the upward axis are drawn from a truncated normal distribution with mean w = and standard deviation σ w = 0.15, see the subsection titled "Formulating the PF model for incomplete root system data". The cost function u j in Eq. 17 corresponding to the parametrisation of the PF model (11), (12) for each of these truncated normal root distributions was minimised to find each optimal facilitation constant c a j , for j = g,w. Since the orientations of segments in these root distributions are drawn from non-uniform distributions, the corresponding parametrisations of the PF model (11), (12) truly incorporate the influence of root-oriented PF.
Question Q3 (ii) was tested by taking the same willow and grass root distributions simulated to test Q3 (i), and considering a single cross-species facilitation constant c a . We then minimised the cost function to find this optimal facilitation constant c a which parametrises the PF model (11), (12) so that its profiles match the pressure profiles from the benchmark model (13), (12) corresponding to both root systems. Since these parametrisations of the PF model use one cross-species facilitation constant c a , any agreement with simulated profiles from the benchmark model (13), (12) must result from incorporating the PF caused by the orientation and position of the segments in the simulated distributions of different species. To provide similar evidence we minimised to determine optimal values for the standard deviations (σ g ,σ w ) of the truncated normal distributions from which segment polar angles from the upward x 3 axis in simulated grass and willow root distributions are drawn, with mean values of g = 0.625 and w = respectively, and also identify the optimal cross-species facilitation constant c a with which to parametrise the PF model (11), (12) so that pressure profiles agree with the profiles from the benchmark model (13), (12).

Effectiveness of the calibration pipeline
First, we tested Q1 (i): whether or not the calibration pipeline was capable of finding a facilitation constant c * a g which optimally parametrised the PF model (11), (12) for a simulated grass root distribution with polar angles from the upward x 3 axis drawn from a truncated normal distribution with mean angle g = 0.625 and standard deviation σ g = 1.0. We minimised cost function u j (17), with j = g, to find the optimal facilitation constant, and the use of Bayesian optimisation provided a moderate advantage over Powell's method (Table 2). Powell's method identified a facilitation constant c *  (Table 2).

3
Vol.: (0123456789) To test Q1 (ii): whether or not the calibration pipeline could identify optimal parameters with which to simulate a root distribution, and also the optimal facilitation constant c * a j with which to parametrise the corresponding PF model, we minimised the cost function û j , given by Eq. 18. This identified the optimal mean value * g and standard deviation * g of segment polar angle from the upward x 3 axis when segment polar angles in a grass root distribution were drawn from a truncated normal distribution, and also the optimal facilitation constant c * a g with which to parametrise the corresponding PF model (11), (12) ( Table 3). Due to the use of random number generators in Algorithm 1, two root distributions that are simulated using the same mean and standard deviation σ, for segment angles from the upward x 3 axis, will never have exactly the same structure. This adds noise to the optimisation process. Furthermore, the increase in the parameter space dimension from 1 to 3, and the nonlinear dependence of solutions to the PF model (11), (12) on ( , ) , makes the existence of multiple local minimisers likely. As a result, Powell's method performed much less efficiently than for the 1-dimensional optimisation problem (Tables 2 and 3). Table 3 shows that Bayesian optimisation identifies parameter values that minimise cost function û j to the same extent as those identified by Powell's method, but in less than a third of the computation time.
In testing Q1 (i) and (ii), the sufficient number of evaluations in the Bayesian optimisation scheme was determined by the value of the cost function at the identified minimiser. If this value was within a tolerance of the cost function value at the minimiser identified by Powell's method, then the number of evaluations was deemed sufficient. Additionally, when searching a 1-dimensional parameter space for an optimal facilitation constant, the Bayesian optimisation scheme would eventually begin to cluster its chosen evaluation points in a small region of the parameter space. These clustered points yielded similarly low cost function values, which suggested that the scheme had converged. This was also used to judge the number of evaluations required (Figure 1 Seeking the minimising facilitation constant c * a g of the cost function u g , see Eq. 17. The minimising facilitation constant c * a g will optimally calibrate the PF model (11), (12) when it is parametrised for a simulated grass root distribution in which root segment polar angles from the upward x 3 axis were drawn from a truncated normal distribution, over the interval [ 2 , ] , with a mean angle of g = 0.625 and a standard deviation of σ g = 1.0. The parameter interval searched was [600, 800] with a starting facilitation constant c a g = 762.33

Method
Bayesian Optimisation Powell's  Table 3 Convergence results for optimisation schemes seeking minimisers of the cost function û g (c a g , g , g ) , see Eq. 18 Minimising the cost function identifies the optimal mean * g and standard deviation * g for the truncated normal distribution from which segment polar angles from the upward x 3 axis are drawn in a simulated grass root distribution, and also the optimal facilitation constant c * a g , so that the corresponding PF model (11)  in supplementary material). Given the results in Tables 2 and 3, it was therefore judged that 25 cost function evaluations was sufficient for the Bayesian optimisation schemes used to obtain the results in Tables 4, 5 and 6, and Figs. 4 and 5.
Both results from this section indicate that our pipeline is effective in efficiently calibrating the PF model (11), (12) using experimental data. In the 3-dimensional parameter space there appear to be multiple local minimisers of the cost function û g (Table 3). However, in the case of searching a 1-dimensional parameter space for an optimal c * a j that minimises a cost function u j , results suggest the existence of a unique global minimiser (Table 2). This is supported by Figure 1 in the supplementary material, which shows the convergence of the Bayesian optimisation scheme that was used here.

An increase in root abundance causes an increase in infiltration rate
We next addressed Q2: whether or not the PF model (11) (12), is capable of replicating experimental results for water infiltration through vegetated soil, when parametrised for simulated root distributions in which the angles of segments from the upward axis are drawn uniformly between 2 (horizontal) and π (vertical). Laterally averaged pore-water pressure profiles from the optimal parametrisations of the PF model (11), (12), for these simulated uniform root distributions, were found to agree strongly with the corresponding profiles from the benchmark model (13), (12) (Figure 4 and columns 2 and 5 of Table 4). Thus supporting the hypothesis that increasing root abundance increases infiltration. Furthermore, from Fig. 4, the minimum porewater pressure is lower in soil vegetated by willow compared to grass, meaning more water infiltrates through soil vegetated by willow.
Drastically different optimal values of facilitation constant c * a g and c * a w were obtained when the PF model (11), (12) was calibrated for willow and grass root Table 5 The value of the cross-species axial facilitation constant c * a which minimises the cost function u, see Eq. 19 This minimising value optimally calibrates the PF model (11), (12) when parametrised for: (i) a simulated grass root distribution in which segment polar angles from the upward x 3 axis were drawn from a truncated normal distribution with mean angle g = 0.625 and standard deviation σ g = 1.0, and (ii) a simulated willow root distribution in which segment polar angles from the upward x 3 axis were drawn from a truncated normal distribution with mean angle g = and standard deviation σ g = 0.15

Species Grass Willow
Parameters in truncated normal distribution for φ i over  Table 6 The optimal parameters (c * a , * g , * w ) found by Bayesian optimisation to minimise the cost function û , see Eq. 20 Here * g is the standard deviation of the truncated normal distribution, with mean g = 0.625 , from which the segment polar angles from the upward x 3 axis were drawn in a simulated grass root system, and * w is the standard deviation of the truncated normal distribution, with mean w = , from which the polar angles from the upward x 3 axis were drawn in a simulated willow root system, and c * a is the facilitation constant, so that the PF model corresponding to either simulated root distribution is optimally calibrated to the experimental data in  distributions in which segment angles from the upward x 3 axis were drawn from uniform distributions. This means that these parametrisations of the PF model (11), (12) do not reflect the impact on infiltration of root-oriented PF, they only show the impact of varying root abundance.

Agreement with experimental infiltration data is improved by incorporating root-oriented PF
The next question to be addressed was Q3 (i): whether or not the PF model (11), (12) effectively replicated experimental results for water infiltration through vegetated soil when parametrised for simulated root distributions in which the polar angles of segments from the upward x 3 axis were drawn from a truncated normal distribution with biologically informed mean. In this case, laterally averaged pore-water pressure profiles obtained from optimal parametrisations of the PF model (11), (12) agreed strongly with the profiles obtained from the benchmark model (13), (12) (Columns 3 and 4 of Table 4). The lower values of the cost function when evaluated at optimal facilitation constant values u j (c * a j ) , indicate that this agreement was stronger than for root distributions in which segment orientations from the upward x 3 axis followed a uniform distribution (Table 4). Furthermore, the optimal facilitation constants c * a g and c * a w were considerably more similar than for root distributions in which segment orientations from the upward x 3 axis were uniformly distributed.
Finally, considering Q3 (ii): whether or not the PF model (11), (12) can effectively replicate experimental results for water infiltration through vegetated soil when parametrised, using a single cross-species facilitation constant c a , for root distributions that Fig. 4 The effect of root abundance on infiltration. (a) (willow), (c) (grass) Laterally averaged pore-water pressure profiles, simulated from the benchmark model (13), (12) and optimal parametrisations of the PF model (11), (12) for simulated root distributions in which segment polar angles from the upward x 3 axis were drawn from a uniform distribution. The values of the facilitation constant c a j that optimally parametrise the PF model (11), (12) for the different root distributions are shown in Table 4 were simulated by using a truncated normal distribution, with biologically informed mean value, to assign segment polar angles from the upward x 3 axis. It was found that pore-water pressure profiles from parametrisations of the PF model (11) (12) for simulated root distributions, where segment polar angles from the upward x 3 axis were drawn from a truncated normal distribution with a biologically informed mean value, agreed more strongly with experimental data than solutions from parametrisations of the PF model (11), (12) for the simulated root distributions in which segment angles from the upward axis were drawn from a uniform distribution (Table 5 and Fig. 5a, c, f and h). Furthermore, pore-water pressure profiles from parametrisations of the PF model (11), (12) for root distributions in which segment polar angles from the upward x 3 axis were drawn from a uniform distribution, and from parametrisations of the PF model for root distributions where segment polar angles were drawn from a truncated normal distribution, differ most at the shallow soil depths, where root density is highest ( Fig. 5e and g).

Discussion
Improving models for soil moisture transport by incorporating root-oriented preferential flow The influence of root systems on soil hydraulics is modelled in many ways. Approaches often involve modifying Richards' equation, for example by introducing a sink term for water uptake that depends on the root biomass distribution (Wu et al. 1999). This oversimplifies root system influence since properties like root diameter distribution and branching angle are neglected.
Richards' equation for moisture transport in the soil, coupled with the Darcy flux for axial moisture transport through root tissue, is often used to solve simultaneously for moisture transport in the soil and in the root system (Arbogast et al. 1993;Doussan et al. 2006). The coupling is achieved through a sink/ source term in each equation which is proportional to the radial conductivity of the root tissue and the difference in matric potential between the soil and root moisture. The influence of root abundance is incorporated through the dependence of the radial hydraulic conductivity of root tissue on the fraction of roots which occupy that layer of soil. However, no dependence of the soil's hydraulic conductivity on root abundance or structure is considered, despite the experimental evidence for this (Song et al. 2017;Leung et al. 2018).
Dual-porosity and dual-permeability models for soil moisture transport are another approach for incorporating flow along preferential channels created by plant roots (Gerke and Van Genuchten 1993;Simunek et al. 2005;Valle et al. 2017). These models either couple different parametrisations of Richards' equation, or couple a transport equation or ordinary differential equation with Richards' equation, in order to account for the flow in the channels and the bulk soil. Dual-porosity or dual-permeability models do not, however, incorporate the influence of root system structure, despite evidence for this determining the layout of preferential channels (Noguchi et al. 1999).
We have addressed the shortcomings of existing approaches with a model for moisture transport in vegetated soil that incorporates root-oriented preferential flow by using a volumetric root density function ψ and a flow-anisotropy matrix H to impose a dependence of the soil moisture flux on root abundance and orientation. Our approach tackles two limitations of existing models. Firstly, the incorporation of root orientation and distribution information into a flow-anisotropy matrix, allows us to account for the influence of a specific root system's architecture on the flow of moisture through the soil which it inhabits. Secondly, since this information is aggregated in the flow-anisotropy matrix H, explicit representation of root systems is not needed, making simulation across a wide range of spatial and temporal scales possible.
Results show that the PF model (11), (12) can be parametrised and calibrated effectively, using a single cross-species facilitation constant c a , for simulated willow and grass root distributions in which root segment orientations are drawn from biologicallyrealistic non-uniform distributions. Simulations from these parametrisations of the PF model agreed more strongly with experimental data (Leung et al. 2018), than those from parametrisations of the PF model that only incorporate the effect of root abundance, and not the root-oriented PF caused by the structural traits of different root architectures (Figs. 4, 5 and Table 4). The capacity to model root-oriented preferential flow appears therefore to improve agreement with Fig. 5 The capacity of the PF model to replicate experimental infiltration results, when incorporating the effect of root system architecture. (a) (willow), (c) (grass) Laterally averaged profiles, simulated from the benchmark model (13), (12) and the PF model (11), (12), where the PF model was optimally parametrised for simulated willow root distributions in which segment polar angles from the upward x 3 axis were drawn from a truncated normal distribution with mean angle w = and standard deviation σ w = 0.15, and grass root distributions simulated using mean polar angle g = 0.625 and standard deviation σ g = 1.0. (b), (d) Laterally averaged profiles of the volumet-ric root densities ψ. (e), (g) Laterally averaged profiles from optimal parametrisations of the PF model (11), (12), the faintest lines correspond to the earliest times and the lines of full opacity show the final profiles. (f), (h) Distance at each time step, � 2 , between profiles from the PF model when optimally parametrised for uniform and truncated normal root distributions, see Eq. 16 for the definitions of p j,c * a and p e j . For all results in this Figure, the PF model (11) (12) was parametrised using the cross-species optimal facilitation constant c * a = 665.78 experimental data on how root systems impact water infiltration through vegetated soil. This provides evidence that root-oriented PF may be a key cause of the variances in hydraulic characteristics that are observed between soils vegetated by different plant species. Due to root activities e.g. elongation and mucilage exudation, the wettability and porosity of the rhizosphere soil in the direct vicinity of roots will differ from that of bulk soil (Carminati et al. 2010;Dexter 1991;Bruand et al. 1996;Feeney et al. 2006). In the PF model, this heterogeneity and anisotropy in the hydraulic properties of a vegetated soil is incorporated using the flow-anisotropy matrix H. Future PF models may be enhanced by encoding some dependence on the volumetric root density into the residual moisture content, saturated moisture content, saturated hydraulic conductivity, and other shape parameters in the Van Genuchten (1980) functions for soil moisture content (3) and hydraulic conductivity (4). In this way, the model would explicitly incorporate the changes to soil hydraulic properties that are observed experimentally in regions very close to root surfaces, whilst still exhibiting the phenomena of root-oriented PF.
Acquiring data with which to parametrise the PF model With the PF model (11), (12), we have developed a pipeline from data on a root system's architecture to simulations of the root-oriented PF it induces (Fig. 3). Assembling this pipeline requires the decomposition of root systems into segments, where the dimensions and position of each segment are available. Such data can be obtained by excavating and digitizing root systems (Danjon et al. 1999). More recently, x-ray computed tomography has allowed 3D architectures of root systems to be imaged without excavation (Zhao et al. 2020), thus providing a potential method of obtaining our required data in a non-invasive manner. Data can also be generated by root growth simulators like CRootBox (Schnepf et al. 2018) and OpenSim-Root (Postma et al. 2017), or by modelling the volumetric root density ψ as the solution to a partial differential equation for root growth (Dupuy et al. 2005). The latter being amenable to applications at different spatial and temporal scales. With data on root system architecture and saturated hydraulic conductivity, our pipeline (Fig. 3) can be applied to complex root systems to simulate the patterns of preferential infiltration that they induce (Fig. 6). Furthermore, if sufficient root system data is provided, then calibrating the PF model (11) (12) requires searching a 1-dimensional parameter space for the facilitation constant c * a that minimises a cost function like that given by Eq. 17. Table 2, and the convergence results in Figure 1 from the supplementary material, suggest that a Bayesian optimisation scheme will locate a unique global minimiser of such a cost function.
Towards new approaches integrating root-oriented preferential flow into soil hydrological models Vegetated soils generally exhibit higher rates of infiltration during rainfall and irrigation than bare soils (Huang et al. 2015;Leung et al. 2018). Below ground however, the impact of root systems varies. Fig. 6 The PF model applied to real root architecture data. a A discretisation of a pine root system into multiple segments (Dupuy et al. 2005). b Volumetric root density function ψ for the discretisation in (a), and the resultant moisture transport from the parametrised PF model (11), (12). Arrow orientation and size indicate flow direction and strength respectively Results obtained by Song et al. (2017) showed that PF caused by Vetivier grass root systems facilitated downward moisture drainage through a landfill site, resulting in harmful chemicals leaching into water supplies. In contrast, Bermuda grass root systems were found to induce a pattern of PF which reduced soil moisture drainage (Song et al. 2017). Furthermore, results of Lourenço et al. (2006) showed that hill slope regions where below ground pore-water pressure is high are more prone to landslide, and that pressure distributions vary significantly with the type of vegetation cover (Ghestem et al. 2011). The idea of harnessing the different PF patterns induced by different vegetation types may therefore have great potential as a method of making soils more resilient to drought, and maintaining desirable conditions for agriculture and land management purposes. Unfortunately, real world implementation of this strategy is limited because it is difficult to predict the effect of vegetation types on the hydraulic properties of soil. Software packages which model flood inundation focus on above ground water flow, with some account taken for the different drying and wetting cycles of affected soil (Teng et al. 2017). The hydraulic properties of cultivated soil affect water distribution and nutrient dynamics, which has consequences for soil irrigation and fertility (van der Heijden et al. 2013). Widely used crop models, such as APSIM (Holzworth et al. 2014) and DSSAT (Jones et al. 2003), predict responses to environmental conditions and integrate groundwater flow models into their simulations. They do not, however, account for PF induced by crop roots. Models like the PF model (11), (12) help to improve the understanding of these PF patterns that are induced by root systems of specific crops. Therefore, they could be used to facilitate the design and implementation of species-specific irrigation schemes, based on the root system architectures of different crops. This work shows that it is possible to incorporate the influence of root architecture while modelling soil hydraulics. Furthermore, results suggest that there may exist a single facilitation constant with which to parametrise the PF model (11), (12), so that it adequately describes the root-oriented PF induced by different root system species at different stages of growth. This indicates the potential of the PF model (11), (12) as a tool to help understand general water transport patterns in vegetated soil, and supplement field scale studies into the impact of specific root systems on water infiltration through soil.

3
Vol:. (1234567890) Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.