# Call for participation: Collaborative benchmarking of functional-structural root architecture models. The case of root water uptake

Andrea Schnepf, Christopher K Black, Valentin Couvreur, Benjamin M

Delory, Claude Doussan, Axelle Koch, Timo Koch, Mathieu Javaux,

Magdalena Landl, Daniel Leitner, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Andrea Schnepf, Christopher K Black, Valentin Couvreur, Benjamin M Delory, Claude Doussan, et al.,

Call for participation: Collaborative benchmarking of functional-structural root architecture models.

The case of root water uptake. 2019. hal-04210395

HAL Id: hal-04210395

https://hal.science/hal-04210395v1

Preprint submitted on 18 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Call for participation: Collaborative benchmarking of
functional-structural root architecture models.
The case of root water uptake.

Andrea Schnepf<sup>1,14,*</sup>, Christopher K. Black<sup>2</sup>, Valentin Couvreur<sup>3</sup>, Benjamin M. Delory<sup>4</sup>, Claude Doussan<sup>5</sup>, Axelle Koch<sup>6</sup>, Timo Koch<sup>7</sup>, Mathieu Javaux<sup>1,3</sup>, Magdalena Landl<sup>1,14</sup>, Daniel Leitner<sup>8</sup>, Guillaume Lobet<sup>1,14</sup>, Trung Hieu Mai<sup>1</sup>, Félicien Meunier<sup>9,10</sup>, Lukas Petrich<sup>11</sup>, Johannes A. Postma<sup>12</sup>, Eckart Priesack<sup>13</sup>, Volker Schmidt<sup>11</sup>, Jan Vanderborght<sup>1,14</sup>, Harry Vereecken<sup>1,14</sup>, and Matthias Weber<sup>11</sup>

<sup>1</sup> Institut für Bio- und Geowissenschaften: Agrosphäre (IBG-3),
Forschungszentrum Jülich GmbH, Wilhelm-Johnen-Str., D-52425 Jülich, Germany.

<sup>2</sup>Department of Plant Science, The Pennsylvania State University, 102 Tyson Building, University Park PA 16802, USA.

<sup>3</sup>Earth and Life Institute, Agronomy, Université catholique de Louvain, Louvain-la-Neuve, Belgium

<sup>4</sup>Institute of Ecology, Leuphana University Lüneburg, Universitätsallee 1, 21335 Lüneburg, Germany.

<sup>5</sup>UMR 1114 EMMAH, INRA/UAPV, 84914, Avignon cedex 9, France.

<sup>6</sup>Earth and Life Institute, Environmental Sciences, Université catholique de Louvain, Louvain-la-Neuve, Belgium

<sup>7</sup>Department of Hydromechanics and Modelling of Hydrosystems, University of Stuttgart, Pfaffenwaldring 61, 70569 Stuttgart, Germany.

<sup>8</sup>Simulationswerkstatt, Ortmayrstrasse 20, A-4060 Leonding, Austria.

<sup>9</sup>CAVELab - Computational and Applied Vegetation Ecology, Ghent University, Ghent, Belgium.

<sup>10</sup>Department of Earth and Environment, Boston University, Boston, USA.

<sup>11</sup>Institute of Stochastics, Ulm University, Helmholtzstr. 18, D-89069 Ulm, Germany.

<sup>12</sup>Institut für Bio- und Geowissenschaften: Plant Sciences (IBG-2),
Forschungszentrum Jülich GmbH, Wilhelm-Johnen-Str., D-52425 Jülich, Germany.

<sup>13</sup>Institute of Soil Ecology, Helmholtz Zentrum München, Ingolstädter Landstr. 1, 85764 Neuherberg, Germany.

<sup>14</sup>International Soil Modelling Consortium ISMC, Jülich, Germany

*Corresponding author, email: a.schnepf@fz-juelich.de

Abstract

Three-dimensional models of root growth, architecture and function are becoming im-
portant 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 archi-
tecture 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 com-
paring 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 towards 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 modelling 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 (Bidel et al., 2000, e.g.), rhizodeposition (Nygren and Perttunen, 2010), competition between species (Dun), 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 modelled 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. (2017) showed that the error made by neglecting root volumes physically present in the soil domain is negligibly small in case of root water uptake. Thus, models of this type are sufficiently accurate and computationally cheaper than explicit 3D. 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 modellers 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 modelling 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. (2017) 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. 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:

