# Collaborative benchmarking of functional-structural root architecture models: Quantitative comparison of simulated root water uptake

Andrea Schnepf, Christopher K Black, Valentin Couvreur, Benjamin M Delory, Claude Doussan, Adrien Heymans, Mathieu Javaux, Deepanshu Khare, Axelle Koch, Timo Koch, 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.. Collaborative benchmarking of functional-structural root architecture models: Quantitative comparison of simulated root water uptake. in silico Plants, 2023, 5 (1), pp.1-21. 10.1093/insili-coplants/diad005. hal-04191593

HAL Id: hal-04191593

https://hal.inrae.fr/hal-04191593

Submitted on 30 Aug 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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>
Original Article

Collaborative benchmarking of functional-structural root
architecture models: Quantitative comparison of simulated
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>, Adrien Heymans<sup>3</sup>, Mathieu Javaux<sup>1,6</sup>, Deepanshu Khare<sup>1,14</sup>, Axelle Koch<sup>6</sup>, Timo Koch<sup>7,8,9</sup>, Christian W. Kuppe<sup>10</sup>, Magdalena Landl<sup>1,14</sup>, Daniel Leitner<sup>1,14</sup>, Guillaume Lobet<sup>1,14</sup>, Félicien Meunier<sup>11,12</sup>, Johannes A. Postma<sup>10</sup>, Ernst D. Schäfer<sup>2,13</sup>, Tobias Selzner<sup>1,14</sup>, Jan Vanderborght<sup>1,14</sup> and Harry Vereecken<sup>1,14</sup>

ABSTRACT

Functional-structural root architecture models have evolved as tools for the design of improved agricultural management practices and for the selection of optimal root traits. In order to test their accuracy and reliability, we present the first benchmarking of root water uptake from soil using five well-established functional-structural root architecture models: DuMux, CPlantBox, R-SWMS, OpenSimRoot and SRI. The benchmark scenarios include basic tests for water flow in soil and roots as well as advanced tests for the coupled soil-root system. The reference solutions and the solutions of the different simulators are available through Jupyter Notebooks on a GitHub repository. All of the simulators were able to pass the basic tests and continued to perform well in the benchmarks for the coupled soil-plant system. For the advanced tests, we created an overview of the different ways of coupling the soil and the root domains as well as the different methods used to account for rhizosphere resistance to water flow. Although the methods used for coupling and modelling rhizosphere resistance were quite different, all simulators were in reasonably good agreement with the reference solution. During this benchmarking effort, individual simulators were able to learn about their strengths and challenges, while some were even able to improve their code. Some now include the benchmarks as standard tests within their codes. Additional model results may be added to the GitHub repository at any point in the future and will be automatically included in the comparison.

**KEYWORDS:** Functional-structural root architecture models; model comparison; benchmark; root water uptake; quantitative comparison.

<footnote n="1"><sup>1</sup>Forschungszentrum Jülich GmbH, Institute of Bio- and Geosciences – Agrosphere (IBG-3), 52425 Jülich, Germany.</footnote>

<footnote n="2">²Department of Plant Science, The Pennsylvania State University, 102 Tyson Building, University Park PA 16802, USA.</footnote>

<footnote n="3">³Earth and Life Institute, Agronomy, Université catholique de Louvain, Louvain-la-Neuve, Belgium</footnote>

<footnote n="4">⁴Institute of Ecology, Leuphana University Lüneburg, Universitätsallee 1, 21335 Lüneburg, Germany.</footnote>

<footnote n="5">⁵INRAE, Avignon Université, EMMAH, F-84000 Avignon, France.</footnote>

<footnote n="6">⁶Earth and Life Institute, Environmental Sciences, Université catholique de Louvain, Louvain-la-Neuve, Belgium</footnote>

<footnote n="7">⁷Department of Mathematics, University of Oslo, Postboks 1053, Blindern, 0316 Oslo, Norway</footnote>

<footnote n="8">⁸SCAN department, Simula Research Laboratory, Kristian Augusts gate 23, 0164 Oslo, Norway</footnote>

<footnote n="9">⁹Department of Hydromechanics and Modelling of Hydrosystems, University of Stuttgart, Pfaffenwaldring 61,
70569 Stuttgart, Germany.</footnote>

<footnote n="10">¹⁰Forschungszentrum Jülich GmbH, Institute of Bio- and Geosciences – Plant Sciences (IBG-2), 52425 Jülich, Germany.</footnote>

<footnote n="11">¹¹CAVELab - Computational and Applied Vegetation Ecology, Ghent University, Ghent, Belgium.</footnote>

<footnote n="12">¹²Department of Earth and Environment, Boston University, Boston, USA.</footnote>

<footnote n="13">¹³School of Mathematical Sciences, University of Nottingham, Mathematical Sciences Building, University Park, Nottingham, NG7 2RD, United Kingdom</footnote>

<footnote n="14">¹⁴International Soil Modelling Consortium ISMC, Jülich, Germany</footnote>

<footnote>*Corresponding author's e-mail address: Andrea Schnepf, Forschungszentrum Jülich GmbH, Institute of Bio- and Geosciences - Agrosphere (IBG-3), 52425 Jülich, Germany; E-mail: a.schnepf@fz-juelich.de</footnote>

<footnote>**Citation:** Schnepf A, Black CK, Couvreur V, Delory BM, Doussan C, Heymans A, Javaux M, Khare D, Koch A, Koch T, Kuppe CW, Landl M, Leitner D, Lobet G, Meunier F, Postma JA, Schäfer ED, Selzner T, Vanderborght J, Vereecken H. 2023. Collaborative benchmarking of functional-structural root architecture models: Quantitative comparison of simulated root water uptake. *In Silico Plants* 2023: diad005; doi: 10.1093/insilicoplants/diad005</footnote>

<footnote>Handling editor: Graeme Hammer</footnote>
# 1. INTRODUCTION

This work responds to the call to participate in the collaborative benchmarking of functional-structural root architecture models (Schnepf et al., 2020). The need for this activity was outlined in Schnepf et al. (2022). Functional-structural root architecture models (FSRMs) have evolved as tools for the design of agricultural management schemes for improved resource efficiency and for the selection of root traits to optimize plant performance in specific environments. Although reliable applications are crucial, this type of models has never been benchmarked. Successful benchmark studies have been conducted in other fields, such as in crop modelling (Agricultural Model Intercomparison and Improvement Project AGMIP, Rosenzweig et al., 2014), reactive transport modelling (Steefel et al., 2015), or modelling of soil water flow and solute transport (Vanderborght et al., 2005). This work sets out to gain an overview of the differences between the outputs of several FSRMs for the benchmark problems described in Schnepf et al. (2020). It also sought to gain an understanding of the underlying reasons, be it differences in the mathematical formulation of the processes and their coupling, in the numerical scheme, or even from coding errors.

The call included benchmark problems for root growth models, soil water flow models, root water flow models, and for water flow in the coupled soil-root system. All the benchmarks and corresponding reference solutions were published in the form of Jupyter Notebooks on the GitHub repository https://github.com/RSA-benchmarks/collaborative-comparison. Fig. 1 shows an overview of the different benchmark problems described in Schnepf et al. (2020). In particular, we distinguish between individual modules that solve flow and transport problems in one spatial domain only, soil or roots, and the coupled scenarios, where flow and transport problems are solved in both domains and include an exchange between domains, which is an important feature of functional-structural root architecture models.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** The benchmark problems described in Schnepf et al. (2020) are divided into five categories: M1 focuses on root growth, M2 on water flow in soil, M3 on water flow in root, C1 on the coupled root-soil system with a static root architecture, and C2 on the coupled root-soil system with a growing root architecture. The categories marked in green are addressed in this paper.</figcaption></figure>

Here, we address the individual modules M2 (water flow in soil) and M3 (water flow in roots), as well as the coupled module C1, i.e., water uptake from a dynamic soil by a static, predefined, root architecture. Thus, three folders within the repository are relevant for this manuscript, “M2 Water flow in soil”, “M3 Water flow in roots” and “C1 Coupled problem, static RSA”. For each scenario, they contain subfolders that follow the same structure. They contain folders named “Numerical results” where each participating simulator could upload their results, a Jupyter Notebook “Benchmark problem.ipynb” that describes the benchmark problem and its reference solution, as well as a Jupyter Notebook “Automated comparison.ipynb” that automatically loads the numerical results and creates plots for comparing simulator results to the reference solution. Where appropriate, additional folders contain the prescribed root grid files.

The simulation domains and sample results for each benchmark scenario is illustrated in Fig. 2. Growing root architectures are addressed in a separate paper by Delory et al. (in prep.). The focus of this paper is to evaluate the effect of various methods used by the different simulators of coupling the soil and the root domains on simulated root water uptake.

Five groups that develop functional-structural root architecture models have contributed to the benchmark problems with their solutions: DuMuₓ (Koch et al., 2018, 2021), CPlantBox-DuMuₓ (coupling CPlantBox, Schnepf et al. (2018); Zhou et al. (2020) and DuMuₓ, Mai et al. (2019)), OPENSIMROOT (Postma et al., 2017), R-SWMS (Javaux et al., 2008) and SRI (Beudez et al., 2013).

# 2. METHODS

Throughout this paper, the term “model” refers to the mathematical equations and parameters that describe the system of interest. The term “simulator” refers to the software code that is an implementation of the numerical schemes used to solve the described model. A brief description of the benchmark problems is given in Section 2.1. For the benchmarks on water flow in either soil or roots, M2 and M3, the simulators solved the same mathematical model but with different numerical methods and software implementations. The mathematical equation in the case of water flow in soil is the Richards equation, while water flow in roots is described by a Darcy-type axial flow model with a root water uptake source term that is proportional to the pressure head difference between the xylem and the root surface. In the benchmarks of water flow in the coupled soil-root system, C1.1 and C1.2, each simulator makes a distinct choice concerning the coupling of the different submodels, specifically the coupling of soil and root domains. The simulators’ sufficiently accurate solutions to the subproblems of water flow in soil, M2, and water flow in roots, M3, were required to interpret the solutions of the coupled problem of root water uptake from a drying soil by a static root architecture, C1.2.

A detailed description of the different simulators is given in Section A and an overview of their main characteristic features is given in Table 2. In addition, each simulator has a different approach to calculate the sink term for root water uptake that results from coupling the two subproblems of soil and root water
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Illustration of the simulation domains and results of the different benchmark problems: (a) Simulated infiltration front in benchmark scenario M2.1 infiltration, (b) simulated soil drying due to evaporation in M2.2 evaporation, simulated xylem water pressure head within a non-growing (c) single root in M3.1 and (d) root system in M3.2, both embedded in a soil with constant soil water pressure head, (e) soil water pressure head as a result of radial water flow towards a root in C1.1, (f) results of the coupled soil-root problem in C1.2 for a root system embedded in a dynamic soil</figcaption></figure>

flow. A detailed description of the available sink term approaches for the different simulators is shown in Fig. 3.

## 2.1. Brief description of the benchmark problems

The benchmark problems, which are elaborated here, combine the communities of soil physics and plant physiology, which commonly follow different standards regarding the dimensions used for the water potential. Throughout this manuscript, we express the water potentials within both the soil and the xylem in dimensions that are common in soil physics, i.e. pressure head (L). Water potential can be converted from dimensions of head to Joule $M^{-1}$ through multiplication with the acceleration constant due to gravity, $g$ (L $T^{-2}$). A further multiplication with the density of water, $\rho_w$ (M L$^{-3}$), results in dimensions of pressure, P. By convention, we consider the atmospheric potential as reference, equal to zero.

### 2.1.1. Soil subproblems

