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.

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

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."

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  x Xylem 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. 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.

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.

How to Participate
The participation includes three steps:  Analytical solution, Equations (11) and (12) Sand, loam, clay (  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." (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.

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) treelike 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.

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(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(Delory et al., , 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).
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. 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).

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.   Vanderborght et al. (2005). 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).
Equation (2) is fitted to the computed Y(d) using a nonlinear least square means fitting procedure. The fitting constant β is used to compare modeled rooting depth, with high β corresponding to deep rooting.

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: where θ is the volumetric soil water content (cm 3 cm −3 ), K is the hydraulic conductivity (cm day −1 ), ψ s is the soil water pressure head (cm), and e 3 = (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.

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 −1 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 where D w is the water diffusivity (defined as D w = K(θ ) ∂ψ s ∂θ ), θ 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), η = |z| − [K(θsur)−K(θi)]t θ sur −θ 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).

Required output.
The following simulation results of participating models are to be uploaded via pull requests

A text file consisting of nine pairs of rows containing comma
separated depth values (cm) in the first, and water content (cm 3 cm −3 ) 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.

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 −1 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.

Reference solution.
The analytical solution to this problem is given by . Figure 4 shows the rate of evaporation over time for the four scenarios soil, loam 1, loam 2, clay. 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.

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 welldescribed, 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), where r root is the root radius (cm), k r is the radial conductivity (d −1 ), ψ 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 3 d −1 ), and ζ is the axial coordinate (cm).

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 = ψ 0 , and no axial flow at the root tips.

Reference solution.
For constant k r and k x we can solve Equation (6) yielding with c = 2r root πk r /k x . The integration constants d 1 and d 2 for above boundary conditions are given by where l seg is the segment length, and d is the determinant of above matrix see Appendix A. Figure 5 shows the analytical solution to this benchmark using the parameters given in   (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.

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  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.

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 −1 . 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. 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.

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.

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.

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.

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 soilroot systems or not. If not, there are approaches to estimate soil water pressure head drop at the root-soil interface from bulk soil

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 steadyrate assumption and using the matric flux potential (h c ) = 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,000cm) are given by nostress (r) = r out + (q root r root − q out r out ) r 2 /r 2 root 2(1 − ρ 2 ) + ρ 2 1 − ρ 2 ln r out r − 1 2 + q out r out ln r r out and stress (r) = rout − rroot + q out r out ln 1 ρ (12) r 2 /r 2 root − 1 + 2ρ 2 lnr root /r ρ 2 − 1 + 2ρ 2 ln1/ρ + q out r out ln r r root + root , where ρ = 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 nonstressed 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 θ being the water content. The initially available water volume in the soil domain is given by Thus, until the onset of stress, the corresponding time at which a given radial profile is reached is given by 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 02cm, 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.
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.

C1.2: Water Uptake by a Root System From Drying Soil
This benchmark scenario considers water uptake by a static 8days-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).
subject to the boundary conditions specified above, where ζ is a scalar parameterization (local axial coordinate) of the root segments. The specific radial fluxq in units (cm 2 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.

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 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 vertexcentered 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 threedimensional 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.

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."

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.

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.

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.

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.5cm 3 d −1 is scaled proportional to the root volume divided by the maximal root volume at maturity.

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."

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.

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 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.

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/collaborativecomparison.