Call for Participation: Collaborative Benchmarking of Functional-Structural Root Architecture Models. The Case of Root Water Uptake

Three-dimensional models of root growth, architecture and function are becoming important tools that aid the design of agricultural management schemes and the selection of beneficial root traits. However, while benchmarking is common in many disciplines that use numerical models, such as natural and engineering sciences, functional-structural root architecture models have never been systematically compared. The following reasons might induce disagreement between the simulation results of different models: different representation of root growth, sink term of root water and solute uptake and representation of the rhizosphere. Presently, the extent of discrepancies is unknown, and a framework for quantitatively comparing functional-structural root architecture models is required. We propose, in a first step, to define benchmarking scenarios that test individual components of complex models: root architecture, water flow in soil and water flow in roots. While the latter two will focus mainly on comparing numerical aspects, the root architectural models have to be compared at a conceptual level as they generally differ in process representation. Therefore, defining common inputs that allow recreating reference root systems in all models will be a key challenge. In a second step, benchmarking scenarios for the coupled problems are defined. We expect that the results of step 1 will enable us to better interpret differences found in step 2. This benchmarking will result in a better understanding of the different models and contribute toward improving them. Improved models will allow us to simulate various scenarios with greater confidence and avoid bugs, numerical errors or conceptual misunderstandings. This work will set a standard for future model development.

1. Introduction

A growing number of different modeling techniques and software libraries are now available to build functional-structural root architecture models. Different available models of root architecture and functions have been discussed and qualitatively compared in Dunbabin et al. (2013). The available models differ in the way they represent different processes, such as root growth, water flow, solute transport are captured and translated into mathematical equations (process-level differences); in how they solve mathematical problems by their choice of analytical or numerical approach, numerical scheme, programming technique (solution-level differences); and in how they couple the different processes to the full model (coupling-level differences). However, the extent of discrepancies is currently unknown. Thus, a framework for quantitatively comparing functional-structural root architecture models is required. In addition to the explanatory or predictive power of a model, it is also important to understand the performance of these models, e.g., in terms of accuracy or computational cost. The most commonly used type of functional-structural root architecture models represent the structure of the root system as a 1-dimensional branched network of discrete segments which is geometrically embedded in a 3-dimensional soil domain (Koch et al., 2018b). The root architecture may either be known from measurements, such as 2D or 3D images, or from root architectural models. Suitable models are then used to simulate the “functions,” such as carbon flow and use in root systems (e.g., Bidel et al., 2000), rhizodeposition (Nygren and Perttunen, 2010), competition between species (Dunbabin, 2007), plant anchorage (Dupuy et al., 2007), water and nutrient uptake (Dunbabin et al., 2006; Javaux et al., 2008). Exchange between soil and root is typically modeled via source/sink terms. From the point of view of the soil domain, roots are often considered as line sources, i.e., it is assumed that their diameter is small compared to the relevant spatial scale of the soil. The advantage of this approach is that it does allow to consider root system architecture (position of each segment in time and 3D space) explicitly while being computationally less expensive than an explicit representation of root volumes in the soil domain. By direct comparison with explicit 3D simulations, Daly et al. (2018) showed for the case of young wheat plants that the error made by neglecting root volumes physically present in the soil domain is negligibly small in case of root water uptake. Thus we may expect that, for plants where the line source assumption holds, models of this type are sufficiently accurate. They are also computationally cheaper than explicit 3D and allow the consideration of older and thus larger root systems. The challenge is now to develop a commonly accepted framework for benchmarking functional-structural root architecture models. This includes defining a set of benchmark problems to test model accuracy and performance. We propose that models should be evaluated against two different kinds of references: First, we will develop simple benchmark scenarios, if possible with analytical solutions, that serve as a reference for model verification. Secondly, we define data sets that can be used as references for the evaluation of more complex models without analytical solution. These data sets should as good as possible describe the system we want to model and contain as little uncertainty as possible (Luo et al., 2012). This benchmark activity focuses on two processes, root architecture development and root water uptake. We propose this benchmarking framework to be used by the community of modelers and other participants to compare their model outputs against those of the reference solutions of benchmarks defined in this paper. The use of this framework thus aims to be a collaborative effort. We will refer to any numerical model that implemented some or all of the benchmark problems as “participating model” or “simulator.”

2. Benchmark Problems for Models of Root Architecture and Function

In order to benchmark models of root architecture and function, we propose a multi-step approach with growing level of complexity. The individual benchmarks refer as much as possible to published work, however, we streamlined the different problems and made the notation consistent throughout this paper. A list of symbols is provided in Table 1. The intrinsic nature of functional-structural root architecture models involves multiple coupled domains and processes. A single process in a single domain (e.g., water flow in soil) is referred to as “module” here. The first set of benchmarks (M1–M3) is about individual modules (M) only, i.e., they either deal with only root growth, water flow in soil or water flow inside roots. The scenarios are simple, possibly have analytical solutions, and the goal is to build trust in the accuracy of the individual participating models and to help interpret potentially diverging results of the coupled benchmark problems. Benchmark problems M1 are about root architecture development. It is known that the representation of growth processes can be very different between different simulators. Thus, the goal is to calibrate each simulator individually to given root image data (reference data). M2 is about modeling water flow in soil. Here all participating models solve the same equation, namely the Richards equation, and differences may occur due to differences in numerical implementation. M3 deals with water flow inside the root system for static soil water conditions. As for M2, differences between models are expected to be mainly due to the numerical implementation of this well-defined process. The second set of benchmarks (C1 and C2) is about coupled root-soil models. Benchmark problems C1 consider a static (non-growing) root system and focus on comparison of numerical representation of agreed-upon equations and process representations as well as on the coupling approach to compute the sink term for root water uptake. For this benchmark, we provide a reference solution that is based on a computational mesh that was generated with consideration of the physical presence of the roots in the soil domain. Thus, root water uptake was simulated not by a sink term but as a boundary condition at the root surface in soil. Our approach is similar to Daly et al. (2018) but in addition couples the soil domain to the root domain so that pressure gradients along the roots are simulated. Benchmark problem C2 compares the water uptake of fully coupled models with growing root systems.

TABLE 1

Table 1. List of notations.

Each benchmark problem is described in a Jupyter Notebook that is publicly available on a github repository. Each Jupyter Notebook has a list of contributing authors at its beginning. We will provide codes for automatic analyses and comparison of different model results with the reference solutions or reference data. This makes the analysis transparent and easily modifiable and facilitates including even future participating models' outputs at any later time.

2.1. Levels of Contribution

Any group using or developing functional-structural root architecture models is invited to participate in this collaborative model comparison. Not every model might be suited for all of the provided benchmark problems. Thus, every participant may decide in which individual benchmark problem they would like to participate. However, to reach a certain level of complexity, the “module” benchmarks should be simulated first before the “coupled” benchmarks. Table 2 gives an overview of the key features of these problems and their implementations. One important aim of this activity is a joint publication that shows and discusses the results of the different participating models in comparison to the reference solutions and reference data provided as well as to gain an overview of the extent of deviations between the different simulators.

TABLE 2

Table 2. Description of benchmark scenarios to be implemented in 3D functional-structural root architecture models^a.

2.1.1. How to Participate

The participation includes three steps:

(1) Registration: Any interested researcher is welcome to contact the corresponding author of this paper, Andrea Schnepf, with the following information: Name, affiliation, name or reference to the participating simulator. Upon signing a letter of agreement confirming that results of other participants will not be published without consent, researchers will be accepted as participants and enabled to include their individual simulation results to the github repository of this benchmark initiative, https://github.com/RSA-benchmarks/collaborative-comparison.

(2) Simulation: Each participant implements all or a selected number of benchmark problems in their respective simulator and makes the results in the prescribed formats available to the github repository through pull requests. Requested formats include the Root System Markup Language, RSML (Lobet et al., 2015) for root architectures and the Visualization Toolkit, VTK (Schroeder et al., 2006) for 3D and 1D simulation outputs. Python scripts to read and write RSML files will be provided on the github repository. Packages to read and write VTK files are for example available at https://pypi.org/project/vtk/.

(3) Analysis and publication: The analysis of results and computation of relevant metrics, such as root mean square error, coefficient of determination or Nash–Sutcliffe efficiency, will be done by the code implemented in the Jupyter Notebooks for each benchmark problem. The final goal is to jointly publish the results.

2.2. Benchmarks for Individual Modules

2.2.1. Module M1: Root System Architecture Models

Root system architecture models (RSA models) are that module within a complex functional-structural plant model that simulates the structure, topology, and 3D placement of the roots. They simulate the growth of root systems as (upside down) tree-like structures based on rules regarding elongation, branching and death. Mostly, they are discrete models and represent the root system by a mathematical graph (i.e., nodes and edges/root segments). Each node or segment may be additionally associated with attributes, such as radius, age or hydraulic properties.

The aim of this first benchmarking exercise is to determine if root architecture models currently available are able to reproduce realistic root architectures when being parameterized on the basis of a common experimental data set (Figure 2A). The particular challenge to benchmark RSA models is to include the stochastic nature of these models. We propose to perform the benchmarking of those models in four steps: (1) Parameterizing the root architecture models based on the provided experimental data, (2) Simulating a set of root systems for a dicotyledonous (Lupinus albus) and a monocotyledonous (Zea mays) plant species following two benchmark scenarios (M1.1 and M1.2), (3) Export and store the simulated root systems as Root System Markup Language (RSML) files (Lobet et al., 2015), and (4) Compare the simulation results using the data analysis pipelines available in the associated Jupyter Notebooks. The analysis pipelines are explained below and illustrated in Figure 1. In particular, we include persistent homology as an approach that augments purely trait-based comparisons, i.e., two root systems with the same total root length could be very different based on the persistent homology approach.

FIGURE 1