These sets of benchmark problems consider infiltration and evaporation scenarios (M2.1 and M2.2) that were originally described in Vanderborght et al. (2005).

For three soils (sand, loam, and clay; see Table 1), the infiltration scenario describes the water flow into initially dry soil. Water flow rates in dry soil are low. When the soil surface is kept at a high saturation level, a sharp front forms between dry and wet soil and moves downward into the dry soil. The different simulators implemented this problem for a box-shaped soil domain with a depth of 200 cm and a length and width of 10 cm. The problem was solved as quasi-1D and each simulator was free to choose the resolution in the z-direction. Fingering or hysteresis was not considered. The main interest of this scenario is the correct shape and velocity of the infiltration front. A visualization of the problem is shown in Fig. 2 (a).

We used the same soil types in the evaporation benchmark scenario M2.2. When the soil is unable to support the evaporation rate prescribed at the top boundary of the soil domain, the boundary condition switches from a Neumann boundary condition, where we define the water flux, to a Dirichlet boundary condition where the soil water pressure head is maintained at a critical soil water pressure head of -10,000 cm. In the related benchmark problem, we are interested in the transition times between sink-limited evaporation (evaporation dictated by the atmospheric boundary) and source-limited evaporation (evaporation dictated by water transfer in soil), given the constant potential
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Overview of the different approaches to computing the sink terms for root water uptake used by the different simulators in comparison to the approach used by the reference solution. **Left panel:** The reference solution solves the Richards equation in a 3D soil domain in which the root volume is explicitly considered, i.e. the root domain is cut out of the soil domain. Root water uptake is computed as a boundary condition at the soil-interface (written in blue). The numerical grid has to be highly refined near the root surfaces. **Middle panel:** The simulators disregard the root volume in the soil domain. Root water uptake is computed as sink term of the Richards equation (indicated in blue). The numerical grid may be coarse compared to the root diameters; this significantly reduces the computational effort involved. **Right panel (top):** Methods of distributing the root water uptake source term. At the continuum level, the root water uptake is either considered as a line, surface, or volume source. Delta or kernel functions are used to distribute the source terms to the centerline, the root surface, or a tubular support region around the root. At the discrete levels, they are then distributed to those df (degrees of freedom) that correspond to the respective support regions (here denoted in blue). **Right panel (bottom):** Root water uptake is computed according to the water pressure head difference between the xylem and the soil water pressure head at the soil-root interface, $\psi_{s,r}$. The different methods include i) averaging the df of the soil element in which the root is located, ii) interpolating the corresponding df to evaluate the values at certain positions, such as the centerline or the root surface, or iii) approximating the rhizosphere gradient. The latter method enables the estimation of the rhizosphere gradients in the order of mm in a coarser (cm scale) soil discretisation. The root and soil illustrations in the left and middle panel are taken from Koch et al. (2022) (license: CC-BY-4.0).</figcaption></figure>

<table><caption>**Table 1.** Soil hydraulic properties taken from Vanderborght et al. (2005); Schnepf et al. (2020), according to the van Genuchten and Mualem-van Genuchten equations ($\theta_r, \theta_s$: residual and saturated volumetric water contents, $\alpha, n, \lambda$: empirical parameters, $\psi_s$: soil water pressure head, $K_s$: saturated hydraulic conductivity)</caption><tr><td></td><td>θ<sub>r</sub></td><td>θ<sub>s</sub></td><td>α</td><td>n</td><td>K<sub>s</sub></td><td>λ</td></tr><tr><td>Soil type</td><td>(-)</td><td>(-)</td><td>(cm<sup>-1</sup>)</td><td>(-)</td><td>(cm d<sup>-1</sup>)</td><td>(-)</td></tr><tr><td>Sand</td><td>0.045</td><td>0.43</td><td>0.15</td><td>3.0</td><td>1,000</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></table>

evaporation rate. A visualization of the problem is shown in Fig. 2 (b).

For M2.1 and M2.2, the reference solutions are analytical equations described in Vanderborght et al. (2005); Schnepf et al. (2020). It should be noted that the analytical solution for the evaporation scenario does not account for gravity. Disregarding

gravity for vertical flow is an option for some - but not all - of the simulators. However, gravity may be assumed to have a negligible influence compared to the large differences in soil water pressure head.

### 2.1.2. Root subproblems

The benchmarks for the root subproblems concern water flow within roots. In this exercise, it is assumed that the medium in which the roots are embedded does not exhibit any resistance to water flow or dynamic changes in water content (such as in hydroponics). This allows for a static pressure head to be defined for this medium. Benchmark scenario M3.1 describes the root water pressure head distribution in the xylem of a single vertical root segment. A pressure head of -1,000 cm is defined as the boundary condition of the root collar, a no-flux boundary condition at the root tip, and a constant pressure head of -200 cm in the medium in which the root is embedded (Fig. 2 (c)). The analytical solution to this problem is given in the appendix of Schnepf et al. (2020). Its formulation is equivalent to that of Landsberg and Fowkes (1978); Meunier et al. (2017a) but uses exponential instead of hyperbolic functions and includes gravity. This benchmark provides a basic test to compare the root
water pressure head distributions along a single root axis for the different simulators.

The benchmark for water uptake by a small root system from a static soil (M3.2) describes root water flow within a prescribed root system of a lupine plant that was manually traced in an MRI (magnetic resonance imaging) image (Fig. 2 (d)). The root architecture is provided in the root system markup language (RSML, Lobet et al., 2015) on the corresponding GitHub repository<fnref n="1" />. It consists of 2,883 segments with lengths ranging from 0.5 mm to 1 mm (mean value 0.8 mm). With such a high spatial resolution, simulators are expected to achieve results close to the reference solution (Meunier et al., 2017b). As for M3.1, the pressure head of the medium outside the roots is predefined at -200 cm. Two scenarios are considered: one in which the root hydraulic properties are constant (M.2a, with root axial conductance $k_x = 4.32 \times 10^{-2}$ cm³ d⁻<fnref n="1" /> and root radial conductivity $k_r = 1.728 \times 10^{-4}$ d⁻<fnref n="1" />) and one in which the root hydraulic properties are age-dependent and different for each root branching order (M3.2b, as shown in Fig. S1). The reference solutions to benchmarks M3.2a,b are computed using the hybrid analytical-numerical approach given by Meunier et al. (2017b), which is described in detail in Schnepf et al. (2020). The derivation of this solution using exponential instead of hyperbolic functions is described in appendix B. The main focus of this benchmark scenario is to compare root water pressure head distributions within the root system.

### 2.1.3. Coupled problems including root-soil interaction

The benchmark scenario C1.1 follows the paper of Schröder et al. (2008). We refer to Schnepf et al. (2020) for a more detailed benchmark description. In brief, it considers the 1-dimensional radially symmetric soil water pressure head distribution around a single root, which takes up water at a prescribed rate and switching to a constant soil water pressure head at the root boundary if the critical value of -15,000 cm is reached (Fig. 2 (e)). Based on the steady-rate assumption, the reference solution is an approximate analytical solution to this problem (Schröder et al., 2008). With this benchmark, we want to test if the different simulators are able to reproduce the soil water pressure head gradients that can occur in dry soil between the bulk soil and the root surface. This pressure head is relevant for upscaling to the root system scale (see C1.2). In this regard, the different simulators still use their volumetric sink term approach for computing root water uptake from soil. They therefore consider a single straight root in the center of a 3-dimensional soil domain. Only horizontal flow between soil elements is considered, meaning that gravity can be disabled in the simulators. The constant root uptake is described as a volumetric sink term in the center of the domain. The horizontal resolution of the soil domain needs to be fine enough to capture the gradients of soil water pressure heads. The simulations start at an initial homogeneous soil water pressure head of -100 cm and are stopped at the onset of water stress, i.e., when the pressure head at the root xylem upper boundary is -15,000 cm. The time of stress onset and the gradients of soil

water pressure heads, are then compared. The gradients of soil water pressure heads from the 3D solution are sampled along a transect normal to the root surface.

Benchmark C1.2 is the most complex, but essential benchmark for this FSRM benchmarking. It represents the water uptake from a drying soil by a static (i.e. constant) root system architecture (see Schnepf et al., 2020). The root architecture is provided in the RSML format on the corresponding GitHub repository<fnref n="2" />. It is a younger version of the root system of M3.2 and consists of only 584 segments with lengths ranging from 0.2 mm to 1.8 mm (mean value: 0.9 mm). The root hydraulic properties are either constant (C1.2a) or vary with root type and age (C1.2b) with the same root hydraulic properties as in M3.2a,b. To challenge the different simulators, we chose an initially dry soil that becomes even drier as water is taken up by the roots, with large gradients of soil water pressure heads expected around the roots. Transpiration is prescribed by a sinusoidal function, with maximum transpiration at noon and no transpiration at midnight. In this scenario, we are interested in comparing actual transpiration and cumulative root water uptake over time, particularly during times of stress onset, as well as the root sink term depth profiles. A visualization of the problem is shown in Fig. 2 (f).

The simulators differ in terms of how soil and root problems are coupled. The approaches used by the different simulators are outlined in Fig. 3 and are described in detail in the following paragraphs.

Since an analytical solution is not available for this benchmark, our reference solution is numerical and already described by (Schnepf et al., 2020; Koch, 2022). In brief, the reference solution is calculated using a soil domain that explicitly considers the roots (i.e. the roots are cut out of the soil domain), and in which the root water uptake is not a volumetric sink term but a (Robin-type, i.e., a linear combination of Neumann and Dirichlet) boundary condition for the Richards equation on the root-soil interface $\Gamma$,

$$ (K(\theta)(\nabla\psi_s + \nabla z)) \cdot n_{\Gamma} = -k_r[\psi_s(x_{\Gamma}) - \psi_x(\Pi x_{\Gamma})], \quad (1) $$

where $\psi_s$ and $\psi_x$ are the water pressure heads in soil and roots, $K(\theta)$ is the unsaturated hydraulic conductivity, $k_r$ is the root radial conductivity, $n_{\Gamma}$ is the unit normal vector on $\Gamma$ pointing from root surface to the soil and $\Pi$ is a surjective mapping of every root-soil surface point to a unique point on the root centerline network (Koch, 2022).

The numerical solution requires a high-resolution mesh of the soil domain (excluding the roots) that is refined near the root-soil boundary and fully coupled to the 1D-model root solution at the interface. We computed this with a DuMu$^x$ implementation, which we refer to as “DuMu$^x$_explicit_interface”.

In this approach, the root-soil interface was explicitly resolved by the computational mesh in the soil domain. The root domain was implicitly defined as a volume, created by moving a ball with a radius $R$ (where $R$ is the constant local radius field per

<footnote n="1">¹ https://github.com/RSA-benchmarks/collaborative-comparison/M3Waterflowinroots/M3.2Rootsystem</footnote>

<footnote n="2">² https://github.com/RSA-benchmarks/collaborative-comparison/C1Coupledproblem, staticRSA/C1.2Rootsystem/root_grid</footnote>
root segment) along the root segment centerlines. The root-soil interface is therefore implicitly given by signed distance functions. The 3D mesh based on this implicit function was created with the software CGAL (The CGAL Project, 2019). For C1.2, we used a tetrahedral mesh with a strong local grid refinement towards the root-soil interface and around 8.3 million elements. In comparison to Schnepf et al. (2020), the mesh was further refined to decrease the discretisation error (the smallest 10% of elements had a mean size of 10 µm). We discretised Darcy's equation in the root domain using a cell-centered finite volume method with two-point flux approximation (TPFA) and discretised the Richards equation in the soil using a vertex-centered finite volume method (Box). The simulation for C1.2 used 256 adaptive time steps; the time-stepping criterion was based on the number of Newton iterations (the smallest and largest time step sizes were 300 s and 1200 s, respectively). The runtime was approximately 31 h on a multi-core machine with 64 cores (the assembly procedure makes use of shared-memory parallelism).