**(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 RSML (Lobet et al., 2015) for root architectures and VTK (Schroeder et al., 2006) for 3D and 1D simulation outputs. Python scripts
153 to read and write RSML files will be provided on the github repository. Packages to read and write
154 VTK files are for example available at https://pypi.org/project/vtk/.

155 (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 156 the code implemented in the Jupyter Notebooks for each benchmark problem. The final goal is to 157 jointly publish the results.
158
<table><caption>Table 1: List of notations</caption><thead><tr><th>Symbol</th><th>Units</th><th>Description</th></tr></thead><tbody><tr><td>d</td><td>cm</td><td>depth</td></tr><tr><td>D<sub>w</sub></td><td>cm<sup>2</sup>d<sup>-1</sup></td><td>water diffusivity</td></tr><tr><td>e<sub>3</sub></td><td>(0,0,1)</td><td>standard unit vector</td></tr><tr><td>J</td><td>cm<sup>3</sup>cm<sup>-2</sup>d<sup>-1</sup></td><td>water flux per unit soil surface area</td></tr><tr><td>k<sub>r</sub></td><td>cm<sup>3</sup>cm<sup>-2</sup>cm<sup>-1</sup>d<sup>-1</sup></td><td>root radial conductivity (defined as volume of water per unit root surface area, pressure head gradient and time)</td></tr><tr><td>k<sub>x</sub></td><td>cm<sup>4</sup>cm<sup>-1</sup>d<sup>-1</sup></td><td>specific root axial conductance</td></tr><tr><td>K(θ)</td><td>cm<sup>3</sup>cm<sup>-2</sup>d<sup>-1</sup></td><td>soil hydraulic conductivity</td></tr><tr><td>K<sub>sat</sub></td><td>cm<sup>3</sup>cm<sup>-2</sup>d<sup>-1</sup></td><td>saturated soil hydraulic conductivity</td></tr><tr><td>l</td><td>cm</td><td>length</td></tr><tr><td>n</td><td>-</td><td>van Genuchten shape parameter</td></tr><tr><td>q</td><td>cm<sup>3</sup>cm<sup>-2</sup>d<sup>-1</sup></td><td>water flux per unit root surface area</td></tr><tr><td>Q</td><td>cm<sup>3</sup>d<sup>-1</sup></td><td>volumetric water flow rate</td></tr><tr><td>Q&#773;</td><td>cm<sup>3</sup>d<sup>-1</sup></td><td>daily average volumetric water flow rate</td></tr><tr><td>Q<sub>r</sub></td><td>cm<sup>3</sup>d<sup>-1</sup></td><td>radial root water flow rate</td></tr><tr><td>Q<sub>x</sub></td><td>cm<sup>3</sup>d<sup>-1</sup></td><td>axial root water flow rate</td></tr><tr><td>r<sub>root</sub></td><td>cm</td><td>root radius</td></tr><tr><td>S<sub>w</sub></td><td>cm d<sup>-0.5</sup></td><td>sorptivity (infiltration) or desorptivity (evaporation)</td></tr><tr><td>t</td><td>d</td><td>time</td></tr><tr><td>v</td><td>(v<sub>1</sub>,v<sub>2</sub>, v<sub>3</sub>)</td><td>normalised direction of the xylem, pointing towards the root tip</td></tr><tr><td>w</td><td>cm</td><td>width</td></tr><tr><td>x, y, z</td><td>-</td><td>spatial coordinates, z-axis pointing upward, soil surface is at z=0</td></tr><tr><td>Y</td><td>-</td><td>cumulative root fraction from surface to depth d</td></tr><tr><td>α</td><td>cm<sup>-1</sup></td><td>van Genuchten shape parameter</td></tr><tr><td>β</td><td>-</td><td>root distribution index</td></tr><tr><td>η</td><td>cm</td><td>position of the infiltration front (eqn. (4))</td></tr><tr><td>λ</td><td>-</td><td>van Genuchten-Mualem parameter</td></tr><tr><td>Λ</td><td>-</td><td>root domain (network of root center-lines)</td></tr><tr><td>Ω</td><td>-</td><td>soil domain</td></tr><tr><td>Φ</td><td>cm<sup>2</sup>d<sup>-1</sup></td><td>matric flux potential</td></tr><tr><td>θ</td><td>cm<sup>3</sup>cm<sup>-3</sup></td><td>volumetric water content</td></tr><tr><td>θ<sub>a</sub></td><td>cm<sup>3</sup>cm<sup>-3</sup></td><td>reference water content</td></tr><tr><td>θ<sub>res</sub></td><td>cm<sup>3</sup>cm<sup>-3</sup></td><td>residual water content</td></tr><tr><td>θ<sub>sat</sub></td><td>cm<sup>3</sup>cm<sup>-3</sup></td><td>saturated water content</td></tr><tr><td>ψ</td><td>cm</td><td>water pressure head, described as potential energy per unit weight of water (i.e. units are cm of water column), given as relative to air pressure of 1020 cm and excluding the gravitational potential local coordinate along root axis</td></tr><tr><td>ζ</td><td>-</td><td>local coordinate along root axis</td></tr><tr><td colspan="3"><strong>sub indices</strong></td></tr><tr><td>collar</td><td>root collar (upper boundary of root system domain)</td><td></td></tr><tr><td>i</td><td>initial</td><td></td></tr><tr><td>pot</td><td>potential</td><td></td></tr><tr><td>r</td><td>radial</td><td></td></tr><tr><td>res</td><td>residual</td><td></td></tr><tr><td>s</td><td>soil</td><td></td></tr><tr><td>sat</td><td>saturation</td><td></td></tr><tr><td>seg</td><td>root segment</td><td></td></tr><tr><td>sim</td><td>simulation</td><td></td></tr><tr><td>sur</td><td>soil surface</td><td></td></tr><tr><td>tip(s)</td><td>root tip(s) (boundaries of root system domain)</td><td></td></tr><tr><td>top</td><td>top, position of the soil surface</td><td></td></tr><tr><td>out</td><td>outer radius of soil cylinder around a single root</td><td></td></tr><tr><td>x</td><td>xylem</td><td></td></tr></tbody></table>
<table><caption>Table 2: Description of benchmark scenarios to be implemented in 3D functional-structural root architecture models.<fnref n="1" /></caption><thead><tr><th></th><th>Benchmark problem</th><th>Domain</th><th>Initial conditions</th><th>Boundary conditions</th><th>Evaluation</th><th>Remarks</th></tr></thead><tbody><tr><td rowspan="2">RSA</td><td>M1.1: RSA calibration</td><td>t<sub>sim</sub>=11 (8) for lupine (maize)</td><td>seed position (0,0,-3)</td><td>n.a.</td><td>Comparison against the measured root systems provided - traits and persistent homology (PH)</td><td>Model parameters are determined from calibration against traced images provided in the github repository in RSML format in the folder in M1.1 RSA calibration/M1.1 Reference data; 100 realisations for each model setup</td></tr><tr><td>M1.2: RSA simulation</td><td>t<sub>sim</sub>=60</td><td>seed position (0,0,-3)</td><td>n.a.</td><td>No reference solution, comparison amongst models - traits, PH, RLD</td><td>RSA model parameters from M1.1; 10 realisations for each model setup</td></tr><tr><td rowspan="2">Soil</td><td>M2.1: Infiltration</td><td>1×w×d=10×10×200, t<sub>sim</sub>=1</td><td>ψ<sub>s,i</sub> = -400</td><td>at z = 0 { J<sub>s</sub> = 100 if ψ<sub>s</sub> &lt; 0 ,<br>ψ<sub>s</sub> = 0 else<br>∂ψ<sub>s</sub>/∂z|<sub>z=200</sub> = 1, no-flux at the sides</td><td>Analytical solution, eqn. (3)</td><td>sand, loam, clay (Table 3)</td></tr><tr><td>M2.2: Evaporation</td><td>1×w×d=10×10×100, t<sub>sim</sub>=10</td><td>ψ<sub>s,i</sub> = -40 for sand and -200 for all other scenarios</td><td>at z = 0 { J<sub>s</sub> = J<sub>s,ref</sub> if ψ<sub>s</sub> &gt; -10,000 ,<br>ψ<sub>s</sub> = -10,000 else<br>no-flux at all other boundaries</td><td>Analytical solution, eqn. (4)</td><td>scenario 1: sand, J<sub>s,ref</sub> = -0.1, scenario 2: loam, J<sub>s,ref</sub> = -0.1, scenario 3: loam, J<sub>s,ref</sub> = -0.3, scenario 4: clay, J<sub>s,ref</sub> = -0.3</td></tr><tr><td rowspan="2">Xylem</td><td>M3.1: Single root</td><td>1 vertical root, L=50</td><td>n.a.</td><td>ψ<sub>x|collar</sub> = -1000, Q<sub>r|tip</sub> = 0</td><td>Analytical solution, Eq. (11)</td><td>k<sub>x</sub>=0.0432, k<sub>r</sub>=1.73×10<sup>-4</sup>, ψ<sub>s</sub> = -200</td></tr><tr><td>M3.2: Root system</td><td>14-day old root system</td><td>n.a.</td><td>ψ<sub>x|collar</sub> = -500, Q<sub>r|tips</sub> = 0</td><td>Hybrid analytical solution (Meunier et al., 2017)</td><td>root hydraulic properties in scenario (a): Table 4, (b): Fig. 7, ψ<sub>s</sub> = -200, static RSA given in the root_grid folder of this benchmark</td></tr><tr><td rowspan="2">Coupled 1</td><td>C1.1: Single RWU</td><td>1D radially symmetric, r<sub>root</sub> = 0.02, r<sub>out</sub> = 0.6, t<sub>sim</sub>=20</td><td>ψ<sub>s,i</sub> = -100</td><td>at r = { r<sub>root</sub> { q<sub>r</sub> = q<sub>root</sub>, if ψ<sub>s</sub> &gt; -15,000 ,<br>ψ<sub>s</sub> = -15,000 else<br>q<sub>r</sub>|<sub>r=r<sub>out</sub></sub> = 0</td><td>Analytical solution, Eqs. (16) and (17)</td><td>sand, loam, clay (Table 3), scenarios 1-3: q<sub>root</sub>=0.1, scenarios 4-6: q<sub>root</sub>=0.05</td></tr><tr><td>C1.2: RWU, static RSA</td><td>static 8-day old root system, soil: 1×w×d = 8x8x15, t<sub>sim</sub>=3</td><td>ψ<sub>s,i</sub> = -659.8 - z</td><td>Q<sub>r|collar</sub> = { Q<sub>r</sub> = 6.4 if ψ<sub>s</sub> &gt; -15,290 ,<br>ψ<sub>s</sub> = -15,290 else<br>Q<sub>r|tips</sub> = 0, no-flux at all soil faces</td><td>Reference solution: explicit 3D simulation</td><td>loam (Table 3), static RSA given in the root_grid folder of this benchmark, root hydraulic properties in scenario (a): Table 4, (b): Fig. 7</td></tr><tr><td>Coupled 2</td><td>C2.1: RWU, dynamic RSA</td><td>growing root system, soil: 1×w×d = 25x25x100, t<sub>sim</sub>=60</td><td>ψ<sub>s,i</sub> = -200</td><td>Q<sub>r|collar</sub> = { Q<sub>r</sub> = 0.5 × rel<sub>LAI</sub> if ψ<sub>x</sub> &gt; -15,000 ,<br>ψ<sub>x</sub> = -15,000 else<br>Q<sub>r|tips</sub> = 0, no-flux at all soil faces</td><td>No reference solution, comparison amongst models</td><td>loam (Table 3), k<sub>x</sub>=0.0432, k<sub>r</sub>=1.73×10<sup>-4</sup>, RSA parameters from M1.1, rel<sub>LAI</sub> scales the potential transpiration</td></tr></tbody></table>

<footnote n="1">¹ All paths are relative to the github repository https://github.com/RSA-benchmarks/collaborative-comparison.git. For other abbreviations and units see Table 1.</footnote>
## 2.1 Benchmarks for individual modules

### 2.1.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 parameterised on the basis of a common experimental data set (Fig. 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) Parameterising 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 Fig. 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, i.e. we are looking at process-level differences between the different models. Thus, 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), modellers simulate root systems for the same duration as the age of the root systems in the reference dataset.

**Reference data set** Although the parameterisation 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 analysed using the semi-automated root image analysis software SmartRoot (Lobet et al., 2011) and root tracings were saved as RSML files for further analysis (Fig. 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 modellers to parameterise their respective RSA model.

**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 realisations 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".
<figure><img src="image_3.png" /><figcaption>Figure 1: Example of root images used for the benchmarking dataset. Panel (A) shows an image of lupin root systems, 11 days old, growing in an aeroponic setup. Panel (B) shows an image of a maize root system growing on filter paper (5 days old). All images were analysed using the semi-automated root image analysis software SmartRoot (Lobet et al., 2011), colours 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".</figcaption></figure>

**Reference data analysis and automated model comparison** Statistical evaluation of a root architecture model has for example been done by Schnepf et al. (2018); Delory et al. (2018). 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) (Fig. 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).

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 Fig. 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, modellers 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 (Fig. 2b). This is of great importance, as parameterisation 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 modelling only, thus effect of soil properties on root growth is not explicitly modelled.

**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 paramter set in RSML format. Due to the stochastic nature of root architecture models, 100 realisations of each model setup are requested. The file file format should be RSML and the file name should be of the form "modelname_replicate", e.g. "CRootBox_1.rsml".

**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 day 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 analysed 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 modelling approach described in Oram et al. (2018). Briefly, relative cumulative root length density (Y(d)) is computed using Eq. (1)

$$Y(d) = \sum_{i=0}^{i=d} RLD(i) / \sum_{i=0}^{\infty} RLD(i) \qquad (1)$$

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

$$Y(d) = 1 - \beta^d, \qquad (2)$$
<figure><img src="image_4.png" /><figcaption>Figure 2: Presentation of the data analysis pipelines used for the benchmarking of root architecture models. Panels a and b show the first (M1.1) and second (M1.2) benchmark scenarios, respectively.</figcaption></figure>

## 2.1.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:

$$ \frac{\partial \theta}{\partial t} = \nabla \cdot (K(\theta)(\nabla \psi_s + \mathbf{e}_3)), \quad (3) $$

where $\theta$ is the volumetric soil water content (cm$^3$cm$^{-3}$), K is the hydraulic conductivity (cm day$^{-1}$), $\psi_s$ is the soil water pressure head (cm), and $\mathbf{e}_3 = (0, 0, 1)$ is the downward 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.

We expect differences between the outputs of different simulators to be mainly numerical solution-level differences, i.e., due to numerical scheme and implementation. Different numerical solutions of the Richards equation have been analysed 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
289 in soil is described. The analytical solutions provided in that paper are related to vertical changes
290 in the soil profile only. As most functional-structural root architecture models have a 3D soil mod-
291 ule, they will prescribe no-flux boundary conditions at the sides of a domain with 25 cm length
292 and width for the numerical implementation of those problems.

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

<table><caption>Table 3: Soil hydraulic taken from Vanderborght et al. (2005). $\theta_{res}$ is the residual water content, $\theta_{sat}$ is the saturated water content, $\alpha$ and $n$ are the van Genuchten parameters, $K_{sat}$ is the saturated soil hydraulic conductivity and $\lambda$ is the van Genuchten-Mualem parameter</caption><thead><tr><th>Soil type</th><th>θ<sub>res</sub><br>(-)</th><th>θ<sub>sat</sub><br>(-)</th><th>α<br>(cm<sup>-1</sup>)</th><th>n<br>(-)</th><th>K<sub>s</sub><br>(cm d<sup>-1</sup>)</th><th>λ<br>(-)</th></tr></thead><tbody><tr><td>sand</td><td>0.045</td><td>0.43</td><td>0.15</td><td>3.0</td><td>1000</td><td>0.5</td></tr><tr><td>loam</td><td>0.08</td><td>0.43</td><td>0.04</td><td>1.6</td><td>50</td><td>0.5</td></tr><tr><td>clay</td><td>0.1</td><td>0.40</td><td>0.01</td><td>1.1</td><td>10</td><td>0.5</td></tr></tbody></table>

295

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 $\psi_s = -400$ cm. All profiles are 200 cm deep, at the top boundary we prescribe a constant influx of $100 \text{ cm d}^{-1}$ as long as the soil is still unsaturated, and a Dirichlet boundary condition of $\psi_s = 0$ cm as soon as the soil is fully saturated. 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).

305

**Reference solution** The analytical solution is given by the travelling wave equation

$$ \Delta\eta(\theta) = \eta(\theta) - \eta(\theta_a) = (\theta_{sur} - \theta_i) \int_{\theta}^{\theta_a} \frac{D_w(\theta)d\theta}{[K(\theta_{sur}) - K(\theta_i)](\theta - \theta_i) - [K(\theta) - K(\theta_i)](\theta_{sur} - \theta_i)}, \quad (4) $$

306 where $D_w$ is the water diffusivity (defined as $D_w = K(\theta) \frac{\partial \psi_s}{\partial \theta}$), $\theta_{sur}$ is the water content at the soil surface, $\theta_i$ is the initial water content, $\theta_a$ is a reference water content (taken to be $\theta_a = (\theta_{sur} + \theta_i)/2$), $\eta = |z| - \frac{[K(\theta_{sur}) - K(\theta_i)]t}{\theta_{sur} - \theta_i}$ and $\Delta\eta(\theta)$ is the distance of the front to the the position of the reference water content. The implementation of this analytical solution, implemented in the Jupyter Notebook M2.1 Benchmark problem.ipynb, reproduces Figure 4abc 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 Fig. (3)).

314

315

316

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

317

318

319

320

1. A text file consisting of two rows containing comma separated depth values (cm) in the first, and water content (cm³cm⁻³) in the second for each time point and infiltration scenario (i.e. 3 (time points) × 3 (scenarios) results = 18 rows). The file name should be of the form "simulatorname.txt", e.g. "DuMux.txt".

321

322

Note that we do not prescribe spatial or temporal resolution of the outputs, as that may depend on the individual numerical schemes.
<figure><img src="image_5.png" /><figcaption>Figure 3: Results of M2.1: Infiltration into three initially dry soils: sand, loam and clay.</figcaption></figure>

**M2.2: Evaporation** This benchmark reproduces Fig. (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 ($\psi_i = -40$cm for the sand scenario and $\psi_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 $\psi_s = -10.000$ cm.

**Reference solution** The analytical solution to this problem is given by

$$J_s(z = 0, t) = \begin{cases} J_{s,pot} & \text{for } t < t_{pot} \\ \frac{S_w(\theta_{sur}, \theta_i)}{2\sqrt{t' + t - t_{pot}}} & \text{for } t \ge t_{pot} \end{cases} \quad (5)$$

where $t' = \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(1 + \left\{1 - \frac{14}{9}\left[1 - \frac{\alpha}{(1-\beta)^2}\right]\right\}^{0.5})}{2(1-\beta)[\frac{\alpha}{(1-\beta)^2} - 1]}$, $\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$. Fig. (4) shows the rate of evaporation over time for the four scenarios soil, loam 1, loam 2, clay.

**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.
<figure><img src="image_6.png" /><figcaption>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</figcaption></figure>

### 2.1.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 Roosevelt and Fowler (2004) and Doussan et al. (1998). 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 Eq. (6),

$$2r_{root}\pi k_r(\psi_s - \psi_x) = -k_x \frac{\partial^2 \psi_x}{\partial \zeta^2}, \qquad (6)$$

where $r_{root}$ is the root radius (cm), $k_r$ is the radial conductivity (d$^{-1}$), $\psi_s$ is the soil water pressure head of the surrounding soil (cm), $\psi_x$ is the root water pressure head inside the xylem (cm), $k_x$ is the axial conductance (cm$^3$ d$^{-1}$), and $\zeta$ 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
358 prescribe the root water pressure head at the root collar as $\psi_x|_{collar} = \psi_0$, and no axial flow at the root tips.

360 **Reference solution** For constant $k_r$ and $k_x$ we can solve Eq. (6) yielding

$$ \psi_x(\zeta) = \psi_s + d_1 e^{\sqrt{c}\zeta} + d_2 e^{-\sqrt{c}\zeta}, \qquad (7) $$

with $c = 2r_{root}\pi k_r/k_x$. The integration constants $d_1$ and $d_2$ for above boundary conditions are given by

$$ d_1 = d^{-1} \left( e^{-\sqrt{c}l_{seg}}(\psi_0 - \psi_s) + 1 \right) \qquad (8) $$

$$ d_2 = -d^{-1} \left( e^{\sqrt{c}l_{seg}}(\psi_0 - \psi_s) + 1 \right), \qquad (9) $$

(10)

where $l_{seg}$ is the segment length, and $d$ is the determinant of above matrix

$$ d = e^{-\sqrt{c}l_{seg}} - e^{\sqrt{c}l_{seg}}, \qquad (11) $$

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

<table><caption>Table 4: Parameters of scenario M3.1.</caption><tr><td>l</td><td>50</td><td>length of a single straight root (cm)</td></tr><tr><td>r<sub>root</sub></td><td>0.02</td><td>radius (cm)</td></tr><tr><td>k<sub>z</sub></td><td>4.32 × 10<sup>-2</sup></td><td>axial conductivity (cm<sup>3</sup> d<sup>-1</sup>)</td></tr><tr><td>k<sub>r</sub></td><td>1.73 × 10<sup>-4</sup></td><td>radial conductivity (d<sup>-1</sup>)</td></tr><tr><td>ψ<sub>s</sub></td><td>-200</td><td>static soil water pressure head (cm)</td></tr><tr><td>ψ<sub>0</sub></td><td>-1,000</td><td>Dirichlet boundary conditions at the root collar (cm)</td></tr></table>

<figure><img src="image_7.png" /><figcaption>Figure 5: Results of M3.1: Root water pressure head distribution within a single vertical root</figcaption></figure>

363 **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/.

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

369 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 (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 visualised in Fig. 6(a,b) with colours denoting root order and root segment age.

<figure><img src="image_8.png" /><figcaption>Figure 6: Visualisation of the root system of M3.2 with colours denoting (a) root order, (b) root segment age, (c) root water pressure head.</figcaption></figure>

**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 parameterisation takes into account that roots get a higher axial conductivity and lower radial conductivity as they are becoming older (see Fig. 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 visualisation of the model output is shown in Fig. 6(c) for the constant root hydraulic properties scenario. Fig. 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.

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.
<figure><img src="image_9.png" /><figcaption>Figure 7: Root hydraulic properties dependency on root type and root segment age.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: Results of M3.2. Left panel: Xylem pressure in each root segment of a root system with constant hydraulic properties. Right panel: Xylem pressure in each root segment of a root system with age-dependent hydraulic properties.</figcaption></figure>

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

**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 Mai et al. (2019); Beudez et al. (2013), see also benchmark C1.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 linearises 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 -15000cm) are given by

$$ \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( \ln \frac{r_{out}}{r} - \frac{1}{2} \right) \right] + q_{out}r_{out} \ln \frac{r}{r_{out}} \quad (12) $$

and
$$ \Phi_{stress}(r) = \left( \Phi_{r_{out}} - \Phi_{r_{root}} + q_{out}r_{out} \ln \frac{1}{\rho} \right) \frac{r^2/r_{root}^2 - 1 + 2\rho^2 \ln r_{root}/r}{\rho^2 - 1 + 2\rho^2 \ln 1/\rho} + q_{out}r_{out} \ln \frac{r}{r_{root}} + \Phi_{root}, \quad (13) $$

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

Given the soil water pressure head at the outer boundary, the solution computes the soil water pressure head profile towards 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

$$ V = \int_{0}^{2\pi} \int_{r_{root}}^{r_{out}} r\theta dr d\phi = 2\pi \int_{r_{root}}^{r_{out}} r\theta dr $$

, $\theta$ being the water content. The initially available water volume in the soil domain is given by

$$ V_i = (r_{out}^2 - r_{root}^2)\pi\theta_i. $$

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

$$ t = \frac{(V_i - V)}{2r_{root}\pi q_{root}}. $$

For the three soils sand, loam, and clay (Table 3), we compute the analytical solution with the following parameters: $r_{root} = 0.02\text{cm}$, $r_{out} = 1\text{cm}$, $q_{root} = 0.5\text{cm/d}$, $\psi_{s,lim} = -15000\text{cm}$, $q_{out} = 0\text{cm/d}$ and for different soil water pressure heads at the outer end of the cylinder. Fig. 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 -15000cm) and the time of its occurrence. The value of the initial water content is taken to be $\theta_i = -100\text{cm}$. This analytical solution is for radial water flow in soil towards 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 \times w \times d = 1 \times 1 \times 1\text{ 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.

**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.
<figure><img src="image_11.png" /><figcaption>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</figcaption></figure>

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.1.5 C1.2: Water uptake by a root system from drying soil

This benchmark scenario considers water uptake by a static 8-day-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 \times w \times d = 8 \times 8 \times 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 $\theta_{\text{top}} = 0.129$, corresponding to a pressure head $\psi_{\text{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

$$ \psi_{s,i} = \psi_{s,\text{top}} - z, \qquad (14) $$

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

$$ Q_{\text{pot}}(t) = \bar{Q} \left[ 1 + \sin \left( 2\pi t - \frac{\pi}{2} \right) \right], \qquad (15) $$

where the mean transpiration rate is $\bar{Q} = 6.4 \text{ cm}^3 \text{ d}^{-1}$, the time $t$ is given in days, and $Q_{\text{pot}}(t = 0) = 0$, that is, the simulation starts at night. The potential transpiration rate $Q_{\text{pot}}$, Eq. (15), is enforced at the root collar (Neumann boundary condition) as long as the root water pressure head at the root collar is above $\psi_{x,crit} = -15.290 \text{ cm}$ (corresponding to $-1.5 \text{ 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 $\psi_{x,crit} = -15.290 \text{ 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^{-2} \text{ cm}^3 \text{ d}^{-1}$ and $k_r = 1.73 \times 10^{-4} \text{ d}^{-1}$ (Table 4). For scenario C1.2b the root hydraulic properties depend on the root type and root age and are shown in Fig. 7.

Given the soil domain $\Omega$ and the network of root center-lines $\Lambda$, we solve the following coupled system of equations

$$ \frac{\partial \theta}{\partial t} - \nabla \cdot (K(\theta)(\nabla \psi_s + e_3)) = q(\psi_x, \psi_s) \quad \text{in } \Omega, \qquad (16) $$

$$ -\frac{\partial}{\partial \zeta} \left( k_x \frac{\partial \psi_x}{\partial \zeta} + \frac{\partial z}{\partial \zeta} \right) = \hat{q}(\psi_x, \psi_s) \quad \text{on } \Lambda, \qquad (17) $$

subject to the boundary conditions specified above, where $\zeta$ is a scalar parametrisation (local axial coordinate) of the root segments. The specific radial flux $\hat{q}$ in units ($\text{cm}^2 \text{ d}^{-1}$) is given by the average soil water pressure head on the root surface. The formulation of $q$ in Eq. (16) 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 discretisation 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. Approaches to estimate soil water pressure head drop at the root-soil interface from bulk soil values can be found in Mai et al. (e.g. 2019); Beudez et al. (e.g. 2013). 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 modelled 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. (2017). 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 $\sigma$ representing the root-soil interface, we integrate over each face

$$ F_r = \int_{\sigma} r_{\text{root}} k_r (\psi_s - \psi_x) dA. \qquad (18) $$
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, $\psi_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., 2019; Koch et al., 2018a). 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 et al., 1997)). The maximum time step size is $\Delta t = 1200$ 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 1200 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 Fig. 10 for C1.2b. Fig. 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 Fig. 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><img src="image_12.png" /><figcaption>Figure 10: C1.2: Root water uptake by a static root system over time. Soil colours denote volumetric water content, root colours denote root water pressure head.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>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.</figcaption></figure>

**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$ d, $1.5$ d, $2.5$ d). For output written every 1200 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 day ($t = 0.5$ d, $1.5$ d, $2.5$ d)

* CSV file with three data points per time step (each 1200 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 in days, e.g., "Dumux_1.csv".
## 2.2 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 modelling.

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

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

**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.2.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 day old root system. The domain size is 25×25×100 cm, the potential

transpiration $Q_{pot} = 0.5 \text{ cm}^3\text{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".

## 2.3 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 Conclusions

Functional-structural root architecture models have been compared qualitatively (Dunbabin et al.,

2013, e.g.), but until now no quantitative benchmarking existed. In other communities, bench-

marking 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 towards 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 towards 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.
## 570 Acknowledgements

571 A.S. acknowledges funding by the German Research Foundation (grant number SCHN 1361/3-1).

572 V.C. was supported by the Belgian Fonds National de la Recherche Scientifique (FNRS, grant

573 FC84104). V.S. acknowledges funding by the German Research Foundation (grant number SCHM

574 997/33-1). This research was institutionally funded by the Helmholtz Association (POF III Pro-

575 gram—Research Fields Key Technologies for the Bioeconomy and Terrestrial Environment).

## 576 References

577 Simulating the role of rooting traits in crop-weed competition. *Field Crops Research*, 104:44–51.

578 Anderson, M. J. (2001). A new method for non-parametric multivariate analysis of variance.

579 *Australian Journal of Ecology*, 26:32–46.

580 Bastian, P., Blatt, M., Dedner, A., Engwer, C., Klöfkorn, R., Ohlberger, M., and Sander, O.

581 (2008). A generic grid interface for parallel and adaptive scientific computing. part I: Abstract

582 framework. *Computing*, 82(2–3):103–119.

583 Beudez, N., Doussan, C., Lefeuve-Mesgouez, G., and Mesgouez, A. (2013). Influence of three root

584 spatial arrangement on soil water flow and uptake. results from an explicit and an equivalent,

585 upscaled, model. *Procedia Environmental Sciences*, 19:37–46.

586 Bidel, L. P. R., Pagès, L., Rivière, L. M., Pelloux, G., and Lorendeau, J. Y. (2000). MassFlowDyn

587 I: A carbon transport and partitioning model for root system architecture. *Annals of Botany*,

588 85(6):869–886.

589 Daly, K. R., Tracy, S. R., Crout, N. M., Mairhofer, S., Pridmore, T. P., Mooney, S. J., and Roosevelt, T. (2017). Quantification of root water uptake in soil using X-ray computed tomography and

590 image-based modelling. *Plant, Cell & Environment*, 41:121–133.

591 image-based modelling. *Plant, Cell & Environment*, 41:121–133.

592 Delory, B. B. M., Li, M., Topp, C. C. N., and Lobet, G. (2018). archiDART v3.0: A new data

593 analysis pipeline allowing the topological analysis of plant root systems. *F1000Research*, 7:22.

594 Delory, B. M., Baudson, C., Brostaux, Y., Lobet, G., du Jardin, P., Pagès, L., and Delaplace, P.

595 (2016). archiDART: An R package for the automated computation of plant root architectural

596 traits. *Plant and Soil*, 398(1):351–365.

597 Doussan, C., Pagès, L., and Vercambre, G. (1998). Modelling of the hydraulic architecture of root

598 systems: An integrated approach to water absorption—model description. *Annals of Botany*,

599 81(2):213–223.

600 Dunbabin, V. M., McDermott, S., and Bengough, A. G. (2006). Upscaling from rhizosphere to

601 whole root system: Modelling the effects of phospholipid surfactants on water and nutrient

602 uptake. *Plant and Soil*, 283:57–72.

603 Dunbabin, V. M., Postma, J. A., Schnepf, A., Pagès, L., Javaux, M., Wu, L., Leitner, D.,

604 Chen, Y. L., Rengel, Z., and Diggle, A. J. (2013). Modelling root-soil interactions using three-

605 dimensional models of root growth, architecture and function. *Plant and Soil*, 372(1/2):93–124.

606 Dupuy, L. X., Fourcaud, T., Lac, P., and Stokes, A. (2007). A generic 3D finite element model of

607 tree anchorage integrating soil mechanics and real root system architecture. *American Journal*

608 *of Botany*, 94:1506–1514.

609 Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B., Stouffer, R. J., and Taylor, K. E.

610 (2016). Overview of the coupled model intercomparison project phase 6 (CMIP6) experimental

611 design and organization. *Geoscientific Model Development*, 9(5):1937–1958.

612 Fitter, A. H. (1987). An architectural approach to the comparative ecology of plant root systems.

613 *New Phytologist*, 106(s1):61–77.

614 Fitter, A. H. and Stickland, T. R. (1991). Architectural analysis of plant root systems 2. Influence

615 of nutrient supply on architecture in contrasting plant species. *New Phytologist*, 118(3):383–389.
616 Flemisch, B., Darcis, M., Erbertseder, K., Faigle, B., Lauser, A., Mosthaf, K., Müthing, S., Nuske, P., Tatomir, A., Wolff, M., and Helmig, R. (2011). DuMuₓ: DUNE for multi-{phase, component, scale, physics,...} flow and transport in porous media. *Advances in Water Resources*, 34(9):1102-1112.

620 Geuzaine, C. and Remacle, J.-F. (2009). Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. *International Journal for Numerical Methods in Engineering*, 79(11):1309-1331.

623 Helmig, R. et al. (1997). *Multiphase Flow and Transport Processes in the Subsurface: A Contribution to the Modeling of Hydrosystems*. Springer-Verlag.

626 Hund, A., Trachsel, S., and Stamp, P. (2009). Growth of axile and lateral roots of maize: I development of a phenotyping platform. *Plant and Soil*, 325(1):335-349.

627 Janssen, P. and Heuberger, P. (1995). Calibration of process-oriented models. *Ecological Modelling*, 83(1):55-66. Modelling Water, Carbon and Nutrient Cycles in Forests.

629 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 Journal*, 7(3):1079-1088.

632 Koch, T. (2019). Benchmark C1.2 - Numerical results reference solution. doi:10.18419/darus-471.

633 Koch, T., Gläser, D., Weishaupt, K., Ackermann, S., Beck, M., Becker, B., Burbulla, S., Class, H., Coltman, E., Emmert, S., Fetzer, T., Grüninger, C., Heck, K., Hommel, J., Kurz, T., Lipp, M., Mohammadi, F., Scherrer, S., Schneider, M., Seitz, G., Stadler, L., Utz, M., Weinhardt, F., and Flemisch, B. (2019). DuMuₓ 3 – an open-source simulator for solving flow and transport problems in porous media with a focus on model coupling. *arXiv e-prints*, page arXiv:1909.05052.

638 Koch, T., Gläser, D., Weishaupt, K., Ackermann, S., Beck, M., Becker, B., Burbulla, S., Class, H., Coltman, E., Fetzer, T., Flemisch, B., Grüninger, C., Heck, K., Hommel, J., Kurz, T., Lipp, M., Mohammadi, F., Schneider, M., Seitz, G., Scholz, S., and Weinhardt, F. (2018a). Dumux 3.0.0.

641 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 Journal*, 17.

644 Landl, M., Schnepf, A., Vanderborght, J., Bengough, A. G., Bauke, S. L., Lobet, G., Bol, R., and Vereecken, H. (2018). Measuring root system traits of wheat in 2d images to parameterize 3d root architecture models. *Plant and Soil*, 425(1):457-477.

647 Landsberg, J. and Fowkes, N. (1978). Water movement through plant roots. *Annals of Botany*, 42(3):493-508.

649 Li, M., Duncan, K., Topp, C. N., and Chitwood, D. H. (2017). Persistent homology and the branching topologies of plants. *American Journal of Botany*, 104(3):1-5.

651 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 Physiology*, 177:1382-1395.

654 Lobet, G., Pagès, L., and Draye, X. (2011). A novel image-analysis toolbox enabling quantitative analysis of root system architecture. *Plant Physiology*, 157:29-39.

656 Lobet, G., Pound, M. P., Diener, J., Pradal, C., Draye, X., Godin, C., Javaux, M., Leitner, D., Meunier, F., Nacry, P., Pridmore, T. P., and Schnepf, A. (2015). Root System Markup Language: Toward a unified root architecture description language. *Plant Physiology*, 167:617-627.

659 Luo, Y. Q., Randerson, J. T., Abramowitz, G., Bacour, C., Blyth, E., Carvalhais, N., Ciais, P., Dalmonech, D., Fisher, J. B., Fisher, R., Friedlingstein, P., Hibbard, K., Hoffman, F., Huntzinger, D., Jones, C. D., Koven, C., Lawrence, D., Li, D. J., Mahecha, M., Niu, S. L., Norby, R., Piao, S. L., Qi, X., Peylin, P., Prentice, I. C., Riley, W., Reichstein, M., Schwalm, C., Wang, Y. P., Xia, J. Y., Zaehle, S., and Zhou, X. H. (2012). A framework for benchmarking land models. *Biogeosciences*, 9:3857-3874.
665 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 and Soil*, 439:273-292.

668 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. *Applied Mathematical Modelling*, 52:648-663.

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

674 Oram, N. J., Ravenek, J. M., Barry, K., Weigelt, A., Chen, H., Gessler, A., Gockele, A., de Kroon, H., van der Paauw, J. W., Scherer-Lorenzen, M., Smit-Tiekstra, A., van Ruijven, J., and Mommer, L. (2018). Below-ground complementarity effects in a grassland biodiversity experiment are related to deep-rooting species. *Journal of Ecology*, 106:265-277.

678 Porter, C. H., Villalobos, C., Holzworth, D., Nelson, R., White, J. W., Athanasiadis, I. N., Janssen, S., Ripoche, D., Cufi, J., Raes, D., Zhang, M., Knapen, R., Sahajpal, R., Boote, K., and Jones, J. W. (2014). Harmonization and translation of crop modeling data to ensure interoperability. *Environmental Modelling & Software*, 62:495-508.

682 Roosevelt, T. and Fowler, A. (2004). A model for water uptake by plant roots. *Journal of Theoretical Biology*, 228(2):155-171.

684 Schnepf, A., Huber, K., Landl, M., Meunier, F., Petrich, L., and Schmidt, V. (2018). Statistical characterization of the root system architecture model CRootBox. *Vadose Zone Journal*, 17(1).

686 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 Journal*, 7.

689 Schroeder, W., Martin, K., and Lorensen, B. (2006). *The Visualization Toolkit*. Kitware, 4th edition.

691 Steefel, C. I., Yabusaki, S. B., and Mayer, K. U. (2015). Reactive transport benchmarks for subsurface environmental simulation. *Computational Geosciences*, 19(3):439-443.

693 The CGAL Project (2019). *CGAL User and Reference Manual*. CGAL Editorial Board, 4.14 edition.

695 Van Genuchten, M. (1980). Closed-form equation for predicting the hydraulic conductivity of unsaturated soils. *Soil Science Society of America Journal*, 44(5):892-898.

697 Vanderborght, J., Kasteel, R., Herbst, M., Javaux, M., Thiéry, D., Vanblooster, M., Mouvet, C., and Vereecken, H. (2005). A set of analytical benchmarks to test numerical models of flow and transport in soils. *Vadose Zone Journal*, 4(1):206-221.

700 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. *Annals of Botany*, 118(4):853-864.

## 703 A Derivation of the analytical solution of water flow inside the root system

705 The axial water flow rate in the xylem $Q_x$ (cm$^3$ day$^{-1}$) is given by

$$Q_x = -k_x \left( \frac{\partial \psi_x}{\partial \zeta} + v \cdot \mathbf{e}_3 \right), \qquad (19)$$

706 where $k_x$ is the axial conductance (cm$^3$ day$^{-1}$), $\psi_x$ is the pressure inside the xylem (cm), $\zeta$ is the local axial coordinate $\mathbf{e}_3$ the unit vector in z-direction, and $v$ the normalised direction of the xylem.
709 The radial water flow rate is given by

$$Q_r = -2r_{\text{root}}\pi l_{\text{seg}}k_r(\psi_s - \psi_x) \qquad (20)$$

710 with units (cm³ day⁻¹), where $r_{\text{root}}$ is the root radius (cm), $l_{\text{seg}}$ is the length of each root segment (cm), $k_r$ is the radial conductivity (day⁻¹), and $\psi_s$ is the soil water pressure head of the surrounding soil (cm). The equation is neglecting osmotic potential and is based on Eq. (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.

715 For each segment of length $l_{\text{seg}}$ mass conservation yields

$$0 = Q_x|_{l_{\text{seg}}} - Q_x|_0 + Q_r \qquad (21)$$

$$-\frac{1}{l_{\text{seg}}}Q_r = -\frac{1}{l_{\text{seg}}}(Q_x|_0 - Q_x|_{l_{\text{seg}}}) \quad \text{and for } l_{\text{seg}} \to 0 \qquad (22)$$

$$2r_{\text{root}}\pi k_r(\psi_s - \psi_x) = -k_x \frac{\partial^2 \psi_x}{\partial \zeta^2} \qquad (23)$$

716 see Eq. (3.4) of Roose and Fowler (2004), where $v_3$ is the z-component of the normalised xylem direction (cm).

717 Integrating this ordinary differential equation leads to an explicit equation for $\psi_x(\zeta)$

$$\psi_x(\zeta) = \psi_s + d_1 e^{\sqrt{c}\zeta} + d_2 e^{-\sqrt{c}\zeta}, \qquad (24)$$

718 where $c := 2a\pi k_r/k_x$, and $d_1$ and $d_2$ are integration constants that are derived from the boundary conditions.

720 To exemplify, we calculate $d_1$, and $d_2$ 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 $\psi_x|_{\text{collar}} = \psi_0$ is inserted into the analytic solution Eq. (24), and yields

$$\psi_s + d_1 + d_2 = \psi_0. \qquad (25)$$

724 The Neumann boundary condition $Q_x|_{l_{\text{seg}}} = 0$ (Eq. 20) leads to

$$\frac{\partial \psi_x}{\partial \zeta}|_{l_{\text{seg}}} = v_3, \qquad (26)$$

725 where $l_{\text{seg}}$ is the length of the root segment. Using the derivation of the analytical solution yields

$$d_1 \sqrt{c} e^{\sqrt{c} l_{\text{seg}}} - d_2 \sqrt{c} e^{-\sqrt{c} l_{\text{seg}}} = v_3. \qquad (27)$$

For a straight downward segment $v_3 = -1$, Eqns (25) and (27) can be summarized as

$$\begin{pmatrix} 1 & 1 \\ \sqrt{c}e^{\sqrt{c}l_{\text{seg}}} & -\sqrt{c}e^{-\sqrt{c}l_{\text{seg}}} \end{pmatrix} \begin{pmatrix} d_1 \\ d_2 \end{pmatrix} = \begin{pmatrix} \psi_0 - \psi_s \\ -1 \end{pmatrix} \qquad (28)$$

Solving this linear equation for $d_1$ and $d_2$ yields

$$d_1 = d^{-1} \left( e^{-\sqrt{c}l_{\text{seg}}} (\psi_0 - \psi_s) + 1 \right) \qquad (29)$$

$$d_2 = -d^{-1} \left( e^{\sqrt{c}l_{\text{seg}}} (\psi_0 - \psi_s) + 1 \right), \qquad (30)$$

where $d$ is the determinant of above matrix

$$d = e^{-\sqrt{c}l_{\text{seg}}} - e^{\sqrt{c}l_{\text{seg}}}. \qquad (31)$$