Figure 1. Example of root images used for the benchmarking dataset. (A) Shows an image of lupin root systems, 11 days old, growing in an aeroponic setup. (B) Shows an image of a maize root system growing on filter paper (5 days old). All images were analyzed using the semi-automated root image analysis software SmartRoot (Lobet et al., 2011), colors distinguish different root orders. The RSML files containing the full information about the root systems are provided on the github repository in the folder “M1.1 RSA calibration\M1.1 Reference data.”

2.2.1.1. M1.1 Root system architecture model calibration

The different available root architecture models (see e.g., Dunbabin et al., 2013) are partly different in the way they represent the growth processes, such that the equations describing these processes are very different. For example, branch emergence is a function of apical zone length in CRootBox (Schnepf et al., 2018b) while it is a function of delay time in Root Typ (Pagès et al., 2004). Root radius is a time-dependent function scaled to distance along the root in OpenSimRoot (Postma et al., 2017) while it is computed according to the pipe model in ArchiSimple (Pagès et al., 2014). Thus, we are looking at process-level differences between the different models, and each participating RSA model will have a different set of parameters that drive root growth. This is the reason why, in this benchmark, we do not prescribe a parameter set as in e.g., M2 or M3, but we let each participating model derive its respective model parameters based on a reference dataset. In this first benchmark (M1.1), modelers simulate root systems for the same duration as the age of the root systems in the reference dataset.

