Comparison of numerical methods for radial solute transport to simulate uptake by plant roots

Highlights

• The maximum spatial and time step are dependent on root and soil parameters.

• Transport problems in the rhizosphere can become stiff.

• The Crank-Nicolson method sometimes oscillated and was often slow.

• Adaptive Runge-Kutta methods were universally applicable.

Abstract

The 1D radial solute transport model with non-linear inner boundary condition is widely used for simulating nutrient uptake by plant roots. When included into an architectural root model, this local model has to be solved for a high number of root segments, e. g. ${10}^5$ – ${10}^6$ segments for large root systems. Each root segment comes with its own local parameter set in heterogeneous root architectural models. Depending on the soil and solute, the effective diffusion coefficient spans over more than six orders (e. g. for N, K, and P). Thus a numerical implementation of this rhizosphere transport model is required to be fast, accurate and stable for a large parameter space. We apply 13 methods to this rhizosphere model with root hairs and compare their accuracy, computational speed, and applicability. In particular, the Crank-Nicolson method is compared to higher-order explicit adaptive methods and some stiff solvers. The Crank-Nicolson method sometimes oscillated and was up to a hundred times slower than an explicit adaptive scheme with similar accuracy. For a given spatial resolution, Crank-Nicolson had about one order lower accuracy as other tested methods. The maximum spatial time step can be estimated from root radius, solute diffusion, advection, and soil buffer power. Although Crank-Nicolson is a viable method and often used as de-facto standard method for rhizosphere models, it was not the best performer in our comparison. While the best method remains problem specific, for general use in root architectural models we recommend adaptive Runge-Kutta with cubic or quadratic upwind for advection.

Introduction

The radial rhizosphere model (Barber and Cushman, 1981) has been widely used for simulating transport of various nutrients and their uptake by plant roots of many species. It was for example applied to phosphorus uptake by maize (Macariola-See et al., 2003), zinc uptake by rice crop (Adhikari and Rattan, 2000), magnesium, phosphorus, and potassium uptake by loblolly pine seedlings (Kelly et al., 1992) and N, P, K, Ca, Mg uptake by black cherry, northern red oak, and red maple (Kelly et al., 2001). Solving an individual instance of such a model typically takes a fraction of a second on modern CPUs. Advanced investigations in root biology require to take the soil heterogeneity, root system architecture and the variation in root parameters along that architecture into account. This can typically be done by replacing the simple root growth model with a 3D root architectural model, which discretizes the root system into millions of root segments, each with its own parameter set and instance of the radial rhizosphere model. This approach is for example implemented in Mai et al. (2019) and in the functional structural plant model OpenSimRoot (Postma et al., 2017), which was used in several studies on the utility of root architectural and anatomical traits for uptake of phosphorus (Postma and Lynch, 2011, 2012; Postma et al., 2014). Thus, for simulations with 3D functional structural plant models such as OpenSimRoot, computational speed of the numerical implementation can become critical. The rhizosphere model for a root segment consists of a 1D radial advection-diffusion equation with a Michaelis-Menten uptake kinetics at the inner boundary (rhizoplane) and zero flux outer boundary condition (Barber and Cushman, 1981). Itoh and Barber (1983) added a reaction term describing uptake by root hairs. Though analytical solutions have been identified for various special cases of transport in the rhizosphere (Cushman, 1979a, 1979b, 1980a, 1980b; Van Genuchten, 1981; Roose et al., 2001; Roose and Kirk, 2009; Ou, 2019), there is no analytical closed-form solution in the general case. Thus, in addition to computational speed, decisions on which method to use for an efficient numerical implementation should also be made on stability and accuracy. A certain minimum accuracy is desired over a range of model parameter values which typically is rather wide. For example, Bouldin (1961) considered diffusion rates from $1\times {10}^{-5}$ to $1\times {10}^{-9}{\text{cm}}^2{\text{s}}^{-1}$ for phosphorus, however, without soil sorption, and root radii from $5\times {10}^{-5}$ to $7.5\times {10}^{-4}$ cm in a pure diffusion model without root competition. Relevant model parameters that determine nutrient depletion are soil parameters like solute diffusion, advection, soil buffer power, along with root morphology, and uptake. Uptake by roots is governed by root radius, $r_0$, and Michaelis-Menten parameters describing the relation between concentration and nutrient uptake. Together, these parameters determine the concentration profiles and thus the suitability of specific numerical methods. How the parameter values relate to required minimum step sizes is hinted at by characteristic dimensionless numbers such as the Péclet number (e. g. Huysmans and Dassargues, 2005), Courant number (CFL: Courant, Friedrichs, and Lewy, 1928) and Fourier number (e. g. Kunes, 2012).

Plants not only absorb but also exude solutes. The model can be adapted for this by changing the inner boundary condition and introducing a reaction term to represent degradation rates of the (organic) exudates. Such adaptations benefit from flexibility in the implementation. Thus besides speed, accuracy, and stability, another requirement on the numerical implementation of the model is that it can be easily extended. Passioura and Frere (1967), Nye and Marriott (1969), Newman and Watson (1977), Nye (1981), Barber and Cushman (1981), Itoh and Barber (1983) used the (implicit) Crank-Nicolson method (Crank and Nicolson, 1947) to solve macroscopic rhizosphere nutrient transport and this method is still popular to solve such parabolic equations in rhizosphere models (e. g. Boghi et al., 2018). Other commonly used methods are explicit forward simulations like a forward Euler method (FTCS: forward time central space) (e. g. Leadley et al., 1997, Roshania et al., 2009), which needs, in general, a much smaller time step size than the implicit Crank-Nicolson method (CN) to be stable and achieve comparable accuracy (for non-smooth solutions). Ou (2019), Roose and Kirk (2009) even used first-order (upwind) discretization for the first derivative in space and time, second-order only for the second derivative, and used the numerical results as a reference for the analytical approximations. Explicit methods might need smaller time steps but are straight forward to implement and easy to understand, which are helpful attributes when extending the model. Numerical methods are a much overlooked topic in rhizosphere modeling in general. To our knowledge, so far no comprehensive testing of numerical methods for a classical rhizosphere model has been presented. Here, we investigated several explicit and implicit methods varying in their numerical order in space and time, applied to the classical rhizosphere model with root hairs by Itoh and Barber (1983), concerning the mentioned criteria: speed, accuracy, and applicability. We compared Crank-Nicolson to various explicit Runge-Kutta methods and implicit methods, BDF (Backward Differentiation Formulas) and IRK (implicit Runge-Kutta), with adaptive time steps. The aim was to identify methods that are most suitable for linking to a root architectural model, by fulfilling the criteria for a wide range of parameters, especially a high variation in the proportion of advection and effective diffusion. Numerical experiments with various parameters were performed to estimate the numerical error and computational time and thus to evaluate the different methods.