How to put plant root uptake into a soil water flow model

The need for improved crop water use efficiency calls for flexible modeling platforms to implement new ideas in plant root uptake and its regulation mechanisms. This paper documents the details of modifying a soil infiltration and redistribution model to include (a) dynamic root growth, (b) non-uniform root distribution and water uptake, (c) the effect of water stress on plant water uptake, and (d) soil evaporation. The paper also demonstrates strategies of using the modified model to simulate soil water dynamics and plant transpiration considering different sensitivity of plants to soil dryness and different mechanisms of root water uptake. In particular, the flexibility of simulating various degrees of compensated uptake (whereby plants tend to maintain potential transpiration under mild water stress) is emphasized. The paper also describes how to estimate unknown root distribution and rooting depth parameters by the use of a simulation-based searching method. The full documentation of the computer code will allow further applications and new development.


Introduction
Increased biomass and yield of agricultural crops hinges on improved efficiency of root water uptake 1,2 . Biophysical mechanisms responsible for the improved uptake efficiency in turn rely on linkages between plant physiological processes and the physics of soil water movement; the development of the latter has been well illustrated in 3,4 and 5. Research literature in this and closely related areas has been profound in past decades, though roughly two categories of approach can be identified: one that is mechanistic and one that is phenomenological in nature. The mechanistic approach has improved our understanding of the soil-root interactions by identifying major resistances and driving forces of water flow from soil to roots 3,4 , as well as the scale-dependent interactions and geometry among soil organisms responsible for water and nutrient uptake 6,7 . Yet effective mechanistic modeling needs a large amount of empirical data support, which can not always be met due to various factors constraining research activities. On the other hand, an intuitive, phenomenological approach can be implemented with less empirical data support. Examples of this latter approach may include constraining root water uptake capacity by soil water potential 8 , by root length distribution 9-11 , and through additional modifications of the root uptake sensitivity to local soil water potential and root density [12][13][14][15][16] , which allows compensated uptake to maintain potential transpiration 17,18 .
One of the major goals of agronomic management for increased water use efficiency is to channel water loss through transpiration 19 . This to some extent can be aided by the compensated root uptake, which has been implemented in software such as HYDRUS (http://www.pc-progress.com/en/Default.aspx? hydrus-1d) through the use of a threshold root uptake adaptation factor, whereby reduced water uptake in water-stressed parts of the root zone is fully compensated for by increased uptake in other soil regions where water is available 20 .
However, root water uptake compensation has been implemented in different studies using different models, which naturally calls for the flexibility of implementation of various forms of the compensation mechanism in root uptake studies. The first objective of this paper is to show that this flexibility can be obtained by using a water flow model extended from an original model as described in 5. This extended model was used in generating numerical solutions of water flow problems in 16.
The second objective of this paper is to provide a tool for quantifying plant water uptake under conditions where information of root depth and distribution pattern is not known. This is relevant to agricultural applications as fine root information remains difficult to obtain despite the progress made in advanced underground sensing 21 . We take advantage of inverse modeling, similar to 22 and 10, to find optimal solutions of unknown root parameters based on good measurement of soil water content, which can be obtained through the use of various soil moisture sensors. Despite excellent documentation of the methods for solving the soil water flow equations relevant to agricultural water management 3,5 , the author has not seen a paper documenting the details of how to put the root uptake component into a soil water flow model, as well as providing the details of a search-based optimization procedure for root parameter estimation. This article provides a further documentation of the methods as used in 16, with the focus of highlighting details related to the computer implementation of the numerical methods.

Implementation
The one dimensional Richards equation of soil water flow in a vegetated soil in the vertical direction can be written as where h is water pressure head / soil matric potential (m), t is time (s), z space (m, positive downward and soil surface as zero), K hydraulic conductivity (m/s), and θ ∂ = ∂ c h soil water capacity, and S(z,t) the rate of plant root water uptake at location and time. The reason why there is a negative sign in Equation 1 is that, when we take the space derivative for the gravity potential head z, namely, in ∂z/∂z, the value in the numerator changes in the opposite direction as that in the denominator: gravity head always decreases downwards but our chosen space z increases downwards. Equation 1 can be solved using the implicit finite difference method as described by 5. Here is a summary of key points with a few intuitive comments. Let's for now assume surface fluxes as precipitation, irrigation, or surface evaporation. As Equation 1 considers soil water dynamics both in space (vertical axis) and time, we need to properly discretize the equation considering both space and time, in order to solve it numerically. Here we let all h's that take space derivatives be approximated using values at the old (t) or the new time (t +1), but let K's be approximated in old time (t). All other variables not taking partial derivatives are to be approximated in the old time. Assume the vertical axis is equally divided by Δz, to approximate Equation 1 at z = z i and t n+1 = t n +Δt, we have the following formula to use: * * * * * where S i,t refers to the root water uptake occurring at the ith block at current time t. In 16, S i,t was calculated as , ,

Amendments from Version 1
The major differences between this version of the article and the previously published version are reflected in (1) an update of Figure, (2) the addition of a sink term for root water uptake in all governing equations, including a further discussion of the contents related to the sink term, and (3) minor changes to correct a few typos and update on literature citations.
Any further responses from the reviewers can be found at the end of the article

REVISED
where α(h i,t ) (dimensionless) is the reduction factor of root uptake as a function of soil water potential at location i and time t, T t and L t are potential transpiration (mm/day) and root depth (m) at time t, z i is the vertical coordinate (m), and β is an empirical factor determining vertical distribution of root length density according to 12. The factor of 86400000 was used to ensure the unit of S i,t is second -1 (there are 86400 seconds in one day and one meter is equivalent to 1000 mm).
In Equation 2, soil matric potential at the new time (t+1) is indicated using a star as superscript. One interesting aspect of Equation 2 is that the space index for K is shifted to the upper edge of each of the blocks, while the index of h or θ is located at the center of the blocks (Figure 1). This is the characteristic of the block-centered finite difference scheme and it will enable the calculation of fluxes between two adjacent soil layers, or blocks, as seen in the right sub-figure of Figure 1, which is based on Figure 5-13 of Warrick's book (on page 195). In Figure 1, we reinforce the idea that, similar to the situation of conductivity, K, the index for flux density also is located at the edge of the blocks. Intuitively, this makes sense, because in the head-based equation of soil water flow, K and J have the same unit (both are in m/s). As a result, we can compare the magnitude of K and J directly. For example, if the water conductivity of the surface soil is similar to the rate of precipitation, we can roughly say that all the rainwater would have time to fully infiltrate into the soil before ponding on the soil surface.
The terms of Equation 2 can be collected to form a tridiagonal system of linear equations, with the water potentials for the new time (starred) arranged at the left side and those for the old time on the right side of the equation. Then, terms for different space segments can be collected to form the following equation: * * * − + which can be written in matrix form as Equation 5 is called tridiagonal because, for the coefficients of the nz by nz square matrix, the elements are all zero except for The complete coefficients of the tridiagonal system of the Richards' equation without a root water uptake sink were summarized in Table 5-7 of 5 (page 194). When the root uptake is considered, as in Equation 2, it turns out that only the D i terms for various upper boundary conditions need to be changed, while other terms (i.e., A i , B i , and C i ) remain the same as in Table 5-7 of 5. Specifically, the D i terms for different nodes (or blocks as shown in Figure 1) with root water uptake are shown below: • For internal nodes (i=2, …, nz-1): • For upper node with J(0,t): • For upper node with h(0,t): • For lower node with unit hydraulic gradient: nz t nz t nz t nz t nz nz t After pressure heads for the new time (h *' s) are computed using the Thomas algorithm, the flux of water between each pair of two neighboring blocks can be computed. However, to update the water content of the ith block, the following equation is used: To meet the requirement of mass balance, the water content of the ith block at the new time ( * i θ ) is caused by (1) the original water content at the previous time step (θ i,t ), (2) the contribution from the adjacent upper block where i =1, …, nz. To understand why some terms are positive and others are negative in Equation 11, we need to keep in mind (1) the mass conservation of the ith block (see Figure 1B), and (2) our choice of the vertical axis (downwards positive). For example, if the flux across the lower edge of the ith block is upward in direction, i.e., * 0.5 i J + is negative (because the upward direction is against our positive z axis), then the corresponding term in Equation 11 becomes positive. This is reasonable because an upward flux from the lower edge will add water to the ith block (again, see Figure 1B). Finally, the reason why there is a Δt in the sink term is that by multiplying S i,t with Δt, we get a soil water content (since the unit of S i,t is T -1 ).
Detailed solution procedures of the above model "Austere-Layered" (hereafter referred to as AL.f) are provided in 5. Here in order to let this soil water flow model simulate root uptake, the following new components have been added: • Assume we use the flux upper boundary condition (B_TYPE= 2). At each new time step, the old hydraulic conductivities and soil water capacities are used to find the needed coefficients for solving the tridiagonal system of equations. However, the surface flux is temporarily not considered (C2=0). Then, these coefficients are provided to the Thomas subroutine to find the updated pressure head (HSTAR). The trial soil water content at the new time (THETA) is found using the newly computed conductivities based on values of HSTAR as well as the surface flux.
• The trial value of THETA is accepted only if the soil water content at a depth interval, or that for intervals immediately adjacent to this current depth interval, is unsaturated. In such a case, the soil water retention curve is used to update the new pressure head H based on the trial value of THETA.
• Otherwise, the trial value of THETA is rejected, the value of HSTAR is accepted as the updated value for H, and the updated THETA is computed from the soil water retention curve based on HSTAR.
Testing different ending soil water potentials for optimal root uptake Detailed information regarding the root uptake sink term is included in 23. In particular, Equation 1 of 23 illustrates the importance of four water pressure heads for root uptake: h 1 marks the anaerobiosis point; h 4 the pressure head at permanent wilting point, and h 2 and h 3 marks the range of pressure head for optimal root uptake (see Figure 3 of reference 16). Here we refer to h 3 as the ending water potential for optimal uptake, because root uptake starts to reduce when the soil water potential becomes lower than this point. As different plants might have different sensitivity to a mild water stress, it is useful to test the model's sensitivity to h 3 .

