Keywords
Computer simulation, Richards’ equation, Root growth, Soil water balance, Transpiration, Uptake compensation, Water stress
This article is included in the Agriculture, Food and Nutrition gateway.
This article is included in the Plant Science gateway.
Computer simulation, Richards’ equation, Root growth, Soil water balance, Transpiration, Uptake compensation, Water stress
Increased biomass and yield of agricultural crops hinges on improved efficiency of root water uptake1,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 roots3,4, as well as the scale-dependent interactions and geometry among soil organisms responsible for water and nutrient uptake6,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 potential8, by root length distribution9–11, and through additional modifications of the root uptake sensitivity to local soil water potential and root density12–16, which allows compensated uptake to maintain potential transpiration17,18.
One of the major goals of agronomic management for increased water use efficiency is to channel water loss through transpiration19. 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 available20.
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 sensing21. 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 management3,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.
The one dimensional Richards equation of soil water flow in the vertical direction can be written as
where h is water pressure head (m), t is time (s), z space (m, positive downward and soil surface as zero), K hydraulic conductivity (m/s) and 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 values in the numerator change in the opposite direction as that in the denominator: gravity head always decreases downwards but our chosen space z increases downwards.
The model described by Equation 1 does not consider a sink of plant root water uptake. It 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. Here we let all h’s that take space derivatives be approximated using values at the old (t = n) and the new (t = tn+1) times, but let K’s be approximated in old time (t = n). 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 = zi and tn+1 = tn+Δt, we have the following formula to use:
Another 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. 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 rain water would have time to fully infiltrate into the soil before ponding on the soil surface.
A. Pressure heads and water conductivities are indexed differently. B. Highlight of one internal block. Modified based on Figure 5–13 of 5.
The terms of the implicit scheme 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 4 is called tridiagonal because, for the coefficient of the nz by nz square matrix, the elements are all zero except for those located along the three diagonals. From Equation 4, we can see that A1 = Cnz = 0. Actually, values of the coefficients on the left side of Equation 4 should be determined considering the boundary conditions at the soil surface and at the bottom of the root zone.
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: (a) dynamic root growth, (b) non-uniform root distribution and water uptake, (c) effect of water stress on plant water uptake, and (d) soil evaporation. Following Warrick’s original model, the extended model (ALS.f) was also programmed using Fortran 77. See https://zenodo.org/record/42663. To facilitate further development, variables (array variables and non-array variables) used in the original program AL.f, as well the extended program ALS.f, are listed in Supplementary material 1 (in four tables). The subroutine SINK() in ALS.f encapsulates various options for root uptake mechanisms. Similar to some other external subroutines, such as THOMAS(), variables inside the subroutine are local, so that they are not listed in the tables of Supplementary material 1. In the following we document major changes made to the program ALS.f:
The model AL.f was modified to include the flux update scheme proposed by 23 and was applied in program a2&3.for of 5. The main idea is summarized in the following (variables referred to are included in ALS.f):
• 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.
Detailed information regarding the root uptake sink term is included in Supplementary material 2. In particular, Equation 1 of Supplementary material 2 illustrates the importance of four water pressure heads for root uptake: h1 marks the anaerobiosis point; h4 the pressure head at permanent wilting point, and h2 and h3 marks the range of pressure head for optimal root uptake. Here we refer to h3 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 h3.
Following 14 and 24, h3 may be assumed as for a slow transpiration of and for a fast transpiration of We tested two more possibilities of the h3 values, namely, and for We also assumed that h3 changed in a linear fashion for intermediate transpiration rates between and using different and values:
where This change is reflected in the new subroutine HTHREE, as shown in ALS.f. This (Equation 5) is another form of Equation 2 of Supplementary material 2 but with variable and .
The extended program ALS.f also allows testing different root uptake functions with or without compensation mechanisms, such as those demonstrated in 14,15,25. In particular, in subroutine SINK(), the model by 12 is used when ICPS= 1; a function by 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 (Lnrd) of 9 has four parameters:
The extended program ALS.f also allows testing different root uptake functions with or without compensation mechanisms, such as those demonstrated in 14,15,25. In particular, in subroutine SINK(), the model by 12 is used when ICPS= 1; a function by 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 (Lnrd) of 9 has four parameters:
where zr = z/Zr(t) is the normalized root density at time t with Zr(t) the maximum rooting depth (m) at time t. The coefficients R1 through R4 represent experimentally fitted values for wheat, with R1 = 2.21, R2 = -3.72, R3 = 3.46, and R4 = -1.87 according to 9. The compensation function proposed by 14 at time t has the following form:
where Si is the actual water uptake from ith soil depth increment, αi is the reduction function for this ith depth increment, Lnrd,i the normalized root density function for the ith depth increment, Δzi the thickness (m) of the ith depth increment, and λ is a fitting parameter varying from 0.01 to 2.0 according to 14.
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 after each run but all the other information from within the subroutine ALS() is lost.
The provided two Fortran programs, ALS.f and ALSS.f, can be run from one of the freely available compilers, such as g95 (http://www.g95.org/). A short introduction of how to install and use g95 to a PC under Window XP or Windows 7 can be found at: http://www.ems.psu.edu/~young/meteo473/Bootcamp/g95Instructions.pdf
To run ALS.f and ALSS.f, the following input information is needed (see Supplementary material 3):
• 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 26.
• 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 27 and 28.
• 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 Supplementary material 3. No interpolation is needed for this file.
• Daily weather data. Daily evapo-transpiration rate, daily average air temperature, daily rainfall, relatively 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 result, as shown in 16.
There are four output files from ALS.f (Supplementary material 3):
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 (Supplementary material 3):
• ARD.DAT – Average relative discrepancy between measured and simulated soil water content throughout the simulation period.
• C_AEVA.DAT – Cumulative soil evaporation throughout the simulation period.
• C_ATRA.DAT – Cumulative transpiration 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).
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. 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:
1. Check to see if the current value of h3 (see Equation 1 of Supplementary material 2) is correct by changing, if necessary, the value of parameter IHTH (current default value is set to 1, meaning See the preamble of subroutine SINK() for two other options.
2. Check to see if the parameter of ICPSN is correctly specified. currently, the default value is 1, meaning to use Ojha’s compensation function. Setting the value of ICPSN to 2 will switch to Li’s compensation function, and a value of 3 will invoke Wu’s uptake function (without compensation).
3. Check the nested DO LOOPs structure to adjust the admissible values for index variable II (currently it is set to have 16 different values for BETA) and JJ (currently it is set to have 22 different values for ZRMAX). At the same time, the following lines in the main program of ALSS.f must be changed to be compatible with the range of index values in the nested DO LOOPS.
BETAS(1)=0.5
DO I=1,15
BETAS(I+1)=BETAS(I)+0.25
END DO
ZRMAXS(1)=0.75
DO I=1,21
ZRMAXS(I+1)=ZRMAXS(I)+0.05
END DO
• Finally, it is recommended to put both the program files and the required input data files in the same folder before running the programs.
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.
ALS.f is available from: https://zenodo.org/record/42663
ALSS.f is available from: https://zenodo.org/record/42702
ALS.f archived source code as at the time of publication: http://dx.doi.org/10.5281/zenodo.4266329
ALSS.f archived source code as at the time of publication: http://dx.doi.org/10.5281/zenodo.4270230
License: Academic free license ("AFL") v.3.0 http://opensource.org/licenses/afl-3.0.php
This paper is part of North Dakota Agricultural Experiment Station project ND6149 and Texas A&M AgriLife Research Hatch project TEX09574, both funded by the USDA NIFA.
I confirm that the funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
The author thanks Bob Patton, Paul and Anne Nyren and Lyle Prunty at North Dakota State University for fruitful collaboration leading to this paper.
Supplementary material 1. Definition of variables used in the computer programs ALS.f.
This file contains four tables documenting the variables used in AL.f (Supplementary Tables 1 and 3) and those added in ALS.f (Supplementary Tables 2 and 4).
Click here to access the data.
Supplementary material 2. More details about the root uptake term.
This file documents details of the water uptake term as used in the extended Fortran program ALS.f. It is useful for those who are curious about how the water extraction term is specified in the model.
Click here to access the data.
Supplementary material 3. Data files for ALS.f and ALSS.f.
This includes (a) four input data files common to ALS.f and ALSS.f, (b) four Output files from ALS.f (OUTS.DAT, OUTS2.DAT, INFS.DAT, WUE.DAT), (c) seven output files from ALSS.f (ARD.DAT, C_AEVA.DAT, C_ATRA.DAT, WBELOW.DAT, WADD.DAT, BIGI.DAT, PCTDIF.DAT).
Views | Downloads | |
---|---|---|
F1000Research | - | - |
PubMed Central
Data from PMC are received and updated monthly.
|
- | - |
Competing Interests: No competing interests were disclosed.
Competing Interests: No competing interests were disclosed.
Competing Interests: No competing interests were disclosed.
Alongside their report, reviewers assign a status to the article:
Invited Reviewers | |||
---|---|---|---|
1 | 2 | 3 | |
Version 2 (revision) 08 Aug 22 |
|||
Version 1 08 Jan 16 |
read | read | read |
Provide sufficient details of any financial or non-financial competing interests to enable users to assess whether your comments might lead a reasonable person to question your impartiality. Consider the following examples, but note that this is not an exhaustive list:
Sign up for content alerts and receive a weekly or monthly email with all newly published articles
Already registered? Sign in
The email address should be the one you originally registered with F1000.
You registered with F1000 via Google, so we cannot reset your password.
To sign in, please click here.
If you still need help with your Google account password, please click here.
You registered with F1000 via Facebook, so we cannot reset your password.
To sign in, please click here.
If you still need help with your Facebook account password, please click here.
If your email address is registered with us, we will email you instructions to reset your password.
If you think you should have received this email but it has not arrived, please check your spam filters and/or contact for further assistance.
Comments on this article Comments (0)