## 2.2. Summary of the approaches of the different simulators

In this section we summarise the numerical and coupling approaches used by the different simulators while the details are described Appendix A. The basic conditions for the implementation of the different benchmark problems included free choice of numerical scheme and spatial discretisation for each simulator. For benchmark scenarios C1, where additional differences between simulators are due to the approaches used for coupling the root and the soil domains, we did not prescribe the coupling approach to be used but let the modellers make their own choices.

Table 2 outlines the different numerical approaches used by the simulators for the soil-, root- and coupled problems. DuMu<sup>x</sup> solves the coupled problem in a monolithic way, in which the

coupled root-soil system is solved as a whole. The temporal resolution is performed at the level of the coupled problem using the implicit Euler method. All other simulators employ a sequential coupling, whereby the two subproblems are solved separately using distinct solvers. Information is exchanged only once in each time step, or with iterations that aim to converge towards the monolithic system.

For spatial discretisation of the soil subproblems either finite volumes (DuMu<sup>x</sup>) or finite elements (SRI, R-SWMS, OPENSIMROOT) are used. Those solvers that use sequential coupling may apply a temporal discretisation within each coupling time step, and rely on the implicit Euler method. Newton's method (DuMu<sup>x</sup>) or a fixed point iteration (R-SWMS, OPENSIMROOT, SRI) is used as nonlinear solver. A weighted Laplacian is used to describe the water flow in the tree graph that represents the root architecture. The corresponding edge fluxes are computed analytically (CPlantBox) or numerically (DuMu<sup>x</sup>, R-SWMS, OPENSIMROOT, SRI).

Fig. 3 sheds further light on the different approaches to computing the sink terms for root water uptake used by the different simulators in comparison to the approach used by the reference solution. The left panel illustrates the approach of our (numerical) reference solution. The volume of the roots is explicitly considered in the soil domain. This means cutting out long and thin cylindrical shapes and thus resulting in a highly resolved computational mesh around those structures. The root water uptake is computed via a boundary condition at the root-soil interfaces.

The middle panel illustrates the general approach used by the different simulators, where the volume of the roots is neglected in the soil domain and the root water uptake is computed via a sink term in the Richards equation. Mathematically, it is thus a different problem that the one solved by the reference solution. However, the aim is to produce results as close as possible to the reference solution while being able to use a coarser mesh

<table><caption>**Table 2.** Overview of the numerical solution methods of the different simulators (FV: finite volumes, FE: finite elements, FD: finite differences, TG-*: Tree graph, edge fluxes computed with method *)</caption><thead><tr><th>Simulators</th><th>DuMu<sup>x</sup></th><th>CPlantBox-DuMu<sup>x</sup></th><th>R-SWMS</th><th>OpenSimRoot</th><th>SRI</th></tr></thead><tbody><tr><td><strong>Coupled problem</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Coupling method</td><td>Monolithic</td><td>Sequential</td><td>Sequential</td><td>Sequential</td><td>Sequential</td></tr><tr><td>Time stepping</td><td>Implicit Euler</td><td>Explicit</td><td>Implicit</td><td>Implicit</td><td>Implicit</td></tr><tr><td>Sink term distribution</td><td>Cylinder surface or tube</td><td>Centerline</td><td>Centerline</td><td>Centerline</td><td>Centerline</td></tr><tr><td>Method for considering rhizosphere resistance</td><td>Perimeter average or centerline with interface value reconstruction by steady-state solution</td><td>Transient numerical solution of 1D radial rhizosphere models</td><td>Steady-rate approach</td><td>Adding soil hydraulic resistance to the root radial resistance in the RWU model</td><td>Superposition 3D analytical solution steady rate approach</td></tr><tr><td><strong>Soil subproblem</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Spatial discretisation</td><td>FV</td><td>FV</td><td>FE</td><td>FE</td><td>FE</td></tr><tr><td>Temporal discretisation within each sequential coupling time step</td><td>- (not applicable)</td><td>Implicit Euler</td><td>Implicit Euler</td><td>Implicit Euler</td><td>Implicit Euler</td></tr><tr><td>Nonlinear solver</td><td>Newton</td><td>Newton</td><td>Fixed point iteration</td><td>Fixed point iteration</td><td>Fixed point iteration</td></tr><tr><td><strong>Root subproblem</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Spatial discretisation</td><td>TG-FV</td><td>TG-analytical</td><td>TG-FD</td><td>TG-FD</td><td>TG-FD</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Numerical solutions of the infiltration benchmark (M2.1) compared to the analytical solution: Depth profiles of the volumetric water content, $\theta$, during infiltration in initially uniform dry soils.</figcaption></figure>

and thus benefit from less costly computations. The different approaches used by each simulator regarding how to distribute the sink terms for root water uptake in the soil domain as well as how to estimate the soil water pressure head at the root-soil interface are illustrated in the right panel. Furthermore, it also provides information which approach was used by which simulator for benchmark problem C1.2. In the small icon-like illustrations, the squares represent soil grid elements, the green dots represent the center of the root, the bold black circles stand for the root perimeter, the dashed/dotted black circles stands for the sphere of influence of the root and the red circles represent positions in the soil domain at which the soil water pressure heads are evaluated. The different methods of how to distribute the root water uptake sink term consider soil grid elements within which the roots are located represented as centerlines, the root cylinder surfaces, or a volumetric support region around each root are located. The approaches to determine the soil water pressure head at the soil-root interface include averaging the degrees of freedom corresponding to the soil grid element within which a root is located (e.g., taking the mean of the 8 nodes of a cubic grid or the one value of a cell-centered finite volume cell), interpolating the corresponding degrees of freedom to evaluate values at certain positions (e.g. the root centerline or surface), or methods to approximate subresolution rhizosphere gradients. For more details about the latter, please refer to Appendix A.

# 3. RESULTS

## 3.1. Infiltration and evaporation benchmarks

Fig. 4 shows the results of the benchmark scenario M2.1, where the infiltration fronts into an initially dry soil ($\psi_{s,i} = -400$ cm) are computed for the three different soil types (Table 1). The four simulators DuMu<sup>x</sup>, R-SWMS, OPENSIMROOT

and SRI each solve the benchmark sufficiently well, both in terms of the shape of the infiltration front as well as its temporal evolution. The only exception is SRI for loam and clay, where the front moves too slowly, resulting in deviations from the reference solution of approximately 3 cm in loam and 20 cm in clay after one-and-a-half days, respectively.

In the evaporation benchmark (M2.2), the four simulators perform well (Fig. 5). In the case of the loam 1 and clay scenarios, the transition time of stage-1 to stage-2 evaporation is slightly overestimated by all simulators in a similar way. At the top of the soil domain, the benchmark defines a switch in the boundary condition from Neumann to Dirichlet when the soil water potential reaches a given threshold. This might be a challenge for the numerical solution as it creates a discontinuity, and is also dependent on the grid resolution. Some of the simulators (DuMu<sup>x</sup>, OPENSIMROOT) offer an implementation of this boundary condition in a form where the transition from potential to actual transpiration is smoothed. A comparison for the effects of grid resolution and smooth transition is provided in Fig. S3 for the example of OPENSIMROOT. The smoothed version provided results that were closer to the benchmark, less sensitive to resolution, and a lot quicker to solve.

## 3.2. Water flow in roots

All four simulators compared well with the analytical solution in a basic test of water pressure distribution in a single vertical root (Fig. 6).

We then compared the root water pressure heads in the branched root structure (Fig. 7), comparing the results of the simulators to the reference solution. This reference solution, based on the hybrid analytical solution of Meunier et al. (2017b, 2022), is implemented within CPlantBox. The CPlantBox solution is thus equal to the reference solution.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Numerical solutions of the evaporation benchmark (M2.2) compared to the analytical solution: Temporal evolution of the actual soil evaporation, *E*<sub>act</sub>, for four different soils.</figcaption></figure>

The scenarios of constant and age-dependent root hydraulic properties differ strongly (Fig. 7 (a) and (b)). The root water pressure heads in the lateral roots is much less negative in the age-dependent scenario. In the case of constant root hydraulic properties, the lateral roots do not show strong gradients between their tip and base, and their root water pressure head is determined by the root water pressure head of the primary root at the branching point. In the case of age-dependent root hydraulic properties, the lateral roots show a significant root water pressure head gradient, and the root water pressure head at the tips is equal to a soil water pressure head of –200 cm. The results of the different simulators all compare well with the reference solution.

**3.3. Root water uptake by a static root system**

This section displays the results of the coupled benchmark problems C1 (SRI only took part in benchmark scenario C1.2a). Fig. 8 shows 1D radially symmetric soil water pressure head gradients around a single root at the onset of stress (i.e. the time

at which the soil water pressure head at the root surface reached –15,000 cm) as computed by all simulators in benchmark scenario C1.1. For both prescribed fluxes at the root-soil boundary, the sand scenario was immediately under stress, a fact recognised by all simulators. In the loam scenario, the numerical solutions of R-SWMS and OPENSIMROOT underestimate the soil water pressure head at the onset of stress, while CPlantBox-DuMuₓ and DuMuₓ_CYL overestimate it for the clay scenario. The difference between the two groups of simulators is that the former are based on finite element schemes while the latter is based on finite volume schemes. The results of DuMuₓ_CYL are only provided up to 0.5 cm. This is due to the fact that this simulator used a cuboid domain with equal volume of the cylindrical domain, with edge lengths of 1.06, 1.06 and 1 cm. The stress onset times differ slightly between simulators, but all are in the correct order of magnitude.

Root water uptake by a static (nongrowing) root system from a drying soil is predicted by the different simulators in
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Numerical solutions of the root water pressure head distribution over soil depth in a single vertical root compared to the analytical solution.</figcaption></figure>

benchmark C1.2. Fig. 9a and 9b show the results of the actual and cumulative actual transpiration for all simulators for both the cases with constant and age-dependent root hydraulic properties. Only when the simulators accounted for gradients in the rhizosphere (i.e. a solution that considers the drop in soil hydraulic conductivity near the roots) did the results come close to those of the reference solution. Otherwise, all simulators largely overestimated the transpiration compared to the reference solution. Different methods have been used to consider the hydraulic conductivity drop (see Fig. 3 for their hydraulic properties) and lead to significant improvements in the predicted actual transpiration. Some results slightly overestimated or underestimated the transpiration compared to the numerical reference solution, but all simulators reached the correct order of magnitude.

Although the distribution of root water pressure heads in the lateral roots is quite different between the scenarios with constant and age-dependent root hydraulic properties (see Fig. 7), the overall amount of transpired water is very similar.

Fig. 10a shows sink term depth profiles of root water uptake for constant root hydraulic properties. The left column shows the root water uptake at noon, while the right column shows soil water redistribution during the night. The sink term profiles for the age-dependent case are shown in Fig. 10b. This figure confirms the need to include the rhizosphere hydraulic conductivity drop to obtain smoother and lower variations as those obtained by the reference solution.

## 4. DISCUSSION

To our knowledge, this is the first benchmarking of functional-structural root architecture models in the context of root water uptake from soil. Regarding the individual modules, soil water flow models have been benchmarked before. The water flow in soil benchmarks were taken from Vanderborght et al. (2005) where solutions of different simulators were compared to