2.2.1.1.1. Reference data set. Although the parameterization of 3D models using a set of parameters derived from 2D images has some limitations, it has been shown to be a simple and efficient strategy allowing the simulation of realistic 3D root systems (Landl et al., 2018). Our reference dataset contains two distinct sets of images: (1) images of lupin roots grown for 11 days in an aeroponic setup (Lobet et al., 2011), and (2) images of maize roots grown for 8 days on filter papers (Hund et al., 2009). All images were analyzed using the semi-automated root image analysis software SmartRoot (Lobet et al., 2011) and root tracings were saved as RSML files for further analysis (Figure 1). These RSML files were then processed using functions of the R package archiDART developed to compute root system- and single root-level metrics (Delory et al., 2016, 2018). These metrics have been made open-access (https://github.com/RSA-benchmarks/collaborative-comparison/tree/master/root_architecture/data) and should help modelers to parameterize their respective RSA model.

2.2.1.1.2. Required output. The following results are to be uploaded via pull requests to this path on the github repository: M1 Root architecture development/M1.1 RSA calibration/M1.1 Numerical results.

1. A text file including the outcome of the calibration step, i.e., the set of model input parameters required for the specific simulator.

2. Simulation output from running the root architecture model using this parameter set in RSML format. Due to the stochastic nature of root architecture models, 100 realizations of each model setup are requested. The file format should be RSML and the file name should be of the form “modelname_replicate,” e.g., “CRootBox_1.rsml.”

2.2.1.1.3. Reference data analysis and automated model comparison. Statistical evaluation of a root architecture model has for example been done by Delory et al. (2018); Schnepf et al. (2018a). This motivated the creation of two data analysis pipelines for the first benchmark (M1.1) that will be used to compare simulation outputs with reference experimental data (reference root systems) (Figure 2A). These two data analysis pipelines are implemented in the Jupyter Notebook RSA calibration.ipynb that can be found on the github repository that contains code that will automatically include every model output in the analysis that is available in the prescribed folder. The analysis relies on the functions available in the R package archiDART (Delory et al., 2016, 2018). In the first pipeline, traits computed at the root system level (e.g., total root system length, number of roots per branching order) are compared between all simulated and reference root systems. This comparison takes place in three steps: (1) identifying the key morphological, architectural, and topological (Fitter indices, Fitter, 1987; Fitter and Stickland, 1991) traits explaining differences between simulated and reference root systems using multivariate data analysis techniques (e.g., discriminant analysis and principal component analysis), (2) looking at the point in time, beyond the time period for which there are measurements, when simulated and reference root systems start to diverge/converge with regard to the key root system traits identified in the previous step and how large these differences are, and (3) assess the degree of dissimilarity between simulated and reference root systems using dissimilarity metrics based on the raw data (Janssen and Heuberger, 1995).

FIGURE 2

Figure 2. Presentation of the data analysis pipelines used for the benchmarking of root architecture models. (A,B) Show the first (M1.1) and second (M1.2) benchmark scenarios, respectively.

In the second pipeline, dissimilarities in architecture between reference and simulated root systems are compared using persistent homology. Persistent homology is a topological framework that has proven to be a very powerful tool for capturing variations in plant morphology at different spatial scales (Li et al., 2017, 2018). The main output of a persistent homology analysis is a persistence barcode recording the appearance and disappearance of each root branch when a distance function traverses the branching structure (see Figure 1 in Delory et al., 2018). The degree of similarity between different root system topologies can be assessed by computing a pairwise distance matrix to compare persistence barcodes. In addition, Delory et al. (2018) showed that both trait-based and persistent homology approaches nicely complement each other and allow root researchers to more accurately describe differences in root system architecture (Delory et al., 2018). In our data analysis pipeline, a persistent homology analysis comprises the following steps: (1) computing a persistence barcode for each simulated and reference root system using a geodesic distance function, (2) computing dissimilarities between persistence barcodes using a bottleneck distance, (3) visualize dissimilarities between root systems using multidimensional scaling, and (4) test specific hypotheses using permutational multivariate analysis of variance (PERMANOVA) (Anderson, 2001).

2.2.1.2. M1.2 Long model simulations

In this benchmark, modelers use the same input parameter set as in M1.1, but simulate root system growth and development for a longer time period (60 days). The aim of this second benchmarking exercise is to assess if the different models diverge (or converge) if simulations are run for a longer time period and extrapolate beyond the provided data set (Figure 2B). This is of great importance, as parameterization of RSA models is often based on relatively young plants, whereas knowledge of RSA of older root systems is scarce. Therefore, for this M1.2 scenario, experimental data are not used as the basis of comparison anymore. It has to be noted that these two benchmark problems focus on root architecture dynamics modeling only, thus effect of soil properties on root growth is not explicitly modeled.

2.2.1.2.1. Required output. The following results are to be uploaded via pull requests to this path on the github repository: M1 Root architecture development/M1.2 RSA simulation/M1.2 Numerical results.

1. A text file including the model input parameters used for the specific simulator.

2. Simulation output from running the root architecture model using this parameter set in RSML format. Due to the stochastic nature of root architecture models, 100 realizations of each model setup are requested. The file format should be RSML and the file name should be of the form “modelname_replicate,” e.g., “CRootBox_1.rsml.”

2.2.1.2.2. Analysis pipeline for M1.2. For the second benchmark (M1.2), three data analysis pipelines are used to compare simulation outputs given by different root architecture models. For this benchmark, the reference experimental data cannot be used as a reference as data of 60 days old plants is not available. The first two data analysis pipelines for M1.2 are very similar to the ones described earlier for the M1.1 benchmark. First, model outputs are compared using morphological, architectural, and topological traits computed at the root system level. Second, differences in root system morphology are analyzed using persistent homology. In addition to these two analysis pipelines, we included a third one to analyse differences in vertical root distribution between root systems simulated with different root architecture models. To do so, we use the modeling approach described in Oram et al. (2018). Briefly, relative cumulative root length density [Y(d)] is computed using Equation (1).

$\begin{array}{ll}Y(d)={\displaystyle \sum _{i=0}^{i=d}}RLD(i)/{\displaystyle \sum _{i=0}^{\infty }}RLD(i)& (1)\end{array}$

Equation (2) is fitted to the computed Y(d) using a non-linear least square means fitting procedure. The fitting constant β is used to compare modeled rooting depth, with high β corresponding to deep rooting.

$\begin{array}{ll}Y(d)=1-{\beta }^d,& (2)\end{array}$

2.2.2. Module 2: Water Flow in Soil Only

In this module, we describe benchmark problems that only relate to water flow in soil. Water flow in soil is most commonly described by the Richards equation in three dimensions:

$\begin{array}{ll}\frac{\partial \theta }{\partial t}=\nabla ·(K(\theta )(\nabla {\psi }_s+e_3)),& (3)\end{array}$

where θ is the volumetric soil water content (cm³cm⁻³), K is the hydraulic conductivity (cm day⁻¹), ψ_s is the soil water pressure head (cm), and e₃ = (0, 0, 1) is the standard unit vector.

The relationship between soil water pressure head and water content is generally described by the water retention curve. In the following we will use the van Genuchten equation (Van Genuchten, 1980) to describe this curve specifying the soil moisture characteristic of specific soils. All participating simulators will solve the exact same equation (i.e., Equation 3), with the same initial and boundary conditions. Therefore, differences between the outputs of different simulators are numerical solution-level differences, i.e., due to numerical scheme and implementation. Different numerical solutions of the Richards equation have been analyzed before, and for some settings analytic solutions exist. We will use the benchmarks presented by Vanderborght et al. (2005) to benchmark the part of the participating functional structural root architecture models where water movement in soil is described. The analytical solutions provided in that paper are related to vertical changes in the soil profile only. As most functional-structural root architecture models have a 3D soil module, they will prescribe no-flux boundary conditions at the sides of a domain with 25 cm length and width for the numerical implementation of those problems.

In the following we will describe the benchmarks for water movement in soil. Table 3 gives an overview of the soil hydraulic properties that will be used throughout all the benchmarks involving water flow in soil.

TABLE 3

Table 3. Soil hydraulic taken from Vanderborght et al. (2005).

2.2.2.1. M2.1: Infiltration

This benchmark scenario is taken from Vanderborght et al. (2005). All parameters, initial and boundary conditions are given in Table 2 and are described below. For each of the soil types, sand, loam and clay, we consider the rate of infiltration into a soil with an initial homogeneous soil water pressure head of ψ_s = −400 cm. All profiles are 200 cm deep, at the top boundary we prescribe a constant influx of 100 cm d⁻¹ as long as the soil is still unsaturated, and a Dirichlet boundary condition of ψ_s = 0 cm as soon as the soil is fully saturated. Note that the prescribed infiltration value is high, such that in most scenarios, the boundary condition will switch to Dirichlet very soon. At the bottom boundary, we prescribe free drainage. Since this problem only produces gradients in the vertical direction, we compare numerical model results with the 1D analytical solution described in Vanderborght et al. (2005).

2.2.2.1.1. Reference solution. The analytical solution is given by the traveling wave equation

$\begin{array}{l}\Delta \eta (\theta )=\eta (\theta )-\eta ({\theta }_a)=({\theta }_{sur}-{\theta }_i)\\ \begin{array}{c}{\displaystyle {\int }_{\theta }^{{\theta }_a}}\frac{D_w(\theta )d\theta }{\left[K({\theta }_{sur})-K({\theta }_i)\right](\theta -{\theta }_i)-\left[K(\theta )-K({\theta }_i)\right]({\theta }_{sur}-{\theta }_i)},\end{array}& (4)\end{array}$

where D_w is the water diffusivity (defined as $D_w=K(\theta )\frac{\partial {\psi }_s}{\partial \theta }$), θ_sur is the water content at the soil surface, θ_i is the initial water content, θ_a is a reference water content (taken to be θ_a = (θ_sur + θ_i)/2), $\eta =\left|z\right|-\frac{\left[K({\theta }_{sur})-K({\theta }_i)\right]t}{{\theta }_{sur}-{\theta }_i}$ and Δη(θ) is the distance of the front to the position of the reference water content. The implementation of this analytical solution, implemented in the Jupyter Notebook M2.1 Benchmark problem.ipynb, reproduces Figures 4a–c from Vanderborght et al. (2005), where the water content is plotted after 0.1, 0.2, and 0.3 days for the sand scenario; 0.2, 0.5, and 1 days for the loam scenario; and 0.1, 0.2, and 0.5 days for the clay scenario (see Figure 3).

FIGURE 3

Figure 3. Results of M2.1: Infiltration into three initially dry soils: sand, loam, and clay.

2.2.2.1.2. Required output. The following simulation results of participating models are to be uploaded via pull requests to this path on the github repository: M2 Water flow in soil/M2.1 Infiltration/M2.1 Numerical results.

1. A text file consisting of nine pairs of rows containing comma separated depth values (cm) in the first, and water content (cm³cm⁻³) in the second row of the pair. The first three pairs represent the three time points for the sand scenario, the second three pairs represent the three time points of the loam scenario, and the last three pairs represent the three time points for the clay scenario. The file name should be of the form “simulatorname.txt,” e.g., “DuMux.txt.”

Note that we do not prescribe spatial or temporal resolution of the outputs, as that may depend on the individual numerical schemes.

2.2.2.2. M2.2: Evaporation

This benchmark reproduces Figure 5 of Vanderborght et al. (2005). We consider four scenarios (sand, loam 1, loam 2, clay) in which we are interested in the actual evaporation over time from an initially moist soil (ψ_i = −40 cm for the sand scenario and ψ_i = −200 cm for all other scenarios). The domain is 100 cm deep with a width and length of 10 cm. At the top boundary, we prescribe a constant efflux of J_s,pot = 0.1 cm d⁻¹ for the sand and loam 1 scenario, and 0.3 cm/day for the loam 2 and clay scenarios, at the bottom we prescribe zero-flux. When the soil reaches a critical soil water pressure head of −10,000 cm at the surface, we switch to a Dirichlet boundary condition with ψ_s = −10,000 cm.

2.2.2.2.1. Reference solution. The analytical solution to this problem is given by

$\begin{array}{ll}{J}_s(z=0,t)=\{\begin{array}{l}{J}_{s,pot}\text{ for }t<t_{pot}\hfill \\ \frac{S_w({\theta }_{sur},{\theta }_i)}{2\sqrt{t^{\prime }+t-t_{pot}}}\text{ for }t\ge t_{pot}\hfill \end{array}& (5)\end{array}$

where $t^{\prime }=\frac{S_w^2({\theta }_{sur},{\theta }_i)}{4J_{wpot}^2}$, $t_{pot}=\frac{S_w^2({\theta }_{sur},{\theta }_i)}{2J_{wpot}^2}$, $S_w({\theta }_i,{\theta }_{sur})=({\theta }_i-{\theta }_{sur})\sqrt{\frac{4}{\mu }{\int }_0^1{D}_w(\Theta )d\Theta }$, $\Theta =\left|\frac{\theta -{\theta }_{sur}}{{\theta }_i-{\theta }_{sur}}\right|$, $\mu =\frac{3\beta \left(1+{\left\{1-\frac{14}{9}\left[1-\frac{\alpha }{{(1-\beta )}^2}\right]\right\}}^{0.5}\right)}{2(1-\beta )\left[\frac{\alpha }{{(1-\beta )}^2}-1\right]}$, $\alpha =\frac{{\int }_0^1{(1-\beta \Theta )}^2{D}_w(\Theta )d\Theta }{{\int }_0^1{D}_w(\Theta )d\Theta }$, and $\beta ={\left[\frac{{\int }_0^1\Theta D_w(\Theta )d\Theta }{{\int }_0^1{D}_w(\Theta )d\Theta }\right]}^2$. Figure 4 shows the rate of evaporation over time for the four scenarios soil, loam 1, loam 2, clay.

FIGURE 4

Figure 4. Results of M2.2: Rate of evaporation with respect to time from sand with J_s,pot = 0 1 cm/d, loam with J_s,pot = 0 1 cm/d, loam with J_s,pot = 0 3 cm/d, and clay with J_s,pot = 0 3 cm/d.

2.2.2.2.2. Required output. The following simulation results of participating models are to be uploaded via pull requests to this path on the github repository: M2 Water flow in soil/M2.2 Evaporation/M2.2 Numerical results.

1. A text file consisting of two rows containing comma separated depth values (cm) in the first, and root pressure head (cm) in the second for each scenario [i.e., 4 (scenarios) × 2 (rows) = 8 rows]. The file name should be of the form “simulatorname.txt,” e.g., “DuMux.txt.”

Note that we do not prescribe spatial or temporal resolution of the outputs, as that may depend on the individual numerical schemes. It is the responsibility of each participant, to upload the best possible solution.

2.2.3. Module 3: Water Flow in Roots

In this benchmark, we consider water flow in xylem with constant and homogeneous soil water pressure head. This problem is well-described, e.g., in Doussan et al. (1998) and Roose and Fowler (2004). Its analytical solution for a single root was already derived by Landsberg and Fowkes (1978). In Appendix A, we present a derivation that is equivalent to the solution of Landsberg and Fowkes (1978) but uses exponential instead of hyperbolic functions. Briefly, conservation of mass in a branched root network with both axial and radial water flow, neglecting plant water storage and osmotic potential, yields Equation (6),

$\begin{array}{ll}2r_{root}\pi k_r({\psi }_s-{\psi }_x)=-k_x\frac{{\partial }^2{\psi }_x}{\partial {\zeta }^2},& (6)\end{array}$

where r_root is the root radius (cm), k_r is the radial conductivity (d⁻¹), ψ_s is the soil water pressure head of the surrounding soil (cm), ψ_x is the root water pressure head inside the xylem (cm), k_x is the axial conductance (cm³ d⁻¹), and ζ is the axial coordinate (cm).

2.2.3.1. M3.1: A single root in static soil with constant root hydraulic properties

In this benchmark problem, we assume a vertical single straight root segment surrounded by a soil with a constant and uniform soil water pressure head (i.e., the soil is not in hydrostatic equilibrium). We prescribe the root water pressure head at the root collar as ψ_x|_collar = ψ₀, and no axial flow at the root tips.

2.2.3.1.1. Reference solution. For constant k_r and k_x we can solve Equation (6) yielding

$\begin{array}{ll}{\psi }_x(\zeta )={\psi }_s+d_1{e}^{\sqrt{c}\zeta }+d_2{e}^{-\sqrt{c}\zeta },& (7)\end{array}$

with c = 2r_rootπk_r/k_x. The integration constants d₁ and d₂ for above boundary conditions are given by

$\begin{array}{ll}{d}_1=d^{-1}\left(e^{-\sqrt{c}{l}_{seg}}({\psi }_0-{\psi }_s)+1\right)& (8)\end{array}$

$\begin{array}{ll}{d}_2=-d^{-1}\left(e^{\sqrt{c}{l}_{seg}}({\psi }_0-{\psi }_s)+1\right),& (9)\end{array}$

where l_seg is the segment length, and d is the determinant of above matrix

$\begin{array}{ll}d=e^{-\sqrt{c}{l}_{seg}}-e^{\sqrt{c}{l}_{seg}},& (10)\end{array}$

see Appendix A. Figure 5 shows the analytical solution to this benchmark using the parameters given in Table 4.

FIGURE 5

Figure 5. Results of M3.1: Root water pressure head distribution within a single vertical root.

TABLE 4

Table 4. Parameters of scenario M3.1.

2.2.3.1.2. Required output. The following simulation results of participating models are to be uploaded via pull requests to this path on the github repository: M3 Water flow in roots/M3.1 Single root/M31 Numerical results/.

1. A text file consisting of two rows containing comma separated depth values (cm) in the first, and root pressure head (cm) in the second. The file name should be of the form “simulatorname.txt,” e.g., “DuMux.txt.”

Note that we do not prescribe spatial resolution of the outputs, as that may depend on the individual numerical schemes.

2.2.3.2. Benchmark M3.2: A small root system in a static soil

In the following benchmark, we extend benchmark M3.1 from a single root to a root system. We consider water flow inside a small static root system of a lupine plant which was grown for 14 days in a soil-filled column of 20 cm depth and 7 cm diameter. The root system was imaged by MRI at Forschungszentrum Jülich; the segmented root structure is provided in RSML, DGF (Dune grid format) (Bastian et al., 2008) and RSWMS (Javaux et al., 2008) formats in the folder M3 Water flow in roots/M3.2 Root system/root_grid on the github repository. It is visualized in Figures 6A,B with colors denoting root order and root segment age.

FIGURE 6

Figure 6. Visualization of the root system of M3.2 with colors denoting (A) root order, (B) root segment age, (C) root water pressure head.

2.2.3.2.1. Reference solution. The reference solution for this problem is given by the hybrid analytical-numerical solution of water flow in the root hydraulic architecture proposed by Meunier et al. (2017). The advantage of this solution is that it is independent of the spatial resolution of the root system (i.e., root segment length).

We consider two scenarios. The first one uses the same constant root hydraulic properties as given in Table 4, i.e., considering the same root hydraulic properties for each root segment. In the second scenario, we consider age-dependent root hydraulic properties for tap root and laterals of lupine as obtained by Zarebanadkouki et al. (2016, exponential function scenario) and converting distance from root tip to root age by assuming a root growth rate of 1 cm d⁻¹. This parameterization takes into account that roots get a higher axial conductivity and lower radial conductivity as they are becoming older (see Figure 7, a table with the actual values is provided on the github repository, in: M3 Water flow in roots/M3.2 Root system/M3.2 Benchmark problem.ipynb.

FIGURE 7

Figure 7. Root hydraulic properties dependency on root type and root segment age.

A sample 3-D visualization of the model output is shown in Figure 6C for the constant root hydraulic properties scenario. Figure 8 shows the effect of constant and age-dependent root hydraulic properties under otherwise same (soil and boundary) conditions.

FIGURE 8

Figure 8. Results of M3.2. (Left) Xylem pressure in each root segment of a root system with constant hydraulic properties. (Right) Xylem pressure in each root segment of a root system with age-dependent hydraulic properties.

2.2.3.2.2. Required output. The following simulation results of participating models are to be uploaded via pull requests to this path on the github repository: M3 Water flow in roots/M3.2 Root system/M32a Numerical results and M3 Water flow in roots/M3.2 Root system/M32b Numerical results for the constant and age-dependent root hydraulic properties cases.

1. A text file consisting of two rows containing comma separated depth values (cm) in the first, and root pressure head (cm) in the second. The file name should be of the form “simulatorname.txt,” e.g., “DuMux.txt.”

Note that we do not prescribe spatial resolution of the outputs, as that may depend on the individual numerical schemes.

2.2.4. Coupled Benchmark Scenarios C1: Root Water Uptake by a Static Root System

The way of coupling can easily introduce differences in simulated results because of numerical errors (especially when there is two way coupling) or because different assumption are made when implementing the coupling. No analytical solutions exists for the coupled problems presented here, but the coupling (C) benchmarks are intended to quantify differences between model outputs of coupled models. We may see differences observed in the non-coupled benchmarks to be amplified, or to be irrelevant for the coupled problem.

2.2.4.1. C1.1: Water uptake by a single root

This benchmark follows the paper of Schröder et al. (2008). Here we aim to see to what extent the different participating models can reproduce the hydraulic conductivity drop near the root surface under different soil conditions and transpiration demands. Thus, it requires the participating line-source based models to strongly increase the spatial resolution of the 3D soil domain. From this benchmark, we will learn, whether the spatial resolution required to reproduce radial soil water pressure head gradients would be in a feasible order of magnitude for larger soil-root systems or not. If not, there are approaches to estimate soil water pressure head drop at the root-soil interface from bulk soil values as e.g., in Schröder et al. (2009), Beudez et al. (2013), and Mai et al. (2019), see also benchmark C1.2.

2.2.4.2. Reference solution

The analytical solution is based on the analytical solutions of the 1D radially symmetric problem of water uptake by a single root, in which root water uptake is described as a boundary condition at the root-soil interface. We consider here two water uptake regimes, a non-stressed condition with maximum root uptake (q_root), and a stressed condition with a limiting plant root water potential constraining uptake. Based on the steady-rate assumption and using the matric flux potential $\Phi (h_c)={\int }_{-\infty }^{h_c}K(h)dh$ that linearizes the Richards equation, the radial soil water pressure head profiles for non-stressed and stressed conditions (stress conditions are given when the soil water pressure head at the root surface reaches −15,000 cm) are given by

$\begin{array}{l}{\Phi }_{nostress}(r)={\Phi }_{r_{out}}+(q_{root}{r}_{root}-q_{out}{r}_{out})\\ \left[\frac{r^2/r_{root}^2}{2(1-{\rho }^2)}+\frac{{\rho }^2}{1-{\rho }^2}\left(\text{ln}\frac{r_{out}}{r}-\frac{1}{2}\right)\right]\\ \begin{array}{c}+q_{out}{r}_{out}\text{ln}\frac{r}{r_{out}}\end{array}& (11)\end{array}$

and

$\begin{array}{l}{\Phi }_{stress}(r)=\left({\Phi }_{r_{out}}-{\Phi }_{r_{root}}+q_{out}{r}_{out}\text{ln}\frac{1}{\rho }\right)\\ \begin{array}{c}\frac{r^2/r_{root}^2-1+2{\rho }^2\text{ln}{r}_{root}/r}{{\rho }^2-1+2{\rho }^{2}ln1/\rho }+q_{out}{r}_{out}\text{ln}\frac{r}{r_{root}}+{\Phi }_{root},\end{array}& (12)\end{array}$

where $\rho =\frac{r_{out}}{r_{root}}$.

Given the soil water pressure head at the outer boundary, the solution computes the soil water pressure head profile toward the root. Due to the steady-rate assumption, the problem has become a stationary boundary value problem. However, under non-stressed conditions, we can calculate the time that corresponds to a given radial soil water pressure head profile by dividing the volume of water removed from the soil domain by the known water flow rate. The water remaining in a 1 cm long hollow cylinder around the root is given by

$\begin{array}{ll}V={\displaystyle {\int }_0^{2\pi }}{\displaystyle {\int }_{r_{root}}^{r_{out}}}r\theta drd\varphi =2\pi {\displaystyle {\int }_{r_{root}}^{r_{out}}}r\theta dr,& (13)\end{array}$

θ being the water content. The initially available water volume in the soil domain is given by

$\begin{array}{ll}{V}_i=(r_{out}^2-r_{root}^2)\pi {\theta }_i.& (14)\end{array}$

Thus, until the onset of stress, the corresponding time at which a given radial profile is reached is given by

$\begin{array}{ll}t=\frac{(V_i-V)}{2r_{root}\pi q_{root}}.& (15)\end{array}$

For the three soils sand, loam, and clay (Table 3), we compute the analytical solution for two different values of q_root (q_root = 0.1 cm/d and q_root = 0.05 cm/d, alternatively), and with the following parameters: r_root = 0 02 cm, r_out = 1 cm, ψ_s,lim = −15,000 cm, q_out = 0.0 cm/d, for different soil water pressure heads at the outer end of the cylinder. Figure 9 shows the soil water pressure head gradients at the onset of stress (i.e., when the soil water pressure head at the root surface reached −15,000 cm) and the time of its occurrence. The value of the initial soil water pressure head is taken to be ψ_s,i = −100 cm. This analytical solution is for radial water flow in soil toward the root only, i.e., not considering gravity or water flow inside the roots. Ideally, in their numerical implementation of this benchmark, the different participating models will turn off gravity effects. The soil domain for this numerical implementation has a size of l × w × d = 1 × 1 × 1 cm. The horizontal spatial resolution is high enough such that hydraulic conductivity drop near root surface can be resolved. The axial and radial conductances are high, such that the pressure inside the root is everywhere the same and the uptake flux is uniform.

FIGURE 9

Figure 9. Results of C1.1: Soil water pressure head gradients around a single, transpiring, root at the onset of stress and the time of its occurrence.

2.2.4.2.1. Required output. The following simulation results of participating models are to be uploaded via pull requests to this path on the github repository: M3 Water flow in roots/M3.2 Root system/M32a Numerical results and M3 Water flow in roots/M3.2 Root system/M32b Numerical results for the constant and age-dependent root hydraulic properties cases.

1. A text file consisting of two rows containing comma separated radial distances from the root surface (cm) in the first, and soil pressure head (cm) in the second for each soil and transpiration rate scenario [i.e., 3 (soils) × 2 (transpiration rates) × 2 = 12 rows]. The file name should be of the form “simulatorname.txt,” e.g., “DuMux.txt.”

Note that we do not prescribe spatial or temporal resolution of the outputs, as that may depend on the individual numerical schemes.

2.2.5. C1.2: Water Uptake by a Root System From Drying Soil

This benchmark scenario considers water uptake by a static 8-days-old lupine root system given in the public data set (Koch, 2019) as RSML or DGF. The root is the same as the one in benchmark M3.2, only younger, in order to reduce the computational cost for the reference scenario. The root system has been segmented from MRI measurements. The lupine is embedded in a soil box of l × w × d = 8 × 8 × 15 cm filled with loam (soil hydraulic properties given in Table 3). The benchmark is to evaluate the accuracy of root water uptake models under conditions of drying soil. To this end, the soil has an initial water content of θ_top = 0.129, corresponding to a pressure head ψ_{s, top} = −659.8 cm at the soil surface (z = 0). The pressure head in the rest of the domain initially follows a hydrostatic distribution

$\begin{array}{ll}{\psi }_{s,\text{i}}={\psi }_{s,\text{top}}-z,& (16)\end{array}$

where z (in cm) denotes the vertical position (upward-pointing axis, zero at soil surface). At all soil boundaries, as well as at the root tips, no-flux boundaries are prescribed. A potential transpiration rate is given as the sinusoidal diurnal function

$\begin{array}{ll}{Q}_{\text{pot}}(t)=\overline{Q}\left[1+\text{sin}\left(2\pi t-\frac{\pi }{2}\right)\right],& (17)\end{array}$

where the mean transpiration rate is $\overline{Q}=$ 6.4³ cm d⁻¹, the time t is given in days, and Q_pot(t = 0) = 0, that is, the simulation starts at night. The potential transpiration rate Q_pot, Equation (17), is enforced at the root collar (Neumann boundary condition) as long as the root water pressure head at the root collar is above ψ_x,crit = −15, 290 cm (corresponding to −1.5 MPa). If this critical root water pressure head at the root collar is reached, the boundary condition is switched to a Dirichlet type boundary condition, enforcing a constant pressure head ψ_x,crit = −15,290 cm at the root collar. This informal description is intentional, as the actual implementation of such a boundary condition may vary from simulator to simulator. We consider two scenarios. In scenario C1.2a the root hydraulic properties are constant. The tap root and lateral root conductivities are $k_x=4.32\times 10^{-23}$ cm d⁻¹ and $k_r=1.73\times 10^{-4}$ d⁻¹ (Table 4). For scenario C1.2b the root hydraulic properties depend on the root type and root age and are shown in Figure 7.

Given the soil domain Ω and the network of root center-lines Λ, we solve the following coupled system of equations

$\begin{array}{ll}\frac{\partial \theta }{\partial t}-\nabla ·(K(\theta )(\nabla {\psi }_s+e_3))=q({\psi }_x,{\psi }_s)\text{in }\Omega ,& (18)\end{array}$

$\begin{array}{ll}-\frac{\partial }{\partial \zeta }\left(k_x\frac{\partial {\psi }_x}{\partial \zeta }+\frac{\partial z}{\partial \zeta }\right)=\widehat{q}({\psi }_x,{\psi }_s)\text{on }\Lambda ,& (19)\end{array}$

subject to the boundary conditions specified above, where ζ is a scalar parameterization (local axial coordinate) of the root segments. The specific radial flux $\widehat{q}$ in units (cm2 d-1) is given by the difference in the average soil water pressure head on the root surface and in the xylem multiplied by the root radial conductivity. The formulation of q in Equation (18) may be different between different participating models. A discussion on singularity issues when evaluating the soil water content at the root center line can be found in Koch et al. (2018b). In many cases, the soil discretization is much larger than the root diameter, and thus the drop in hydraulic conductivity near the root surface in dry soils may not be sufficiently resolved in the soil domain. Different approaches for the determination of the sink term for root water uptake are likely to differ most in dry soil. The reference solution to this benchmark is designed to evaluate possible differences between the models in that regard.

2.2.6. Reference solution

As no analytical solutions exist for this problem of coupled water flow in the soil-root system, we designed a reference solution with a numerical model that explicitly considers the physical presence of roots in the soil domain, i.e., the soil mesh is highly refined around all roots and water uptake is modeled via boundary conditions at all the root surfaces. Thus, this reference solution does not make any assumptions that are inherent in the definition of the sink terms for root water uptake in the line source-based models. An explicit 3D soil grid is also used in Daly et al. (2018). However here, the soil is additionally coupled to the xylem flow in the root. The root is still modeled as a network of one-dimensional segments (center-line representation). Each segment has a specific radius as specified in the RSML grid file to this benchmark. A three-dimensional representation of the root system is implicitly given by the union of all spheres along the root center-lines. Using this implicit representation a soil grid excluding the root system was generated using the C++ geometry library CGAL (The CGAL Project, 2019). In order to reduce the number of vertices in the mesh, the mesh is locally refined around the root-soil interface. The resulting mesh is available in the Gmsh format (Geuzaine and Remacle, 2009) in the data set. For the evaluation of the radial flux, which is a coupling condition on the soil faces σ representing the root-soil interface, we integrate over each face

$\begin{array}{ll}{F}_r={\displaystyle {\int }_{\sigma }}{r}_{\text{root}}{k}_r({\psi }_s-{\psi }_x)\text{d}A.& (20)\end{array}$

While the soil water pressure head is defined on the face, the corresponding root xylem water pressure head has to be found by a mapping. To this end the integration point is first mapped onto the root surface using its implicit representation. Then the point is mapped onto the corresponding root center-line (a line segment) by finding the closest point on the line segment. There, ψ_x is evaluated. The flux is added as a source term in the corresponding segment in the root. The model is implemented in the open-source porous media simulator DuMu^x (Flemisch et al., 2011; Koch et al., 2018a; Koch et al., 2019). The coupled system is solved with a fully coupled manner, using Newton's method, and monolithic linear solver (block-preconditioned stabilized bi-conjugate gradient solver) in each Newton iteration. The equations are discretized in time with an implicit Euler method, and in space with a locally mass conservative vertex-centered finite volume method (BOX method Helmig, 1997). The maximum time step size is Δt = 1,200 s. The actual time step size may be sometimes chosen smaller, depending on the convergence speed of the Newton method. Output files are produced in regular intervals every 1,200 s starting with the initial solution. The simulation time is 3 d.

Soil water content and root water pressure head in a three-dimensional plot is shown in Figure 10 for C1.2b. Figure 11A shows the potential and actual transpiration rates for both scenarios, with constant and age-dependent root hydraulic properties. The curves hardly differ since the water pressure head drop is dominated by the low conductivity of the dry soil. In Figure 11B, the differences between scenarios are more clearly visible in terms of the minimal and maximal root water pressure head with respect to time.

FIGURE 10

Figure 10. C1.2: Root water uptake by a static root system over time. (Left) Visualization of the volumetric soil water content on vertical and horizontal slices through the soil domain and along the root surfaces. (Right) Root water pressure head.

FIGURE 11

Figure 11. Results of C1.2 for two scenarios, constant and age-dependent root hydraulic properties. (A) Actual transpiration of reference solution. (B) Root water pressure head distributions inside the root system.

2.2.6.1. Required outputs

To compare the results between the participating models, the desired outputs are

• VTK files (3D) of soil water pressure head and water content on the first, second, and third day (t = 0.5, 1.5, 2.5 d). For output written every 1,200 s this means the output files with the number 36, 108, and 180

• VTK files (lines in 3D) of root water pressure head in the first, second, and third day (t = 0.5, 1.5, 2.5 d)

• CSV file with three data points per time step (each 1,200 s starting with t = 0): time and actual transpiration rate

• CSV file with three data point per time step: time and minimum and maximum root water pressure head.

File names of the VTK files should indicate the simulator name, the state variable, the domain, and the output time in days, e.g., “DuMux_soil_theta_1d.vtk.” File names of the CSV files should indicate the simulator name and output time it days, e.g., “Dumux_1.csv.”

2.3. Coupled Benchmark Scenarios C2: Root Water Uptake by a Dynamic Root System

In this benchmark, we wish to explore differences caused by the approach of root growth modeling. We assess how the differences in root architecture parameters resulting from M1.2 propagate (or not) in the computation of the root water uptake from soil. In this example, we do not consider the effect of soil properties on root growth, but only the differences that arise from the different root systems according to M1.2.

2.3.1. C2.1: Water Uptake by a Single Root

Before looking at the root system, we look at how the implementation of the growth itself affects computed root water uptake for a single root. This scenario is analogous to C1.1, but with a single root growing at an elongation rate of 2 cm/d from 1 to 10 cm length.

2.3.1.1. Required outputs

The required outputs for model intercomparison are

• VTK files of 3D soil water pressure head and water content in soil at a temporal resolution of 1 day up until 60 days (point data)

• VTK files of xylem water pressure head (point data)

• Text files with two lines: time and corresponding actual transpiration.

2.3.2. C2.2: Water Uptake by a Root System

This scenario is the same as C1.2b, but replacing the static root system with a growing root system. The root growth parameters are for each model the results of M1.2; simulations start from a seed and run until a 60 days old root system. The domain size is 25 × 25 × 100 cm, the potential transpiration Q_pot = 0.5 cm³ d⁻¹ is scaled proportional to the root volume divided by the maximal root volume at maturity.

2.3.2.1. Required outputs

• VTK files of 3D soil water pressure head and water content in soil at a temporal resolution of 1 day up until 60 days (point data)

• VTK files of xylem water pressure head (point data)

• Text files with two lines: time and corresponding actual transpiration.

File names of the VTK files should indicate the simulator name, the state variable, the domain, and the output time in days, e.g., “DuMux_soil_theta_1d.vtk.” File names of the CSV files should indicate the simulator name and output time it days, e.g., “Dumux_1.csv.”

2.4. Automated Comparison Within All Benchmark Problems

Each benchmark folder on the github repository contains a Jupyter Notebook named “Automated comparison.” It provides the analytical solution of the respective benchmark and in addition includes Python code that automatically loads all the outputs of participating models that are provided in the “Numerical results” folder of that benchmark. As soon as new outputs are provided, they are automatically included in the analysis. Currently, different model outputs are already available. We envision more participating models' outputs to be provided in this way. Future analysis will include graphical and quantitative approaches.

3. Discussion

The benchmark problems considered here cover the basic processes of root water uptake from soil by 3D root architectures and focus on the coupling of root and soil domains. Root water capacity may be important in the case of trees (Janott et al., 2011) or when plants are under extreme water stress (Fang et al., 2019). Cavitation may induce a reduction of the specific root axial conductance (Sperry et al., 2003; Janott et al., 2011; Ahmad et al., 2018). Soil conditions can strongly affect root development (Schnepf et al., 2018b; de Moraes et al., 2019). At a later stage, these processes may be included in the benchmarking initiative by adding suitable benchmarking problems, e.g., including data from field studies, such as that of (de Moraes et al., 2019).

Root water uptake and evapotranspiration are a major factor in larger scale models, such as crop or land surface models (Kimball et al., 2019). They usually consider only the vertical soil dimension, thus have a 1-dimensional soil module. Thus, the functional-structural root architecture models considered here are not directly applicable. However, several examples have shown how the information of the 3-dimensional root hydraulic architecture can be implicitly considered in those models to compute root water uptake from 1D soils (Janott et al., 2011; Couvreur et al., 2014). In analogy to the electric circuit model, Couvreur et al. (2012) introduced a reduction of the 3-dimensional root hydraulic architecture to modular macroscopic equations and parameters operational for land surface and crop models (Baram et al., 2016; Sulis et al., 2019). Such a multiscale approach offers to connect the dots between models and measurable hydraulic and geometrical properties from the cell to the plant scale (Couvreur et al., 2018; Passot et al., 2018; Meunier et al., 2019) thus integrating essential processes and functional-structural properties for large scale models.

Additional processes, such as root water capacity and cavitation or the reduction of considered soil dimensions are out of the scope of this first initiative. However, we hope that it will function as a seed to initialize additional individual studies that consider those processes. We welcome such contributions in the Research Topic “Benchmarking 3D-Models of Root Growth, Architecture and Functioning” of “Frontiers in Plant Science.”

Numerical results of the different simulators will be compared to reference solutions or data where possible. For the root architecture models, measured root systems are available for comparison. The analysis pipelines for the RSA model outputs are outlined for the M1 module. The results of the different RSA simulators will be analyzed using both univariate and multivariate methods on root system traits as well as persistent homology. The soil and root water flow modules M2 and M3 have analytical solutions against which simulator results are compared. For the coupled problem with static root system, we offer a reference solution based on an explicit 3D simulation in which the root volume in the soil domain is accounted for. Quantitative comparison between different simulator results and reference solutions will rely on both residual-based and association-based goodness of fit measures (Bellocchi et al., 2010). The only benchmark problem without reference solution is the coupled problem with dynamic root architecture development, C2. Thus, for this problem the outcomes of the different simulators will be compared with each other. The aim is to obtain information about how diverse the different simulators are, and to quantify how the differences that arise from the RSA model choices (M1), the numerical implementation of soil and root water flow (M2-3) as well as the domain coupling choices (C1) propagate into the root water uptake computations. Based on these quantitative results, model users will be able to decide which model is suitable for a given application.

4. Conclusions

Functional-structural root architecture models have been compared qualitatively (e.g., Dunbabin et al., 2013), but until now no quantitative benchmarking existed. In other communities, benchmarking has been done or is ongoing, e.g., AgMIP (Porter et al., 2014) for crop models, CMIP (Eyring et al., 2016) for climate models, subsurface reactive transport models (Steefel et al., 2015). With this paper, we propose a framework for collaborative benchmarking of functional-structural root architecture models that allows quantitative comparison of the outputs of different simulators with reference solutions and with each other. This framework is presented using Jupyter Notebooks. Behind every “module” benchmark, there is a working code that explains and implements the reference solution or analysis of reference data. For both, “module” and “coupled” benchmarks, Jupyter Notebooks facilitate the automated comparison of simulator simulation outputs that are stored in specified folders of a public github repository. In this way, new numerical simulators that may be developed in the future may still be added to the automated comparison. All the analysis that is done in the Jupyter Notebooks is freely available so that the comparisons and analysis of uploaded model outputs will be transparent and repeatable. Future efforts will aim at extending the benchmarks from water flow in root and soil systems to further processes, such as solute transport, rhizodeposition, etc. We expect that this benchmarking will result in a better understanding of the different models and contribute toward improved models, with which we can simulate various scenarios with greater confidence. It will set standards for future model developments, ensuring that bugs, numerical errors or conceptual misunderstandings do not affect the value of future work. This is a step toward developing those models into the much-needed aid in the design of agricultural management schemes and model-guided crop breeding. These models may also be useful in ecology, e.g., to study species complementarity.

Data Availability Statement

The datasets generated for this study can be found in the github repository https://github.com/RSA-benchmarks/collaborative-comparison.

Author Contributions

AS initiated this benchmark initiative. AS, CB, VC, BD, CD, AK, TK, MJ, DL, GL, TM, FM, LP, JP, EP, VS, JV, HV, and MW together designed the benchmark problems. AS, VC, BD, AK, TK, MJ, DL, GL, FM, and JP contributed to the implementation of benchmark problems in Python or R. ML provided the MRI measurements of root architecture. All authors have contributed to the writing of the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer DH declared a past co-authorship with one of the authors, BD to the handling Editor.

Acknowledgments

AS acknowledges funding by the German Research Foundation (grant number SCHN 1361/3-1). VC was supported by the Belgian Fonds National de la Recherche Scientifique (FNRS, grant FC84104). VS acknowledges funding by the German Research Foundation (grant number SCHM 997/33-1). This research was institutionally funded by the Helmholtz Association (POF III Program—Research Fields Key Technologies for the Bioeconomy and Terrestrial Environment). CB acknowledges funding by the Advanced Research Projects Agency-Energy (ARPA-E), U.S. Department of Energy, under Award Number DE-AR0000821. This manuscript has been released as a pre-print at bioRxiv (Schnepf et al., 2019).

References

Ahmad, H. B., Lens, F., Capdeville, G., Burlett, R., Lamarque, L. J., and Delzon, S. (2018). Intraspecific variation in embolism resistance and stem anatomy across four sunflower (Helianthus annuus L.) accessions. Physiol. Plant. 163, 59–72. doi: 10.1111/ppl.12654

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, M. J. (2001). A new method for non-parametric multivariate analysis of variance. Aust. J. Ecol. 26, 32–46. doi: 10.1111/j.1442-9993.2001.01070.pp.x

CrossRef Full Text | Google Scholar

Baram, S., Couvreur, V., Harter, T., Read, M., Brown, P. H., Kandelous, M., et al. (2016). Estimating nitrate leaching to groundwater from orchards: comparing crop nitrogen excess, deep vadose zone data-driven estimates, and hydrus modeling. Vadose Zone J. 15. doi: 10.2136/vzj2016.07.0061

CrossRef Full Text | Google Scholar

Bastian, P., Blatt, M., Dedner, A., Engwer, C., Klöfkorn, R., Ohlberger, M., et al. (2008). A generic grid interface for parallel and adaptive scientific computing. Part I: abstract framework. Computing 82, 103–119. doi: 10.1007/s00607-008-0003-x

CrossRef Full Text | Google Scholar

Bellocchi, G., Rivington, M., Donatelli, M., and Matthews, K. (2010). Validation of biophysical models: issues and methodologies. A review. Agron. Sustain. Dev. 30, 109–130. doi: 10.1051/agro/2009001

CrossRef Full Text | Google Scholar

Beudez, N., Doussan, C., Lefeuve-Mesgouez, G., and Mesgouez, A. (2013). Influence of three root spatial arrangement on soil water flow and uptake. Results from an explicit and an equivalent, upscaled, model. Proc. Environ. Sci. 19, 37–46. doi: 10.1016/j.proenv.2013.06.005

CrossRef Full Text | Google Scholar

Bidel, L. P. R., Pagès, L., Rivière, L. M., Pelloux, G., and Lorendeau, J. Y. (2000). MassFlowDyn I: a carbon transport and partitioning model for root system architecture. Ann. Bot. 85, 869–886. doi: 10.1006/anbo.2000.1149

CrossRef Full Text | Google Scholar

Couvreur, V., Faget, M., Lobet, G., Javaux, M., Chaumont, F., and Draye, X. (2018). Going with the flow: multiscale insights into the composite nature of water transport in roots. Plant Physiol. 178, 1689–1703. doi: 10.1104/pp.18.01006

PubMed Abstract | CrossRef Full Text | Google Scholar

Couvreur, V., Vanderborght, J., Beff, L., and Javaux, M. (2014). Horizontal soil water potential heterogeneity: simplifying approaches for crop water dynamics models. Hydrol. Earth Syst. Sci. 18, 1723–1743. doi: 10.5194/hess-18-1723-2014

CrossRef Full Text | Google Scholar

Couvreur, V., Vanderborght, J., and Javaux, M. (2012). A simple three-dimensional macroscopic root water uptake model based on the hydraulic architecture approach. Hydrol. Earth Syst. Sci. 16, 2957–2971. doi: 10.5194/hess-16-2957-2012

CrossRef Full Text | Google Scholar

Daly, K. R., Tracy, S. R., Crout, N. M. J., Mairhofer, S., Pridmore, T. P., Mooney, S. J., et al. (2018). Quantification of root water uptake in soil using X-ray computed tomography and image-based modelling. Plant Cell Environ. 41, 121–133. doi: 10.1111/pce.12983

PubMed Abstract | CrossRef Full Text | Google Scholar

de Moraes, M. T., Debiasi, H., Franchini, J. C., Bonetti, J. A., Levien, R., et al. (2019). Mechanical and hydric stress effects on maize root system development at different soil compaction levels. Front. Plant Sci. 10:1358. doi: 10.3389/fpls.2019.01358

PubMed Abstract | CrossRef Full Text | Google Scholar

Delory, B. M., Baudson, C., Brostaux, Y., Lobet, G., du Jardin, P., Pagès, L., et al. (2016). archiDART: an R package for the automated computation of plant root architectural traits. Plant Soil 398, 351–365. doi: 10.1007/s11104-015-2673-4

CrossRef Full Text | Google Scholar

Delory, B. M., Li, M., Topp, C. N., and Lobet, G. (2018). archiDART v3.0: a new data analysis pipeline allowing the topological analysis of plant root systems. F1000Res. 7:22. doi: 10.12688/f1000research.13541.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Doussan, C., Pagès, L., and Vercambre, G. (1998). Modelling of the hydraulic architecture of root systems: an integrated approach to water absorption—model description. Ann. Bot. 81, 213–223. doi: 10.1006/anbo.1997.0540

CrossRef Full Text | Google Scholar

Dunbabin, V. (2007). Simulating the role of rooting traits in crop-weed competition. Field Crops Res. 104, 44–51. doi: 10.1016/j.fcr.2007.03.014

CrossRef Full Text | Google Scholar

Dunbabin, V. M., McDermott, S., and Bengough, A. G. (2006). Upscaling from rhizosphere to whole root system: modelling the effects of phospholipid surfactants on water and nutrient uptake. Plant Soil 283, 57–72. doi: 10.1007/s11104-005-0866-y

CrossRef Full Text | Google Scholar

Dunbabin, V. M., Postma, J. A., Schnepf, A., Pagès, L., Javaux, M., Wu, L., et al. (2013). Modelling root-soil interactions using three-dimensional models of root growth, architecture and function. Plant Soil 372, 93–124. doi: 10.1007/s11104-013-1769-y

CrossRef Full Text | Google Scholar

Dupuy, L. X., Fourcaud, T., Lac, P., and Stokes, A. (2007). A generic 3D finite element model of tree anchorage integrating soil mechanics and real root system architecture. Am. J. Bot. 94, 1506–1514. doi: 10.3732/ajb.94.9.1506

PubMed Abstract | CrossRef Full Text | Google Scholar

Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B., Stouffer, R. J., et al. (2016). Overview of the coupled model intercomparison project phase 6 (CMIP6) experimental design and organization. Geosci. Model Dev. 9, 1937–1958. doi: 10.5194/gmd-9-1937-2016

CrossRef Full Text | Google Scholar

Fang, Y., Yabusaki, S. B., Ahkami, A. H., Chen, X., and Scheibe, T. D. (2019). An efficient three-dimensional rhizosphere modeling capability to study the effect of root system architecture on soil water and reactive transport. Plant Soil 441, 33–48. doi: 10.1007/s11104-019-04068-z

CrossRef Full Text | Google Scholar

Fitter, A. H. (1987). An architectural approach to the comparative ecology of plant root systems. New Phytol. 106, 61–77. doi: 10.1111/j.1469-8137.1987.tb04683.x

CrossRef Full Text | Google Scholar

Fitter, A. H., and Stickland, T. R. (1991). Architectural analysis of plant root systems 2. Influence of nutrient supply on architecture in contrasting plant species. New Phytol. 118, 383–389. doi: 10.1111/j.1469-8137.1991.tb00019.x

CrossRef Full Text | Google Scholar

Flemisch, B., Darcis, M., Erbertseder, K., Faigle, B., Lauser, A., Mosthaf, K., et al. (2011). DuMu^x: DUNE for multi-$phase, component, scale, physics, …$ flow and transport in porous media. Adv. Water Resour. 34, 1102–1112. doi: 10.1016/j.advwatres.2011.03.007

CrossRef Full Text | Google Scholar

Geuzaine, C., and Remacle, J.-F. (2009). Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. Int. J. Numer. Methods Eng. 79, 1309–1331. doi: 10.1002/nme.2579

CrossRef Full Text | Google Scholar

Helmig, R. (1997). Multiphase Flow and Transport Processes in the Subsurface: A Contribution to the Modeling of Hydrosystems. Berlin; Heidelberg: Springer-Verlag.

Google Scholar

Hund, A., Trachsel, S., and Stamp, P. (2009). Growth of axile and lateral roots of maize: I development of a phenotying platform. Plant Soil 325, 335–349. doi: 10.1007/s11104-009-9984-2

CrossRef Full Text | Google Scholar

Janott, M., Gayler, S., Gessler, A., Javaux, M., Klier, C., and Priesack, E. (2011). A one-dimensional model of water flow in soil-plant systems based on plant architecture. Plant Soil 341, 233–256. doi: 10.1007/s11104-010-0639-0

CrossRef Full Text | Google Scholar

Janssen, P., and Heuberger, P. (1995). Calibration of process-oriented models. Ecol. Modell. 83, 55–66. doi: 10.1016/0304-3800(95)00084-9

CrossRef Full Text | Google Scholar

Javaux, M., Schröder, T., Vanderborght, J., and Vereecken, H. (2008). Use of a three-dimensional detailed modeling approach for predicting root water uptake. Vadose Zone J. 7, 1079–1088. doi: 10.2136/vzj2007.0115

CrossRef Full Text | Google Scholar

Kimball, B. A., Boote, K. J., Hatfield, J. L., Ahuja, L. R., Stockle, C., Archontoulis, S., et al. (2019). Simulation of maize evapotranspiration: an inter-comparison among 29 maize models. Agric. For. Meteorol. 271, 264–284. doi: 10.1016/j.agrformet.2019.02.037

CrossRef Full Text | Google Scholar

Koch, T. (2019). Benchmark C1.2–Numerical Results Reference Solution. doi: 10.18419/darus-471

CrossRef Full Text

Koch, T., Gläser, D., Weishaupt, K., Ackermann, S., Beck, M., Becker, B., et al. (2018a). Dumux 3.0.0. doi: 10.5281/zenodo.2479595

CrossRef Full Text

Koch, T., Gläser, D., Weishaupt, K., Ackermann, S., Beck, M., Becker, B., et al. (2019). DuMu^x 3-an open-source simulator for solving flow and transport problems in porous media with a focus on model coupling. arXiv 1909.05052.

Google Scholar

Koch, T., Heck, K., Schröder, N., Class, H., and Helmig, R. (2018b). A new simulation framework for soil-root interaction, evaporation, root growth, and solute transport. Vadose Zone J. 17. doi: 10.2136/vzj2017.12.0210

CrossRef Full Text | Google Scholar

Landl, M., Schnepf, A., Vanderborght, J., Bengough, A. G., Bauke, S. L., Lobet, G., et al. (2018). Measuring root system traits of wheat in 2d images to parameterize 3d root architecture models. Plant Soil 425, 457–477. doi: 10.1007/s11104-018-3595-8

CrossRef Full Text | Google Scholar

Landsberg, J., and Fowkes, N. (1978). Water movement through plant roots. Ann. Bot. 42, 493–508. doi: 10.1093/oxfordjournals.aob.a085488

CrossRef Full Text | Google Scholar

Li, M., Duncan, K., Topp, C. N., and Chitwood, D. H. (2017). Persistent homology and the branching topologies of plants. Am. J. Bot. 104, 1–5. doi: 10.3732/ajb.1700046

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, M., Frank, M. H., Coneva, V., Mio, W., Chitwood, D. H., and Topp, C. N. (2018). The persistent homology mathematical framework provides enhanced genotype-to-phenotype associations for plant morphology. Plant Physiol. 177, 1382–1395. doi: 10.1104/pp.18.00104

PubMed Abstract | CrossRef Full Text | Google Scholar

Lobet, G., Pagès, L., and Draye, X. (2011). A novel image-analysis toolbox enabling quantitative analysis of root system architecture. Plant Physiol. 157, 29–39. doi: 10.1104/pp.111.179895

PubMed Abstract | CrossRef Full Text | Google Scholar

Lobet, G., Pound, M. P., Diener, J., Pradal, C., Draye, X., Godin, C., et al. (2015). Root system markup language: toward a unified root architecture description language. Plant Physiol. 167, 617–627. doi: 10.1104/pp.114.253625

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, Y. Q., Randerson, J. T., Abramowitz, G., Bacour, C., Blyth, E., Carvalhais, N., et al. (2012). A framework for benchmarking land models. Biogeosciences 9, 3857–3874. doi: 10.5194/bg-9-3857-2012

CrossRef Full Text | Google Scholar

Mai, T. H., Schnepf, A., Vereecken, H., and Vanderborght, J. (2019). Continuum multiscale model of root water and nutrient uptake from soil with explicit consideration of the 3D root architecture and the rhizosphere gradients. Plant Soil 439, 273–292. doi: 10.1007/s11104-018-3890-4

CrossRef Full Text | Google Scholar

Meunier, F., Draye, X., Vanderborght, J., Javaux, M., and Couvreur, V. (2017). A hybrid analytical-numerical method for solving water flow equations in root hydraulic architectures. Appl. Math. Modell. 52, 648–663. doi: 10.1016/j.apm.2017.08.011

CrossRef Full Text | Google Scholar

Meunier, F., Heymans, A., Draye, X., Couvreur, V., Javaux, M., and Lobet, G. (2019). MARSHAL, a novel tool for virtual phenotyping of maize root system hydraulic architectures. in silico Plants. diz012. doi: 10.1093/insilicoplants/diz012

CrossRef Full Text | Google Scholar

Nygren, P., and Perttunen, J. (2010). “Rhizodeposition: a carbon efflux often neglected in functional-structural plant models,” in Proceedings of the 6th International Workshop on Functional-Structural Plant Models (Davis, CA).

Google Scholar

Oram, N. J., Ravenek, J. M., Barry, K., Weigelt, A., Chen, H., Gessler, A., et al. (2018). Below-ground complementarity effects in a grassland biodiversity experiment are related to deep-rooting species. J. Ecol. 106, 265–277. doi: 10.1111/1365-2745.12877

CrossRef Full Text | Google Scholar

Pagès, L., Bécel, C., Boukcim, H., Moreau, D., Nguyen, C., and Voisin, A.-S. (2014). Calibration and evaluation of archisimple, a simple model of root system architecture. Ecol. Modell. 290, 76–84. doi: 10.1016/j.ecolmodel.2013.11.014

CrossRef Full Text | Google Scholar

Pagès, L., Vercambre, G., Drouet, J., Lecompte, F., Collet, C., and Le Bot, J. (2004). Root typ: a generic model to depict and analyse the root system architecture. Plant Soil 258, 103–119. doi: 10.1023/B:PLSO.0000016540.47134.03

CrossRef Full Text | Google Scholar

Passot, S., Couvreur, V., Meunier, F., Draye, X., Javaux, M., Leitner, D., et al. (2018). Connecting the dots between computational tools to analyse soil–root water relations. J. Exp. Bot. 70, 2345–2357. doi: 10.1093/jxb/ery361

PubMed Abstract | CrossRef Full Text | Google Scholar

Porter, C. H., Villalobos, C., Holzworth, D., Nelson, R., White, J. W., Athanasiadis, I. N., et al. (2014). Harmonization and translation of crop modeling data to ensure interoperability. Environ. Modell. Softw. 62, 495–508. doi: 10.1016/j.envsoft.2014.09.004

CrossRef Full Text | Google Scholar

Postma, J. A., Kuppe, C., Owen, M. R., Mellor, N., Griffiths, M., Bennett, M. J., et al. (2017). Opensimroot: widening the scope and application of root architectural models. New Phytol. 215, 1274–1286. doi: 10.1111/nph.14641

PubMed Abstract | CrossRef Full Text | Google Scholar

Roose, T., and Fowler, A. (2004). A model for water uptake by plant roots. J. Theor. Biol. 228, 155–171. doi: 10.1016/j.jtbi.2003.12.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Schnepf, A., Black, C. K., Couvreur, V., Delory, B. M., Doussan, C., Koch, A., et al. (2019). Call for participation: collaborative benchmarking of functional-structural root architecture models. The case of root water uptake. bioRxiv. doi: 10.1101/808972

CrossRef Full Text | Google Scholar

Schnepf, A., Huber, K., Landl, M., Meunier, F., Petrich, L., and Schmidt, V. (2018a). Statistical characterization of the root system architecture model CRootBox. Vadose Zone J. 17. doi: 10.2136/vzj2017.12.0212

CrossRef Full Text | Google Scholar

Schnepf, A., Leitner, D., Landl, M., Lobet, G., Mai, T. H., Morandage, S., et al. (2018b). CRootBox: a structural–functional modelling framework for root systems. Ann. Bot. 121, 1033–1053. doi: 10.1093/aob/mcx221

PubMed Abstract | CrossRef Full Text | Google Scholar

Schröder, T., Javaux, M., Vanderborght, J., Körfgen, B., and Vereecken, H. (2008). Effect of local soil hydraulic conductivity drop using a three-dimensional root water uptake model. Vadose Zone J. 7. doi: 10.2136/vzj2007.0114

CrossRef Full Text | Google Scholar

Schröder, T., Javaux, M., Vanderborght, J., Körfgen, B., and Vereecken, H. (2009). Implementation of a microscopic soil–root hydraulic conductivity drop function in a three-dimensional soil–root architecture water transfer model. Vadose Zone J. 8, 783–792. doi: 10.2136/vzj2008.0116

CrossRef Full Text | Google Scholar

Schroeder, W., Martin, K., and Lorensen, B. (2006). The Visualization Toolkit. 4th Edn. Kitware Inc.

Sperry, J. S., Stiller, V., and Hacke, U. (2003). Xylem hydraulics and the soil–plant–atmosphere continuum. Agron. J. 95, 1362–1370. doi: 10.2134/agronj2003.1362

CrossRef Full Text | Google Scholar

Steefel, C. I., Yabusaki, S. B., and Mayer, K. U. (2015). Reactive transport benchmarks for subsurface environmental simulation. Comput. Geosci. 19, 439–443. doi: 10.1007/s10596-015-9499-2

CrossRef Full Text | Google Scholar

Sulis, M., Couvreur, V., Keune, J., Cai, G., Trebs, I., Junk, J., et al. (2019). Incorporating a root water uptake model based on the hydraulic architecture approach in terrestrial systems simulations. Agric. For. Meteorol. 260–270, 28–45. doi: 10.1016/j.agrformet.2019.01.034

CrossRef Full Text | Google Scholar

The CGAL Project (2019). CGAL User and Reference Manual. 4.14 Edn. CGAL Editorial Board. Available online at: https://www.cgal.org/index.html

Van Genuchten, M. (1980). Closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898. doi: 10.2136/sssaj1980.03615995004400050002x

CrossRef Full Text | Google Scholar

Vanderborght, J., Kasteel, R., Herbst, M., Javaux, M., Thiéry, D., Vanclooster, M., et al. (2005). A set of analytical benchmarks to test numerical models of flow and transport in soils. Vadose Zone J. 4, 206–221. doi: 10.2136/vzj2005.0206

CrossRef Full Text | Google Scholar

Zarebanadkouki, M., Meunier, F., Couvreur, V., Cesar, J., Javaux, M., and Carminati, A. (2016). Estimation of the hydraulic conductivities of lupine roots by inverse modelling of high-resolution measurements of root water uptake. Ann. Bot. 118, 853–864. doi: 10.1093/aob/mcw154

PubMed Abstract | CrossRef Full Text | Google Scholar

A. Appendix

A.1. Derivation of the Analytical Solution of Water Flow Inside the Root System

The axial water flow rate in the xylem Q_x (cm³ day⁻¹) is given by

$\begin{array}{ll}{Q}_x=-k_x\left(\frac{\partial {\psi }_x}{\partial \zeta }+v·e_3\right),& (\mathrm{A1})\end{array}$

where k_x is the axial conductance (cm³ day⁻¹), ψ_x is the pressure inside the xylem (cm), ζ is the local axial coordinate e₃ the unit vector in z-direction, and v the normalized direction of the xylem.

The radial water flow rate is given by

$\begin{array}{ll}{Q}_r=-2r_{root}\pi l_{seg}{k}_r({\psi }_s-{\psi }_x)& (\mathrm{A2})\end{array}$

with units (cm³ day⁻¹), where r_root is the root radius (cm), l_seg is the length of each root segment (cm), k_r is the radial conductivity (day⁻¹), and ψ_s is the soil water pressure head of the surrounding soil (cm). The equation is neglecting osmotic potential and is based on Equation (3.3) of Roose and Fowler (2004). Note that around the root a homogeneous soil water pressure head is assumed, therefore there is actually no hydrostatic equilibrium.

For each segment of length l_seg mass conservation yields

$\begin{array}{ll}0=Q_x{|}_{l_{seg}}-Q_x{|}_0+Q_r& (\mathrm{A3})\end{array}$

$\begin{array}{ll}-\frac{1}{l_{seg}}{Q}_r=-\frac{1}{l_{seg}}\left(Q_x{|}_0-Q_x{|}_{l_{seg}}\right)\text{ and for }{l}_{seg}\to 0& (\mathrm{A4})\end{array}$

$\begin{array}{ll}2r_{root}\pi k_r({\psi }_s-{\psi }_x)=-k_x\frac{{\partial }^2{\psi }_x}{\partial {\zeta }^2}& (\mathrm{A5})\end{array}$

see Equation (3.4) of Roose and Fowler (2004), where v₃ is the z-component of the normalized xylem direction (cm).

Integrating this ordinary differential equation leads to an explicit equation for ψ_x(ζ)

$\begin{array}{ll}{\psi }_x(\zeta )={\psi }_s+d_1{e}^{\sqrt{c}\zeta }+d_2{e}^{-\sqrt{c}\zeta },& (\mathrm{A6})\end{array}$

where c: = 2aπk_r/k_x, and d₁, and d₂ are integration constants that are derived from the boundary conditions.

To exemplify, we calculate d₁, and d₂ for a Dirichlet boundary condition at the root collar, and no-flux boundary conditions at the tip. The Dirichlet boundary conditions at the collar of the root system ψ_x|_collar = ψ₀ is inserted into the analytic solution (Equation A6), and yields

$\begin{array}{ll}{\psi }_s+d_1+d_2={\psi }_0.& (\mathrm{A7})\end{array}$

The Neumann boundary condition Q_x|_lseg = 0 (Equation A2) leads to

$\begin{array}{ll}\frac{\partial {\psi }_x}{\partial \zeta }{|}_{l_{seg}}=v_3,& (\mathrm{A8})\end{array}$

where l_seg is the length of the root segment. Using the derivation of the analytical solution yields

$\begin{array}{ll}{d}_1\sqrt{c}{e}^{\sqrt{c}{l}_{seg}}-d_2\sqrt{c}{e}^{-\sqrt{c}{l}_{seg}}=v_3.& (\mathrm{A9})\end{array}$

For a straight downward segment v₃ = −1, Equations (A7) and (A9) can be summarized as

$\begin{array}{ll}\begin{array}{ll}\left(\begin{array}{cc}1& 1\\ \sqrt{c}{e}^{\sqrt{c}{l}_{seg}}& -\sqrt{c}{e}^{-\sqrt{c}{l}_{seg}}\end{array}\right)\left(\begin{array}{c}{d}_1\\ d_2\end{array}\right)\hfill & =\left(\begin{array}{c}{\psi }_0-{\psi }_s\\ -1\end{array}\right)\hfill \end{array}& (\mathrm{A10})\end{array}$

Solving this linear equation for d₁ an d₂ yields

$\begin{array}{ll}{d}_1=d^{-1}\left(e^{-\sqrt{c}{l}_{seg}}({\psi }_0-{\psi }_s)+1\right)& (\mathrm{A11})\end{array}$

$\begin{array}{ll}{d}_2=-d^{-1}\left(e^{\sqrt{c}{l}_{seg}}({\psi }_0-{\psi }_s)+1\right),& (\mathrm{A12})\end{array}$

where d is the determinant of above matrix

$\begin{array}{ll}d=e^{-\sqrt{c}{l}_{seg}}-e^{\sqrt{c}{l}_{seg}}.& (\mathrm{A13})\end{array}$