Testing different functions for root water uptake
The extended program ALS.f also allows for testing different root uptake functions with or without compensation mechanisms, such as those demonstrated in 14,15,26. In particular, in subroutine SINK(), the model of 12 is used when ICPS= 1; a function of 9 is used when ICPS= 3; and the compensation function of 14 (in combination with Wu's function) is used when ICPS= 2. The normalized root density function (L nrd ) of 9 has four parameters: where z r = z/Z r (t) is the normalized root density at time t with Z r (t) the maximum rooting depth (m) at time t. The coefficients R 1 through R 4 represent experimentally fitted values for wheat, with R 1 = 2.21, R 2 = -3.72, R 3 = 3.46, and R 4 = -1.87 according to 9. The compensation function of 14 and 15 at time t has the following form: where S i is the actual water uptake from ith soil depth increment, α i is the reduction function for this ith depth increment, L nrd,i the normalized root density function for the ith depth increment, Δz i the thickness (m) of the ith depth increment, and λ is a fitting parameter varying from 0.01 to 2.0 according to 14.

Details of simulation-based searching program
In order to cope with situations in which the information of rooting depth and root distribution pattern is unknown, the program ALS.f was converted into a subroutine and is callable by a searching main program, ALSS.f, which is available at https://zenodo.org/record/42702. This conversion was quite straight forward in Fortran 77 in that we only did three things: (a) we deleted parameters BETA and ZRMAX from the parameter listing; (b) we converted the program ALS.f from a main program to a subroutine using 2 input dummy variables BETA and ZRMAX and 7 output dummy variables; and (c) we let this subroutine be called multiple times using different values of true BETA and true ZRMAX from arrays BETAS() and ZRMAX() within the main program unit and made sure that only the values of selected 7 variables were returned to the main program but all the other information from within the subroutine ALS() is lost after each run.
To run ALS.f and ALSS.f, the following input information is needed (see 23): • Soil physical properties and in particular soil water retention curve and saturated hydraulic conductivity. The exact format of data structure is shown in the supplementary data file ALS.DAT, which follows the template provided by 5. A method for predicting soil water retention curves based on easily measured soil texture and bulk density data is illustrated in 27.
• Initial soil water content profile. This must be provided and put to exact format as shown in supplementary data file INITF.DAT, following the template provided by 5. The measured soil water content must be interpolated so that every depth increment of soil has a water content value in it. One method of data interpolation is outlined in 16 using routines provided in 28 and 29.
• Measured terminating soil water content. This can be the measured soil water content on last day of simulation and must be entered in the format as shown in TERMF.DAT of 23. No interpolation is needed for this file.
• Daily weather data. Daily evapo-transpiration rate, daily average air temperature, daily rainfall, relative humidity, as well as leaf area index, must be entered in data file WEA2.DAT in the order as shown in the data file. Daily leaf area index can be interpolated based on field measured values, as shown in 16.
There are four output files from ALS.f 23 : • OUTS.DAT -Mirrors output from AL.f of 5.
• OUTS2.DAT -Daily soil water potential and water contents for interested soil depths.
• WUE.DAT -Summary of daily cumulative water loss.
There are seven output files from ALSS.f. They provide the following information based on particular realization of combinations of root distribution shape factor BETAS and the maximum rooting depth ZRMAXS 23 : • ARD.DAT -Average relative discrepancy between measured and simulated soil water content throughout the simulation period.
• WBELOW.DAT -Amount of water to deep drainage out of the specified root zone. This was inherited from AL.f • WADD.DAT -Amount of water added to the system. This was also inherited from AL.f • BIGI.DAT -Cumulative infiltration. This was also inherited from AL.f • PCTDIF.DAT -Percent difference in water balance, defined in 5 as (BIGI-WADD)/WADD*100.
The main difference between the output files of ALS.f and those of ALSS.f is that the outputs from the former are a summary of water balance from the simulation of soil water balance based on one set of root parameters (rooting depth and shape factors), while those from the latter (ALSS.f) contain the simulation outcomes using many sets of root parameters (see next section for further details).

Use cases
Programs ALS.f and ALSS.f in their current forms were customized to simulate the soil water dynamics for 111 days from May 14 to September 1, 2009 in a mixed-grass prairie ecosystem 16 . The model was also applied to simulate soil water dynamics of a same grassland under long-term heavy grazing by cattle 30 . This is reflected both in the programs and the input data files. To use the model in other locations and other durations, appropriate changes must be made both in a few parameters of the programs and input data files. The two potential uses of the programs are outlined in the following: • If root distribution parameter BETA and maximum rooting depth ZRMAX are known, then program ALS.f may be used to simulate soil water dynamics with root uptake. Make sure to provide updated values for the major parameter values, such as JMATUR (Julian day number of growth maturity), JTHAW (Julian day number of growth initiation), BETA, and ZRMAX.
• If root distribution parameter BETA and maximum rooting depth ZRMAX are unknown, then it is recommended to first use ALSS.f to find the optimal values for BETA and ZRMAX before using ALS.f to simulate soil water dynamics. Before running ALSS.f, three things must be done: • Finally, it is recommended to put both the program files and the required input data files in the same folder before running the programs.

Summary
While plant root water uptake has been implemented in commercial software packages, very few researchers have found the time to document important details of how root uptake can be included in a soil water flow model. This paper tries to fill this knowledge gap by using a simple computer program aided with sample input data sets. Specifically, we modified and extended a soil water flow model of 5 to include plant root water uptake capability. We showed that the extended program has the flexibility of considering different sensitivity of root uptake, as well as different water uptake mechanisms. We also modified the model to search for optimal root parameters based on measured soil water content data. This is relevant for new application in that field data of root development usually is expensive and sometimes difficult to obtain, and the utilization of a searching-based modeling approach will help to identify root parameters for improved water balance characterization, which in turn may help to identify ways to improve crop water use efficiency under field conditions. No menu-driven user interface is provided but the author believes that, under some circumstances of knowledge discovery, it is more efficient to work directly with a well-structured computer program, such as the one described in this paper, than relying on a packaged one.  Tables 1 and 3) and those added in ALS.f (Supplementary Tables 2 and 4).

Hirotaka Saito
Graduate School of Agriculture, Tokyo University of Agriculture and Technology, Fuchu, Japan The article entitled "How to put plant root uptake into a soil water flow model" by Dong is in general well written and technically sound. However when I read this manuscript, I just found it a bit difficult to read as most relevant information such as Feddes function was in a supplemental file. There were no examples in the manuscript either. In my opinion, if the author is providing a tool to evaluate root water uptake, several examples (more than one) need to be included in the manuscript; for example how the choice the model affect soil water dynamics. This fact really makes it difficult to evaluate the manuscript. There should be a sink term, S, that is used to account for root water uptake in Eq (1).
○ Figure 1 needs to be upgraded. They get blurred.

○
It is an interesting topic and the text is well written but it is requested some program experience. It is a technical report not fully a scientific paper. It gave sufficient detail information about this model and also some relative programs. If a plant root uptake model with plant growth process which can report LAI or other parameters for ET calculation was put to Soil water flow model, this is more useful. If some cases can be tested and compared with other models, it will make user to use it easily. I fully agreed that it can be published as so far and give more opportunities to test it.