1D analytical solutions of the Richards equation. Further, more complex benchmarks for water flow in soil include e.g. dual permeability models (Bachmair et al., 2010) have not been considered here. In our benchmark activity, we learned that all simulators could well predict the results of the reference solution. As in Vanderborght et al. (2005), the largest discrepancies between analytical and numerical solutions of the infiltration scenario were found for the clay soil, and the differences between simulators were in a similar order of magnitude. This is due to the strongly nonlinear hydraulic functions in this scenario such that the infiltration fronts are very steep. Grid size, convergence criteria and method of evaluating the soil hydraulic conductivities strongly influence the results. As in Vanderborght et al. (2005), simulators tend to predict a slightly larger time period for the potential transpiration in the evaporation scenario. This error could be decreased by reducing the grid size, increasing the maximum number of iterations in the Richards equation solver or implementing the top boundary condition as a smooth function.

Models of water flow within roots have to our knowledge not yet been benchmarked before, except of course the of the comparison of the hybrid analytical solution of Meunier et al. (2017a). Meunier et al. (2022) showed that the appropriate grid size (root segment length) for the numerical solution is a function of the desired accuracy and the ratio of the root radial conductivity and the root axial conductance, and that a segment length of 1 cm could lead to substantial errors of 30%. As the mean segment lengths in our benchmark scenario M3.2 was 0.8 mm, all simulators obtained very accurate results, as can be expected. The comparison of rhizosphere models so far has had some focus on nutrient uptake or pH. Numerous authors compared rhizosphere model outputs to experimental data (e.g. Custos et al., 2020; Samal et al., 2010; Kirk, 1999; Kelly et al., 1992). Nowack et al. (2006) compared different available simulators in their ability to predict nutrient gradients in the rhizosphere as well as nutrient uptake. The reference solution there was the approximate analytical solution of Roosevelt et al. (2001). The benchmark scenario C1.1 now provides a similar comparison for water potential gradients and water uptake by a single root with the steady-rate approach of Schröder et al. (2008) as the reference solution.

There are several publications comparing root water uptake models, (e.g. de Jong van Lier et al., 2013; Cai et al., 2018), but they do not explicitly consider the root architecture. Dunbabin et al. (2013) provided a review about different simulators and Janott et al. (2011) a sensitivity analysis, for root architecture models, but they do not compare different simulators. de Willigen et al. (2012) compared root water uptake models of different complexity (1-, 2-, 3D), where one involved a root architecture model (R-SWMS). They found that the differences between modes was larger in dry than in wet scenarios. Our results show that models may strongly overestimate root water uptake in dry soils. Benchmark C1.2 was designed to challenge the coupled root-soil models. As the soil only and root only modules were solved well by the different simulators, the differences in the results of C1.2 can be mainly attributed to the coupling approaches and how well the soil water potential gradients that develop around the roots are recognised. In Koch (2022), using a coupling approach with steady-state rhizosphere gradient
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Numerical solutions of the root water pressure head distributions over soil depth in the branched root structure (M3.2) compared to the analytical solution ($k_r$: root radial conductivity, $k_x$: root axial conductance).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Simulated soil water pressure head gradients and times of stress-onset of the single root scenario C1.1. The stress onset time in the sand scenario was zero for the reference solution and all the simulators. (Ref: Reference, CPB: CPlantBox-DuMu<sup>x</sup>, DC: DuMu<sup>x</sup>_CYL, OSR: OpenSimRoot, R: R-SWMS.)</figcaption></figure>

in each time step, the modelling error of implicitly representing the root-soil interface instead of explicitly resolving the interface in the 3D mesh was determined to affect the predicted transpiration rate by 5% or less depending on the modelled scenario. Previous work has shown that the root water uptake from dry soil

was overestimated when the rhizosphere hydraulic conductivity drop was not considered, and that the extent of the overestimation was dependent on the soil grid size (Khare et al., 2022; Koch et al., 2022; Koch, 2020; Schröder et al., 2008). If the grid size is too large, the rhizosphere soil hydraulic conductivity gradients
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Numerical solutions of root water uptake by a static root system architecture compared to the solution according to the DuMu<sup>x</sup> explicit interface implementation. The black solid line depicts the potential transpiration.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Depth profiles of root sink term for water uptake at noon and midnight at three time points.</figcaption></figure>
that may develop in dry soil are underestimated, with the soil water pressure head at the root-soil interface and root water uptake thus being overestimated. We have to keep in mind that the reference solution of C1.2 is a numerical solution itself, albeit a more complex one, as it explicitly resolves the 3D geometry of the root and soil. The reference is thus dependent on the accuracy of its numerical scheme and on the implementation. However, the simulators with a rhizosphere version were in good agreement with the reference solution, while other simulators widely overestimated water availability. We thus concluded that it is necessary to account for the drop in soil hydraulic conductivity in dry soil (in absence of any other processes that modify the rhizosphere) in this scenario.

All other benchmark problems are basic tests that should be solved reasonably well before performing benchmark C1.2. In this study, we learned that all simulators are able to solve benchmarks M2 and M3. For M3, this is to be expected, given the small segment lengths of the root system architecture (Meunier et al., 2017a). We can conclude from this exercise that, given the sufficiently small spatial resolution of the root architecture, the simulators did not have other bugs or mistakes that could have rendered their solutions inaccurate. We also learned that certain parameterisations, such as the clay in M2 or C1.1, are more challenging than others. These new insights led to improvements - and may in future lead to further improvements - of the particular modules of the simulators and help to interpret the results of benchmark C1.2 in a more informed way.

Fig. 7 shows that scenarios (a) and (b) exhibit opposite behaviour for the distributions of the root water pressure heads of the lateral roots. For constant root hydraulic properties, the root water pressure head along the lateral roots is always close to the value of the primary root at the branch point. In contrast, the distribution of root water pressure heads of the lateral roots in the scenario for age-dependent root hydraulic properties range between the value of the primary root at the branch point and the local soil water pressure head at the lateral root tip. This is due to the fact that in the age-dependent root hydraulic properties, root radial conductivities near the tips are large and decline towards the root base (i.e., the older parts of the root), while axial conductances are smaller at the root tip and increase towards the root base. This behaviour was well reproduced by all the simulators. Most of the numerical results of the simulators match the reference solution remarkably well.

Although it is an approximation to the true solution, the steady-rate solution (Schröder et al., 2008) offers a suitable way of improving the root water uptake simulations for our benchmark scenario C1.1, i.e., situations in which soil water potential gradients develop around the root. For all coupled benchmark scenarios, we created an overview of the different ways of coupling the soil and the root domains as well as the different methods used to account for rhizosphere resistance to water flow.

In general, the implementation of functional-structural root architecture models in simulation frameworks requires numerous technical and numerical decisions. To provide an overview of these choices, we created Figure 3 and Table 2, with the details being described in the main text. A more detailed justification of the specific choices made for each simulator is beyond the scope

of this overview, but can be found in the documentation for each simulator.

This benchmarking activity was a challenge from the point of view of communication and organisation. During our meetings, we had to address semantic aspects, technical aspects (e.g., the degree of freedom for each simulator to implement the benchmarks), looking at inconsistencies, proposing or suppressing scenarios, etc. For the final set of benchmark problems, none of the modellers needed to make any substantial changes to the code itself to make the implementation of the benchmarks possible. However, issues regarding spatial grids, convergence criteria, time stepping and computation time had to be addressed (see appendix A for details). The benchmark problems are now part of the continuous development of the different simulators.

By bringing together different research groups and encouraging them to share results obtained with their current model implementations, we created a platform to compare and exchange different modelling solutions. In particular, the willing participation in intentionally challenging problems such as benchmark C1.2 allowed the modellers to evaluate and refine their current modelling approaches. Overall, the benchmarking exercise offers an exclusive study into the current state of the art in functional-structural root architecture modelling and how it is applied over a wide range of modelling groups with different focuses of research. Furthermore, FSRM with a mechanistic representation of dynamic soil processes is a relatively new kind of modelling, at least compared to crop models. We hope that this benchmark will help to stimulate new modelling initiatives and ideas in this field. In future, we plan to expand the benchmark framework to enable further cooperation and exchange between different simulators.

The reference and different numerical solutions are stored on the GitHub repository (https://github.com/RSA-benchmarks/collaborative-comparison), and the automated analysis is available in Jupyter Notebooks. Additional model results may be added to the GitHub repository at any point in the future and will be automatically included in the comparison.

Future efforts may aim to extend the benchmarks from water flow to further processes, such as solute transport or rhizodeposition. A further topic for extended benchmarking would be to investigate growing root systems, where small errors on root uptake evaluation may (or may not) result in an error growing with time. Testing bigger and more complex root systems is another challenge for which a reference solution might not be available. However, an ensemble comparison of different simulators or the use of reference data may be used for such an evaluation.

## SUPPORTING INFORMATION

Appendix A provides the detailed description of the numerical solution methods applied by the different simulators. Appendix B provides a derivation of the hybrid analytical solution of water flow inside a branched root system using exponential functions. Supplementary figure S1 shows a flow chart of the SRI simulator. Supplementary figure S3 illustrates the dependence of M2.2 results on implementation of the top boundary condition as well as soil grid resolution.
# ACKNOWLEDGMENTS

# SOURCES OF FUNDING

This work has partially been funded by the German Research Foundation under Germany's Excellence Strategy, EXC-2070 – 390732324 – PhenoRob and by the Helmholtz Association, Germany—POF4-899. AS was partially funded by the German Federal Ministry of Education and Research (BMBF) in the framework of the funding initiative “Rhizo4Bio”, project “RhizoWheat” (grant 031BO910B). CD benefited from a funding of French ANR for OutLAbMRI project (grant No ANR-19-CE04-0006). GL was partially funded by the German Federal Ministry of Education and Research (BMBF) in the framework of the funding initiative “Rhizo4Bio”, project “CROP” (grant 031BO909A). ML was funded through the priority program 2089 “Rhizosphere spatiotemporal organization – a key to rhizosphere functions” funded by the German Research Foundation DFG under the project number 403641034. TK has been financially supported by the European Union's Horizon 2020 Research and Innovation programme under the Marie Skłodowska-Curie Actions Grant agreement No 801133 and by the German Research Foundation, within the Collaborative Research Center on Interface-Driven Multi-Field Processes in Porous Media (SFB 1313, grant 327154368). TS was funded by the BMBF in the framework of the funding initiative ‘Soil as a Sustainable Resource for the Bioeconomy BonaRes’, project ‘BonaRes (Module A): Sustainable Subsoil Management – Soil3; subproject 3’ (grant 031B1066C). VC is a Research Associate of the Fund for Scientific Research (FNRS), co-funded by the European Union (ERC grant 101043083).

# CONFLICT OF INTERESTS

None declared.

# CONTRIBUTIONS BY THE AUTHORS

AS initiated this benchmark initiative. ML provided the RSML files for the M3 and C1.2 benchmarks. TK implemented and ran the simulations for M2, M3 (DuMu<sup>x</sup>), C1.1 (DuMu<sup>x</sup><sub>cyl</sub>), C1.2 (DuMu<sup>x</sup>-ks) and the C1.2 reference solution (Dumux_explicit_interface). DL implemented and ran the simulations for M3 (CPlantBox), C1.1, and C1.2 (CPlantBox-DuMux), with support from DK, TK, and AS. AK, AH and MJ implemented and ran all benchmark problems with R-SWMS. CK, ES and JP implemented and ran all benchmark problems with OpenSimRoot. CD implemented and ran all benchmark problems with SRI. ML, TK, GL and JV contributed to the creation of overview figures. All authors discussed the results and contributed to the final manuscript.

# DATA AVAILABILITY

All simulation results and reference solutions are available through the GitHub repository https://github.com/RSA-bench

marks/collaborative-comparison. The data should be cited using the relevant doi on Zenodo, doi: 10.5281/zenodo.7844443.

The codes used to implement the benchmark problems are available from the individual developers upon request and/or via the following repositories:

## DuMu<sup>x</sup>

M2: https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/tree/master/test/porousmediumflow/richards/benchmarks

M3: https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/tree/master/test/porousmediumflow/1p/rootbenchmark

C1.1 (rotationally-symmetric): https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/tree/master/test/porousmediumflow/richards/annulus

C1.2 (kernel): https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/tree/master/test/multidomain/embedded/1d3d/root_soil_benchmark

C1.1 (1D-3D) and C1.2 (explicit interface): https://git.iws.uni-stuttgart.de/timok/dumux-rootsoilbenchmarking

## CPlantBox-DuMu<sup>x</sup>

M2: https://github.com/Plant-Root-Soil-Interactions-Modelling/dumux-rosi/tree/master/python/soil

M3: https://github.com/Plant-Root-Soil-Interactions-Modelling/dumux-rosi/tree/master/python/roots

C1.1 and C1.2: https://github.com/Plant-Root-Soil-Interactions-Modelling/dumux-rosi/python/coupled/

C1.2<sub>rhizo</sub>: https://github.com/Plant-Root-Soil-Interactions-Modelling/dumux-rosi/python/coupled_rhizo/

## OpenSimRoot

https://gitlab.com/rootmodels/OpenSimRoot/-/tree/OpenSimRoot_v2_benchrun/OpenSimRoot/tests/benchmarks2020

# SRI

The code is available from the author upon request (claude.doussan@inrae.fr).

# A. DETAILED DESCRIPTION OF THE NUMERICAL SOLUTION METHODS APPLIED BY THE DIFFERENT SIMULATORS

In this section, we describe the numerical approaches used by the different simulators in each of the benchmarks in detail.

## A.1. DuMu<sup>x</sup>

DuMu<sup>x</sup> is a generic modular C++ framework specializing in flow and transport in porous media (Flemisch et al., 2011; Koch et al., 2021) and is based on the Distributed Unified Numerics Environement (DUNE) (Bastian et al., 2021). DuMu<sup>x</sup> implements various finite volume and control-volume finite element discretisation methods in space, implicit and
explicit Euler discretisation in time (see below for the discretisation schemes employed for each benchmark), and a generic Newton method based on a numerical approximation of the Jacobian matrix of the possibly nonlinear discretised porous flow equation. Moreover, it provides a framework to couple PDEs posed on different domains or using different discretisation methods (Koch et al., 2021) and to solve coupled problems monolithically, i.e., in an approach where all equations are solved simultaneously. An overview of implemented mathematical models and processes in the context of root-soil interaction is given in Koch et al. (2018). The generic DUNE grid interface (Bastian et al., 2008) supports both structured and unstructured locally adaptive grids (Alkämper et al., 2016) as well as embedded network grids (Sander et al., 2017).

**Implementation of benchmark problems M2** The Richards equation was solved using a cell-centered finite volume scheme with two-point flux approximation (TPFA) as spatial discretisation (Koch et al., 2021) using the arithmetic mean of the hydraulic conductivity in the approximation of the numerical flux. In time, we discretised with an implicit Euler scheme. The nonlinear discrete equation is solved with a quasi-Newton method where the Jacobian is approximated by numerical differentiation of the residual.

In benchmarks M2.1 and M2.2, we used a structured Cartesian grid with a vertical resolution of 1,000 cells and only 1 cell in the horizontal directions (quasi-1D). The grid was locally refined towards the surface. Automated time stepping was based on the number of Newton iterations. To mimic the analytical solution, gravity was considered in benchmark M2.1 and disregarded in benchmark M2.2. The switch from Neumann to Dirichlet boundary conditions at the top of the soil domain in the evaporation benchmark scenario M2.2 was implemented by only prescribing Neumann boundary conditions, but taking the minimum of the potential and critical evaporation rates (where the critical evaporation rate is the rate needed to maintain the soil water pressure head at the top at the critical soil water pressure head).

**Implementation of benchmark problems M3** Darcy's equation was solved on the network with a vertex-centered finite volume scheme (also referred to as Box method or control-volume finite element scheme with linear Lagrangian basis functions (Huber and Helmhig, 2000)). We used the root mesh provided and refined each element once. The numerical solution was verified numerically to converge against the analytic solution with a second order grid refinement in $L^2$-norm, as expected.

**Implementation of benchmark problem C1.1** Benchmark C1.1 was implemented as a coupled problem (using a cuboid soil domain with equal volume discretisation by an unstructured grid that is locally refined around the root but does not conform to the root-soil interface) with negligible resistance in the root and Neumann boundary condition on the root collar. As monolithic coupling method, we used a generic (discretisation-agnostic) method to couple tubular network domains to their embedding domain described in Köppl et al. (2018) (DuMu<sup>x</sup>-CYL). The sink term is distributed on the implicit root surface

(surface sink term), and the soil water pressure head in the coupling term is evaluated as the perimeter-average of the soil water pressure head. We discretized Darcy's equation in the root with a cell-centered finite volume method (TPFA) and the Richards equation in the soil with a vertex-centered finite volume method (Box). Adaptive time stepping based on the number of Newton iterations was used. The time steps at the end of the simulation were chosen in such a way that the critical root collar pressure head is reached within a relative error of $10^{-3}$. In each time step, the coupled system is solved using Newton's method. The Jacobian of the monolithic system contains entries for the degrees of freedom (df) in both the root and the soil domain, and is approximated by numerical differentiation of the residual.

**Implementation of benchmark problem C1.2** In addition to the reference solution described above, DuMu<sup>x</sup> contributed to benchmark C1.2 using the “implicit interface with kernel support” (DuMu<sup>x</sup>-ks) method. This method disregards the roots in the soil domain and instead of a boundary condition, the coupling term appears as a macroscopic sink term in the Richards equation where the soil water pressure head used in the coupling is conceptually given by the perimeter-average of the soil water pressure head on the root-soil surface (Koch et al., 2022). Moreover, the sink term is distributed around each root segment in a fixed local cylinder whose radius is three times the radius of the root (independent of the mesh size). The distribution of the source in a fixed region regularizes the soil water pressure head solution locally in the vicinity of the root. The perimeter-average of the soil water pressure head is thus a bad approximation of the pressure head at the soil-root interface. Therefore, it is instead reconstructed based on the soil water pressure head evaluated at the root centerline and the local analytic solution of the steady Richards equation on a disc (Koch et al., 2022). The solution obtained using the “implicit interface with kernel support” method is an approximation of the solution obtained using the “explicit interface” method. The model error was numerically estimated by Koch et al. (2022) (about 3% differences in local source terms in challenging (dry soil) problems in steady state). We discretized Darcy's equation in the root domain and the Richards equation in the soil domain with a cell-centered finite volume method (TPFA). The structured Cartesian regular soil grid had 7680 cells (diameter 0.86 cm). The runtime for the cases in C1.2 was about 2 minutes on a standard laptop.

Most of the benchmark cases are now also part of the continuous integration pipelines of DuMu<sup>x</sup> (see code availability).

## A.2. CPlantBox-DuMu<sup>x</sup>

In this benchmark, CPlantBox (Zhou et al., 2020; Schnepf et al., 2018) is used as functional-structural modelling framework. The soil subproblems are solved using DuMu<sup>x</sup>, while the root subproblems and the coupling are solved with CPlantBox.

CPlantBox, written in C++ and Python, is a functional-structural plant architecture model that creates dynamic 3D plant geometries (Zhou et al., 2020; Schnepf et al., 2018). In this benchmark, we focus exclusively on static root architectures. CPlantBox supports various grid or parameter input formats and offers auxiliary functions to support functional-structural plant modelling and coupling to external solvers. CPlantBox can
solve water flow on root grids using either the Doussan model (Doussan et al., 2006), i.e. a Darcy-type axial flow model in which the axial flow is proportional to the root water pressure head gradient between two nodes, and with a root water uptake source term that is proportional to the pressure head difference between the xylem and the root surface. Alternatively, the hybrid analytical-numerical solution described in appendix B by Meunier et al. (2017b) can be used. At known soil water pressure head, the solution is independent from the spatial discretisation (root segment length). Furthermore, CPlantBox offers an interface to couple the root architecture with soil grids and enables the use of external solvers for the soil model via Python bindings. In this application, DuMu<sup>x</sup> was used to solve the soil subproblems.

*Implementation of benchmark problems M2* The soil problems were solved with DuMu<sup>x</sup> (see M2 description of DuMu<sup>x</sup>).

*Implementation of benchmark problems M3* The root-only problems were solved using the CPlantBox implementation of the analytical solution for a single root (M3.1) and the hybrid analytical solution for a root system (M3.2), respectively. The CPlantBox solution is thus equal to the reference solution.

We used the root architecture provided as the numerical grid that determined the nodes and edges of the tree graph representation of the root system (i.e. we did not restructure the numerical grid). The flow along each edge is computed by the analytical solution for a single root. As described in Meunier et al. (2017b), the graph Laplacian is then used to compute the water volume flowing into or out of each node.

*Implementation of benchmark problem C1.1* Although the reference solution only considers the rhizosphere from the root surface to the bulk soil and is, therefore, essentially a soil-only problem with the root defining one boundary condition, we implemented benchmark C1.1 as a coupled problem of a soil cylinder with a height of 1 cm and a radius of 0.6 cm, and a single vertical root in the middle. Therefore, we first determined the root collar boundary condition as a volumetric flow rate, $Q_{x, collar}$, from the given flux at the root boundary, $q_{root}$, as $Q_{x, collar} = 2r_{root}\pi q_{root}l_{root}$, where $r_{root}$ and $l_{root}$ are the root radius and length, respectively. To simulate the horizontal water flow towards the root, the cylindrical soil domain was discretised as follows: Using PyGmsh (Schlömer, 2022), a triangular mesh of a disk with a mesh size of 0.5 mm was extruded in the z-direction with one layer to obtain a pentahedral mesh of a quasi-2D soil domain.

The soil subproblem was then solved with DuMu<sup>x</sup>, using a cell-centered finite volume scheme with TPFA in space and an implicit Euler scheme in time. The root subproblem was solved using the analytical solution for a single root, automatically switching from Neumann to Dirichlet boundary condition when the xylem water pressure head at the root collar reached the wilting point. The root water uptake sink term for each soil control element is computed as the sum of the water uptake by each root segment whose centerline is located within it. The root water uptake by each root segment is computed based on the water pressure head difference between the xylem and the respective

soil control element. As we use a finite volume scheme, there is only one value per soil control element.

For the sequential coupling, we used explicit time stepping with a coupling time step of 2.4 h. In each coupling time step, xylem water pressure heads are computed based on the current soil water pressure heads and the root collar boundary condition. The root water uptake for each root segment is then computed, obtaining the sink terms for the soil control elements. They are defined for the subsequent solution of the soil subproblem using DuMu<sup>x</sup>. The solution of the soil subproblem itself uses automated time stepping, as described above. The next time step starts with these updated soil water pressure heads.

At the onset of stress (i.e. -15,000 cm at the root collar), we sampled the soil water pressure heads along a radial line between the center axis and the outer edge of the soil cylinder at 40 sampling points.

*Implementation of benchmark problem C1.2* To solve benchmark C1.2, we used a regular cubic grid with an edge length of 1 cm to discretise the soil domain and provided an RSML-file for the discretisation of the root architecture as tree graph. As in the previous benchmarks, the soil subproblem was solved by DuMu<sup>x</sup>, and the root subproblem was solved according to the CPlantBox implementation of the hybrid analytical-numerical solution of Meunier et al. (2022).

For the sequential coupling, we used a time step of 360 s and the same coupling scheme described in benchmark C1.1. We implemented two setups: In the first setup, we set the soil water pressure head at the soil-root interface to be equal to the soil water pressure head of the soil control element in which the axis of the corresponding root segment was located. In the second setup, we used a sink term definition that explicitly considers subresolution drops in soil water pressure head in the rhizosphere under dry soil conditions. This approach is described in detail by Mai et al. (2019). In brief, we solve a 1D radially symmetric rhizosphere model for each root segment in a mass conservative way: At the start of the simulation, we define hollow soil cylinders around each root segment in such a way that the sum of volumes of the hollow soil cylinders equals the volume of the soil control element, in this case 1 cm<sup>3</sup>. The inner radius of each hollow cylinder is equal to the root radius. The outer radius is computed by dividing the volume of the soil control element between all the root segments within it, proportional to their root volumes. The initial conditions are set to be homogeneous and equal to the value of the soil control element. For each 1D radial model, we used a logarithmic discretisation of 10 grid points, with more grid points near the root surface. In each time step, all the 1D radial models are solved with DuMu<sup>x</sup>, using a model capable of radial symmetry and the results of the previous time step as initial condition. The rhizosphere models thus link the root subproblem with the (macroscopic) soil subproblem. For each time step, the net influx (positive or negative) into each soil control element is divided among all rhizosphere models proportionally to their root segment volumes and prescribed as the outer boundary conditions to the rhizosphere models.

When using the rhizosphere models, the sequential coupling scheme is adapted in such a way that we first solve the root subproblem, which provides us with the xylem pressure heads. We
then solve the local rhizosphere models, which provides us with the root water uptake values for each root segment, and finally, we solve the soil subproblem with the given root water uptake sink terms.

### A.3. OPENSIMROOT

OPENSIMROOT is written in C++ and is a functional-structural plant model used to describe the functionality and growth of root systems (Postma et al., 2017). The models and interactions used in a simulation are specified in XML-formatted input files (Schäfer et al., 2022). The code was written to be modular, meaning that little knowledge - beyond that of the API - is needed to implement new models. Plant growth is determined by growth rates specified in the input file and can be altered based on plant nutrient and carbon status as well as local soil conditions.

To implement coupled systems of partial, ordinary, and algebraic equations, OPENSIMROOT uses a loosely coupled approach. The general coupling scheme for simulating over time uses the predictor-corrector method, Runge-Kutta 4 by default. The individual ordinary differential equations (ODEs) and partial differential equations (PDEs) are each solved with individual time steps and the exchange of information using linear interpolation. To avoid interpolation errors, the time steps are frequently synchronized.

The benchmarks problems are fully defined in the input files, which can be downloaded at https://jugit.fz-juelich.de/rootmodels/OpenSimRoot/-/tree/OpenSimRoot_v2_benchrun/OpenSimRoot/tests/benchmarks2020. We added them to our standard code testing and wrote scripts (R/Python/Shell) for each benchmark, converting the output of OPENSIMROOT to the format required by the benchmarks.

**Implementation of benchmark problems M2** The Richards equation solver is a C++ implementation of the SWMS3D code (Simunek et al., 1995), but with many simplifications. It still uses the same finite element scheme, which is implicit in time, using a fixed point iteration until convergence is reached (Huyakorn and Pinder, 1983). The time step is reduced (and repeated) if the system fails to converge. From the benchmarks, we learned that allowing more iterations is a faster strategy than aggressively reducing the time step. The equation system is solved with a PCG (preconditioned conjugate gradient) solver, and the benchmarks required a relatively high level of precision ($1 \times 10^{-15}$). In rare cases, the convergence of the PCG solver was slow, and we programmed the code to restart the time step with a smaller delta.

To run the benchmarks, we had to make some changes to the code. These changes were related to (1) the water retention curve implementing the standard van Genuchten curve, (2) increasing the precision of the PCG solver, (3) increasing the maximum permitted iterations of the time solver of the Richards equation, and (4) less aggressive scaling of the time step of the solver of the Richards equation. The corresponding code is found in the “benchmark branch” of the git repository https://gitlab.com/rootmodels/OpenSimRoot/-/tree/73ea68dd45e3c1fef5f43e0f38b5472807069a76. This branch has a macro that can compile with or without the

zero-gravity assumption, depending on the requirements of the different benchmarks.

M2.2 was a difficult benchmark to solve, as the switch from Neumann to Dirichlet causes a discontinuity. To represent the switch from Neumann to Dirichlet boundary conditions in the evaporation benchmark scenario M2.2, OPENSIMROOT defines a smooth function at the top of the soil domain, which is intended to reduce the potential evaporation ahead of reaching the condition for the Dirichlet switch. A transition period thus smoothes the problem between potential evaporation and the switching to Dirichlet, which speeds up the simulation considerably and makes it less sensitive to spatial resolution.

**Implementation of benchmark problems M3** The root system is represented by nodes with associated root lengths. Since OPENSIMROOT has no notion of non-growing static root systems, the root systems were implemented to grow (quickly) to size and then stop growing. The benchmarks have corresponding time offsets. The OPENSIMROOT input files contain tables with the time-dependent location of each root tip, ensuring that the root system grows as defined by the benchmarks. The root system is represented by a hydraulic tree, and the corresponding system of linear equations is solved with a PCG solver.

**Implementation of benchmark problem C1.1** In benchmark C1.1, we aggregated the nodal values in 3D to 1D radial distances by linear interpolation.

**Implementation of benchmark problem C1.2** The root system is represented by nodes with associated root length. During the growth of the root system, the creation of the nodes is broadcasted and registered by the Richards equation solver, which maps their location to the nearby nodes in the rectangular finite element grid. Using an inverse-distance-weighted average, the values in the grid and the values on the nodes can be (bidirectionally) mapped to set both the sink term and the hydraulic head at the root surface. Both are updated during every iteration of the time solver of the Richards equation solver. Two parameters can be set: the maximum distance for the neighborhood search and the exponent of the (inverted) distance used to determine the weight of each nearby node. Spreading the sink terms over several nodes thus stabilises the solution. OPENSIMROOT treats individual differential equations as “mini models” (independent class), which communicate over a common, request-driven, API. For example, the Richards equation solver requests the water uptake of all root segments, and these, in return, request this to be computed by the module solving the hydraulic root architecture model. For any given soil water pressure heads at the root surfaces, this model is a steady-state model as described by Doussan et al. (1998). The root hydraulic architecture model requests the soil water pressure heads from the Richards equation solver, which resolves them based on the current (predicted) root water pressure heads. The original request of the Richards equation is thus resolved, and the process is repeated for the next fixed-point-iteration of the coupled problem.

Benchmark C1.2 creates a strong difference between the radial conductivity of the roots and the rhizosphere when the
soil is drying, making it grid-resolution sensitive and very diffi-
cult to solve. We implemented a “rhizo” version, which assumes
that the radial hydraulic resistance includes both the root radial
resistance (i.e., the inverse of the root radial conductivity k<sub>r</sub>)
and a rhizosphere resistance. We define the radial resistance as:
1/K<sub>total</sub> = 1/((S<sub>root</sub>k<sub>r</sub>) + w/((K<sub>soil</sub>S<sub>root</sub>))), where S<sub>root</sub> (cm<sup>2</sup>) is the
root segment surface area, and w (cm) is an arbitrary weight that
worked well with a value of 4.

A.4. R-SWMS

R-SWMS (Javaux et al., 2008) is a simulator written mainly
in Fortran, which solves the 3D water flow and solute transport
equations in the soil, in the plant root system (considered as a
1D connected network, i.e., a hydraulic tree), and the exchanges
between them via a source/sink term approach. In other appli-
cations, root development and maturation may be simulated,
potentially in response to soil mechanics, heat, and the distribu-
tion of solutes in the soil (Landl et al., 2017; Somma et al., 1998;
Clausnitzer and Hopmans, 1994).

The numerical solution for the water flow equation in R-SWMS is almost identical to the SWMS_3D code proposed by Simunek et al. (1995). The Richards equation is solved using the Galerkin finite element method. The discretisation of the flow domain is realized by dividing it into cubical elements and subdividing each cube into 5 tetrahedral subelements. This is different than in the original code of Somma et al. (1998), where half cubes were subdivided into 3 elements each. The corners of the subelements are the domain nodes. This procedure leads to a system of ordinary differential equation. The original code was also adapted to solve the Richards equation based on the mixed form of the Richards equation, enabling an accurate mass balance (Celia et al., 1990) with an implicit method for time solution. The resulting nonlinear system of equations is solved by Picard iteration.

The soil module can either be coupled with a static (i.e., no root growth) or with a dynamically growing root system. Multiple independent growing plants can be simulated. The steady-state water flow equations in the plant-root system described by Doussan et al. (1998) are solved with a linear solver. The coupling between the soil and the root water fluxes within a time step is realized by sequentially solving both systems of equations until convergence is reached for the water pressure head in roots and soils and for fluxes. The sink term for the Richards equation is calculated based on the volumetric average of the uptake fluxes simulated by the Doussan equations. The water pressure at the soil-root surface of the segment used in the Doussan equation was derived from the water potentials at the soil nodes using the inverse distance-weighted average of the pressure heads of the surrounding soil grid nodes. If the segment crosses multiple voxels, the segment is divided into subsegments and a specific water pressure at the soil-root surface is calculated per subsegment and then averaged.

*Implementation of benchmark problems M2* In benchmarks M2.1 and M2.2, we used a structured Cartesian grid with a constant vertical resolution of 0.2cm resulting in 1,000 cells, and with only 1 cell (1 cm x 1 cm) in the horizontal directions

(quasi-1D). In the event of non-convergence, the tolerance on
the soil pressure head was increased (up to 0.1 cm).

*Implementation of benchmark problems M3* Water flow in the root system was solved for one time step in a uniform soil with a soil conductivity of 0 cm/d and a soil hydraulic capacity of 0 (van Genuchten parameter $\alpha$ set to infinity). The root nodes of the RSML file were transformed into the input file format of R-SWMS via a Matlab routine.

*Implementation of benchmark problems C1.1 and C1.2* To calculate the soil resistance at the infra-voxel scale due to the soil conductivity drop around roots, an additional module was used in R-SWMS (R-SWMS-Rhizo). The local water flow within the voxel is represented by the water potential difference between the xylem and the voxel edge and two hydraulic resistances in series: the radial root resistance and the effective infra-voxel soil resistance. The radial root resistance of each root node of the hydraulic tree is then replaced by the sum of the two hydraulic resistances in series. The infra-voxel soil hydraulic resistance is calculated at each iteration from the steady-rate axisymmetrical solution of the Richards equation towards the soil-root interface from Schröder et al. (2008). When multiple roots were present, the superposition approach of Beudez et al. (2013) was used. The sequential coupling between the soil and the root systems within a time step is realized with iterations between soil finite element solution for soil water potential and the uptake flux rates, resulting in the simultaneous solution of soil grid water potential, root xylem potential and root water uptake rates. To avoid instability in the code, the maximum number of iterations was set at 5 to calculate the radial resistance. This infra-voxel solution can be omitted if it is not required.

A.5. SRI

SRI is a modular framework for investigating soil water flow and root uptake. It is based on C++ object-oriented codes that can be run independently or coupled. The 2D/3D soil water flow model is based on the generic finite element framework FAFEMO (Mesgouez and Lefeuve-Mesgouez, 2009) and adapted to solve the Richards equation based on the mixed form of the Richards equation, enabling accurate mass balance (Celia et al., 1990) with an implicit time solution method. Regular or unstructured tetrahedral meshes can be used. The root architecture models that can be currently used are Pmais (Pagès et al., 1989), RootTyp (Pagès et al., 2004), ArchiSimple (Pagès et al., 2014), or an RSML file as root input. The root system function for water uptake and transport within the root system is modelled following Doussan et al. (1998). The soil model can either be coupled with a static (i.e., no root growth) root model, or with a dynamically growing root model. The general algorithm for the solution is given in the flowchart S2.

*Implementation of benchmark problems M2* The FAFEMO-Richards finite element solver was used without adding the root components. For M2.1, the soil was meshed with a uniform cell grid size of 1 cm, formed with tetrahedral elements. For the evaporation case of M2.2, we needed to decrease the z mesh size
down to 0.5 mm at the soil surface, with a geometric increase with depth to 8 mm (mesh size is 4 mm in x, y direction). The 1 cm regular meshing was found to result in large errors compared to the reference solution. In both M2.1 and M2.2, the maximum time step was restricted to 100 s.

*Implementation of benchmark problems M3* The root nodes of the RSML file were used in the “hydraulic tree model” of root architecture, without adding any additional nodes for computation. A code was developed to translate RSML nodes encoded in the description of the root architecture and graph for root flow calculations according to Doussan et al. (1998).

*Implementation of benchmark problem C1.1* Did not participate.

*Implementation of benchmark problem C1.2* In order to solve the fine infra grid size local gradient of water potential near roots, within a larger mesh element size (i.e. about 1 cm), the analytical 3D spherical solution of steady-rate unsaturated flow is used for each segment within the element with an equivalent inner sphere of the same surface area as the cylindrical root segment, and an external radius of the same size as the mesh element. The mean soil mesh voxel water potential is imposed on the external radius while the root uptake flux rate is imposed on the internal radius. The influence of other root segments located within the external radius (i.e., within the mesh cell or in neighbouring cells) on the water potential drop can be accounted for using a superposition approach described (in 2D) by Beudez et al. (2013). The water potential at the root segment surface can be determined from the non-linear system of equations (root flow-Kirchhoff transform of soil water analytical solution) using a Picard iterative method, an under-relaxation and vector acceleration algorithm derived from the Aitken method (Ramière and Helfer, 2015). With this solution, the total root uptake is introduced as a distributed root sink in the finite element solution. A global soil-plant solution within the time step is found from iterations between the soil finite element solution for soil water potential and the uptake flux rates (see flowchart in Fig. S2) resulting in the simultaneous solution of soil grid water potential, root xylem potential and root water uptake rates. The infra-voxel rhizosphere solution can be turned off, and the solution is only obtained from the convergence of uptake incorporated in the distributed root sink term. Here, the water potential at the root surface, which is used to calculate the water uptake of root segments in a grid element, is simply the average of the soil water potential of nodes of the grid element. Similar to M3, the root RSML file was used without adding additional nodes for root flow calculations. Soil mesh elements were 1 cm in size with or without infra-voxel rhizosphere calculations. The maximum number of fixed point iterations used for the infra-voxel rhizosphere solution had to be increased (from about 100 to more than 5,000) to ensure convergence combined with small time steps, although this was at the expense of calculation time. This was more evident for C12b where the simulation over the entire 3 days was not possible with a sufficiently reasonable computation time. Therefore, SRI results are only available for benchmark C12a.

## B. DERIVATION OF THE HYBRID ANALYTICAL SOLUTION OF WATER FLOW INSIDE A BRANCHED ROOT SYSTEM

We describe the water flow within the root xylem as Darcy-type flow, and with a root water uptake source term that is proportional to the pressure difference between xylem and root surface (Roose and Fowler, 2004; Doussan et al., 2006).

The axial volumetric flow rate $q_x$ [cm$^3$ day$^{-1}$] in an individual root segment is given by

$$q_x = -k_x \left( \frac{\partial \psi_x}{\partial l} + e_z \cdot v \right). \quad (B1)$$

The parameter $k_x$ is the axial conductance [cm$^3$ day$^{-1}$], $\psi_x$ is the pressure inside the xylem [cm], $e_z$ the unit vector in z-direction [1], and $v$ the normed direction of the xylem [1]. Thus above equation can be expressed as

$$q_x = -k_x \left( \frac{\partial \psi_x}{\partial l} + v_3 \right), \quad (B2)$$

where $v_3$ is the z-component of the normed xylem direction.

The radial volumetric flow rate $q_r$ [cm$^3$ day$^{-1}$] into the root is given by

$$q_r = 2a\pi k_r d_l (\psi_s - \psi_x), \quad (B3)$$

where $a$ is the root radius [cm], $d_l$ is an infinitesimal length [cm], $k_r$ is the radial conductivity [day$^{-1}$], and $\psi_s$ is the soil water pressure head at the soil root interface [cm].

For a constant $k_x$ along each segment mass conservation yields

$$2a\pi k_r (\psi_s - \psi_x) = -k_x \frac{\partial^2 \psi_x}{\partial l^2}. \quad (B4)$$

### B.1. Analytical solution for a single root segment

For constant $\psi_s$, $k_r > 0$ and $k_x > 0$ we can solve above equation (Eqn. B4) yielding the analytical solution

$$\psi_x(l) := \psi_s + d_i e^{\tau l} + d_j e^{-\tau l} \quad (B5)$$

with

$$\tau := \sqrt{2a\pi k_r / k_x} \quad (B6)$$

[cm$^{-1}$]. The constants $d_i$ and $d_j$ [cm] can be calculated from the boundary conditions (equals Meunier et al. (2017b), B1).

Note that for $k_r = 0$ or $a = 0$ the solution of Eqn. B4 simplifies to

$$\psi_x(l) := d_i z + d_j. \quad (B7)$$

First we apply Dirichlet boundary conditions at the beginning ($l = 0$), and the end ($l = L$) of a single segment, i.e.

$$\psi_x(0) = \psi_{x,i} \quad (B8)$$

$$\psi_x(L) = \psi_{x,j}, \quad (B9)$$
where $\psi_{x,i}$ and $\psi_{x,j}$ denote the root water pressure heads at node $i$ and $j$ of a segment $s_{ij}$. They are inserted into the analytic solution and yield the two equations

$$ \begin{pmatrix} 1 & 1 \\ e^{\tau L} & e^{-\tau L} \end{pmatrix} \begin{pmatrix} d_i \\ d_j \end{pmatrix} = \begin{pmatrix} \psi_{x,i} - \psi_{s,r} \\ \psi_{x,j} - \psi_{s,r} \end{pmatrix}. \quad (B10) $$

From these equations we can calculate the constants $d_i$, and $d_j$ [cm]:

$$ \begin{pmatrix} d_i \\ d_j \end{pmatrix} = \delta^{-1} \begin{pmatrix} e^{-\tau L} & -1 \\ -e^{\tau L} & 1 \end{pmatrix} \begin{pmatrix} \bar{\psi}_i \\ \bar{\psi}_j \end{pmatrix} \quad (B11) $$

with

$$ \delta := e^{-\tau L} - e^{\tau L}, \quad (B12) $$

and

$$ \bar{\psi}_i := \psi_{x,i} - \psi_{s,ij} \quad (B13) $$

$$ \bar{\psi}_j := \psi_{x,j} - \psi_{s,ij}, \quad (B14) $$

where $\psi_{s,ij}$ is the soil water pressure head around the segment connecting node $i$ and $j$ is constant. For the simpler case $k_r = 0$ or $a = 0$ we obtain

$$ \begin{pmatrix} d_i \\ d_j \end{pmatrix} = \begin{pmatrix} (\psi_{x,j} - \psi_{x,i})/L \\ \psi_{x,i} \end{pmatrix}. \quad (B15) $$

## B.2. Analytical solution for connected segments

Multiple segments are connected to represent the root system. We first use the solution for a single segment with node indices $i$ and $j$ for the case of two Dirichlet boundary conditions (see Eqn. B10) and with the constants $d_i$ and $d_j$ (see Eqn. B11) we can write down the axial volumetric flow rate, as

$$ q_x = -k_x \left( \frac{\partial \psi_x}{\partial z} + v_3 \right) \quad (B16) $$

$$ q_{x,ij}(l) = -k_x (d_i \tau e^{\tau l} - d_j \tau e^{-\tau l} + v_3) \quad (B17) $$

(Meunier et al. (2017b), B2). Inserting the constants $d_i$ and $d_j$ (see Eqn B11) yields the explicit axial volumetric flow rate equation

$$ q_{x,ij}(l) = -k_x \left( \frac{\tau}{\delta} \left[ (e^{-\tau L} \bar{\psi}_i - \bar{\psi}_j) e^{\tau l} - (-e^{\tau L} \bar{\psi}_i + \bar{\psi}_j) e^{-\tau l} \right] + v_3 \right) \quad (B18) $$

Evaluation of the axial volumetric flow rate at node $i$, i.e. $l = 0$ yields

$$ q_{x,ij}(0) = -k_x \left( \frac{\tau}{\delta} \left[ (e^{-\tau L} + e^{\tau L}) \bar{\psi}_i - 2 \bar{\psi}_j \right] + v_3 \right). \quad (B19) $$

For conservation of mass in each node $i$ the sum of all axial volumetric flow rates have to cancel out:

$$ \sum_{j \in N(i)} q_{x,ij}(0) = 0, \quad (B20) $$

where $N(i)$ are the indices of the segments $s_{ij}$ starting from node index $i$ to $j$.

The explicit axial volumetric flow rate $q_{i,j}(0)$ can be expressed as linear equation in $\bar{\psi}_i$ and $\bar{\psi}_j$, and for fixed $i$ conservation of mass yields

$$ \sum_{j \in N(i)} c_{ii} \bar{\psi}_i + c_{ij} \bar{\psi}_j - k_x v_3 = 0 \quad (B21) $$

and if we solve for $\psi_{x,i}$ we obtain

$$ \sum_{j \in N(i)} c_{ii} \psi_{x,i} + c_{ij} \psi_{x,j} - b_i = 0 \quad (B22) $$

with

$$ c_{ii} := -k_x \frac{\tau}{\delta} (e^{-\tau L} + e^{\tau L}) \quad (B23) $$

$$ c_{ij} := 2k_x \frac{\tau}{\delta} \quad (B24) $$

$$ b_i := k_x v_3 + c_{ii} \psi_{s,ij} + c_{ij} \psi_{s,ij}, \quad (B25) $$

with units $c_{ii}$, and $c_{ij}$ [cm$^2$/day], and $b_i$ [cm$^3$/day].

The soil water potential only enters on the right hand side. To solve for $\psi_x$ we need to solve

$$ C \psi_x = b. \quad (B26) $$

The results $\psi_x$ depend on the pressure head at the root-soil interface $\psi_s$ and the root collar boundary condition. For a Neumann boundary condition the load vector $b$ is altered accordingly, for Dirichlet boundary condition the row of the root collar is altered in matrix C and vector $b$.

## B.3. Implementation

The model and corresponding solutions are implemented in CPlantBox in Python. For performance reason the assembling of the matrix is additionally implemented in C++. The solving of linear system is done in Python using SciPy.

## LITERATURE CITED

Alkämper, M., Dedner, A., Klöfkorn, R., and Nolte, M. (2016). The dune-alugrid module. *Archive of Numerical Software*, Vol 4.

Bachmair, S., Weiler, M., and Nuetzmann, G. (2010). Benchmarking of two dual-permeability models under different land use and land cover. *Vadose Zone Journal*, 9(2):226–237.

Bastian, P., Blatt, M., Dedner, A., Dreier, N.-A., Engwer, C., Fritze, R., Gräser, C., Grüninger, C., Kempf, D., Klöfkorn, R., Ohlberger, M., and Sander, O. (2021). The dune framework: Basic concepts and recent developments. *Computers & Mathematics with Applications*, 81:75–112. Development and Application of Open-source Software for Problems with Numerical PDEs.

Bastian, P., Blatt, M., Dedner, A., Engwer, C., Klöfkorn, R., Ohlberger, M., and Sander, O. (2008). A generic grid interface for parallel and adaptive scientific computing. part i: abstract framework. *Computing*, 82(2-3):103–119.

Beudez, N., Doussan, C., Lefeuve-Mesgouez, G., and Mesgouez, A. (2013). Influence of three root spatial arrangement on soil water flow and uptake. results from an explicit and an equivalent, upscaled, model. *Procedia Environmental Sciences*, 19:37–46. Four Decades of Progress in Monitoring and Modeling of Processes in the Soil-Plant-Atmosphere System: Applications and Challenges.
Cai, G., Vanderborght, J., Couvreur, V., Mboh, C. M., and Vereecken, H. (2018). Parameterization of root water uptake models considering dynamic root distributions and water uptake compensation. *Vadose Zone Journal*, 17(1):160125.

Celia, M. A., Bouloutas, E. T., and Zarba, R. L. (1990). A general mass-conservative numerical solution for the unsaturated flow equation. *Water Resources Research*, 26(7):1483-1496.

Clausnitzer, V. and Hopmans, J. W. (1994). Simultaneous modeling of transient three-dimensional root growth and soil water flow. *Plant and Soil*, 164(2):299-314.

Custos, J.-M., Moyne, C., and Sterckeman, T. (2020). How root nutrient uptake affects rhizosphere ph: A modelling study. *Geoderma*, 369:114314.

de Jong van Lier, Q., van Dam, J. C., Durigon, A., dos Santos, M. A., and Metselaar, K. (2013). Modeling water potentials and flows in the soil-plant system comparing hydraulic resistances and transpiration reduction functions. *Vadose Zone Journal*, 12(3):vzj2013.02.0039.

de Willigen, P., van Dam, J. C., Javaux, M., and Heinen, M. (2012). Root water uptake as simulated by three soil water flow models. *Vadose Zone Journal*, 11(3):vzj2012.0018.

Doussan, C., Pagès, L., and Vercambre, G. (1998). Modelling of the hydraulic architecture of root systems: An integrated approach to water absorption—model description. *Annals of Botany*, 81(2):213-223.

Doussan, C., Pierret, A., Garrigues, E., and Pagès, L. (2006). Water uptake by plant roots: Ii - modelling of water transfer in the soil root-system with explicit account of flow within the root system - comparison with experiments. *Plant and Soil*, 283(1):99-117.

Dunbabin, V. M., Postma, J. A., Schnepf, A., Pagès, L., Javaux, M., Wu, L., Leitner, D., Chen, Y. L., Rengel, Z., and Diggle, A. J. (2013). Modelling root-soil interactions using three-dimensional models of root growth, architecture and function. *Plant and Soil*, 372(1/2):93-124.

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.

Huber, R. and Helmig, R. (2000). Node-centered finite volume discretizations for the numerical simulation of multiphase flow in heterogeneous porous media. *Computational Geosciences*, 4(2):141-164.

Huyakorn, P. S. and Pinder, G. F. (1983). *The Computational Methods in Subsurface Flow (Third Edition)*. Academic Press, third edition edition.

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

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.

Kelly, J., Barber, S., and Edwards, G. (1992). Modeling magnesium, phosphorus and potassium uptake by loblolly - pine seedlings using a barber-cushman approach. *Plant and Soil*, 139(2):209-218.

Khare, D., Selzner, T., Leitner, D., Vanderborght, J., Vereecken, H., and Schnepf, A. (2022). Root system scale models significantly overestimate root water uptake at drying soil conditions. *Frontiers in Plant Science*, 13.

Kirk, G. (1999). A model of phosphate solubilization by organic anion excretion from plant roots. *European Journal of Soil Science*, 50(3):369-378.

Koch, T. (2020). *Mixed-dimension models for flow and transport processes in porous media with embedded tubular network systems*. PhD thesis, University of Stuttgart.

Koch, T. (2022). Projection-based resolved interface 1d-3d mixed-dimension method for embedded tubular network systems. *Computers & Mathematics with Applications*, 109:15-29.

Koch, T., Gläser, D., Weishaupt, K., et al. (2021). DuMu\* 3 - an open-source simulator for solving flow and transport problems in porous media with a focus on model coupling. *Computers & Mathematics with Applications*, 81:423-443.

Koch, T., Heck, K., Schröder, N., Class, H., and Helmig, R. (2018). A new simulation framework for soil-root interaction, evaporation, root growth, and solute transport. *Vadose Zone Journal*, 17.

Koch, T., Wu, H., and Schneider, M. (2022). Nonlinear mixed-dimension model for embedded tubular networks with application to root water uptake. *Journal of Computational Physics*, 450:110823.

Köppl, T., Vidotto, E., Wohlmuth, B., and Zunino, P. (2018). Mathematical modeling, analysis and numerical approximation of second-order elliptic problems with inclusions. *Mathematical Models and Methods in Applied Sciences*, 28(05):953-978.

Landl, M., Huber, K., Schnepf, A., Vanderborght, J., Javaux, M., Bengough, A. G., and Vereecken, H. (2017). A new model for root growth in soil with macropores. *Plant and Soil*, 415(1-2):99-116.

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

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.

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(1-2):273-292.

Mesgouez, A. and Lefeuve-Mesgouez, G. (2009). Study of transient poroviscoelastic soil motions by semi-analytical and numerical approaches. *Soil Dynamics and Earthquake Engineering*, 29(2):245-248.

Meunier, F., Couvreur, V., Draye, X., Zarebanadkouki, M., Vanderborght, J., and Javaux, M. (2017a). Water movement through plant roots – exact solutions of the water flow equation in roots with linear or exponential piecewise hydraulic properties. *Hydrology and Earth System Sciences*, 21(12):6519-6540.

Meunier, F., Couvreur, V., Vanderborght, J., and Javaux, M. (2022). *How to Define the Appropriate Spatial Resolution of Root Segments When Solving Water Flow in Root System Hydraulic Architectures*. Springer.

Meunier, F., Draye, X., Vanderborght, J., Javaux, M., and Couvreur, V. (2017b). A hybrid analytical-numerical method for solving water flow equations in root hydraulic architectures. *Applied Mathematical Modelling*, 52:648-663.

Nowack, B., Mayer, K. U., Oswald, S. E., van Beinum, W., Appelo, C. A. J., Jacques, D., Seuntjens, P., Gérard, F., Jaillard, B., Schnepf, A., and Roosevelt, T. (2006). Verification and intercomparison of reactive transport codes to describe root-uptake. *Plant and Soil*, 285:305-321.

Pagès, L., Jordan, M., and Pcard, D. (1989). A simulation model of the three-dimensional architecture of the maize root system. *Plant Soil*, 119:147-154.

Pagès, L. L., Vercambre, G., Drouet, J.-L., Lecompte, F., Collet, C. C., and Le Bot, J. (2004). Root Typ: a generic model to depict and analyse the root system architecture. *Plant and Soil*, 258:103-119.

Pagès, L., Bécel, C., Boukcim, H., Moreau, D., Nguyen, C., and Voisin, A.-S. (2014). Calibration and evaluation of archisimple, a simple model of root system architecture. *Ecological Modelling*, 290:76-84.

Postma, J. A., Kuppe, C., Owen, M. R., Mellor, N., Griffiths, M., Bennett, M. J., Lynch, J. P., and Watt, M. (2017). Opensimroot: widening the scope and application of root architectural models. *New Phytologist*, 215(3):1274-1286.

Ramière, I. and Helfer, T. (2015). Iterative residual-based vector methods to accelerate fixed point iterations. *Computers & Mathematics with Applications*, 70(9):2210-2226.

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

Roose, T., Fowler, A., and Darrah, P. (2001). A mathematical model of plant nutrient uptake. *J Math Biol*, pages 347-360.

Rosenzweig, C., Elliott, J., Deryng, D., Ruane, A. C., Müller, C., Arneth, A., Boote, K. J., Folberth, C., Glotter, M., Khabarov, N., Neumann,
K., Piontek, F., Pugh, T. A. M., Schmid, E., Stehfest, E., Yang, H., and Jones, J. W. (2014). Assessing agricultural risks of climate change in the 21st century in a global gridded crop model inter-comparison. *Proceedings of the National Academy of Sciences*, 111(9): 3268–3273.

Samal, D., Kovar, J. L., Steingrobe, B., Sadana, U. S., Bhadoria, P. S., and Claassen, N. (2010). Potassium uptake efficiency and dynamics in the rhizosphere of maize (*Zea mays* l.), wheat (*Triticum aestivum* l.), and sugar beet (*Beta vulgaris* l.) evaluated with a mechanistic model. *Plant and Soil*, 332:105–121.

Sander, O., Koch, T., Schröder, N., and Flemisch, B. (2017). The dune foamgrid implementation for surface and network grids. *Archive of Numerical Software*, Vol 5.

Schäfer, E. D., Owen, M. R., Postma, J. A., Kuppe, C., Black, C. K., and Lynch, J. P. (2022). Simulating crop rootrootssystems using opensim-root. In Lucas, M., editor, *Plant Systems Biology: Methods and Protocols*, pages 293–323, New York, NY. Springer New York.

Schlömer, N. (2022). pygmsh: A python frontend for gmsh.

Schnepf, A., Black, C. K., Couvreur, V., Delory, B. M., Doussan, C., Koch, A., Koch, T., Javaux, M., Landl, M., Leitner, D., Lobet, G., Mai, T. H., Meunier, F., Petrich, L., Postma, J. A., Priesack, E., Schmidt, V., Vanderborght, J., Vereecken, H., and Weber, M. (2020). Call for participation: Collaborative benchmarking of functional-structural root architecture models. the case of root water uptake. *Frontiers in Plant Science*, 11:316.

Schnepf, A., Leitner, D., Bodner, G., and Javaux, M. (2022). Editorial: Benchmarking 3d-models of root growth, architecture and functioning. *Frontiers in Plant Science*, 13.

Schnepf, A., Leitner, D., Landl, M., Lobet, G., Mai, T. H., Morandage, S., Sheng, C., Zörner, M., Vanderborght, J., and Vereecken, H. (2018). CRootBox: a structural-functional modelling framework for root systems. *Annals of Botany*, 121(5):1033–1053.

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.

Simunek, J., Huang, K., and Van Genuchten, M. T. (1995). The swms_3d code for simulating water flow and solute transport in three-dimensional variably-saturated media. *US Salinity Laboratory Agricultural Research Service*, 139.

Somma, F., Hopmans, J., and Clausnitzer, V. (1998). Transient three-dimensional modeling of soil water and solute transport with simultaneous root growth, root water and nutrient uptake. *Plant and Soil*, 202(2):281–293.

Steefel, C., Yabusaki, S., and Mayer, K. (2015). Reactive transport benchmarks for subsurface environmental simulation. *Computational Geosciences*, 19.

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

Vanderborght, J., Kasteel, R., Herbst, M., Javaux, M., Thiéry, D., Vanclilloster, 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.

Zhou, X.-R., Schnepf, A., Vanderborght, J., Leitner, D., Lacointe, A., Vereecken, H., and Lobet, G. (2020). CPlantBox, a whole-plant modelling framework for the simulation of water- and carbon-related processes. *in silico Plants*, 2(1). diaa001.