Root water uptake

Author

Affiliation

Henning Kage

Institute of Crop Science and Plant Breeding, Agronomy and Crop Science, Christian-Albrechts-University, Kiel, Germany

Keywords

Root water uptake, Single root model, Sinkt term

1 Introduction

Water flows from the soil to the root surface along a gradient of total water potential. Assuming that besides the matric potential other components of the total potential are constant and the soil is homogeneous, it follows that the water content at the surface of a root taking up water has to be lower than the water content of the bulk soil. According to Darcy’s Law, the difference in water content between the bulk soil and the root surface is enlarging with an increasing transpiration rate and with a decreasing hydraulic conductivity. The hydraulic conductivity itself is a function of the soil water content and is influenced by the texture and structure of the soil. Since a crop with a low total root length has to achieve higher water uptake rates per unit root length at a given evaporative demand, it is also likely that the draw-down in water content near the roots is getting larger with decreasing root length. The water content distribution around roots is therefore a function of the soil hydraulic properties, the root length density, and the transpiration rate. The water content distribution around roots can be simulated with different approaches. The most simple approach is to assume a constant soil water diffusivity with respect to the radial distance to the root surface.

It has been shown, that for low root length densities a limitation of water extraction rates from a finite volume of soil may exist under conditions of high water uptake rates per unit root length and for low soil hydraulic conductivities (Kage and Ehlers (1996)) even in cases were the average soil water content is well above the permanent wilting point (). Besides the others, the pre crop can have substantial effects on total root length and root length density (). Higher susceptibility to drought stress may the consequence.

One common simplifying assumption in most models of root water uptake is that the soil water diffusivity is constant with respect to the radial distance to the root surface. This assumption is, however, not correct, as the soil water diffusivity decreases with decreasing soil water content. The decrease in soil water diffusivity near the root surface leads to a self-reinforcing effect, as the decrease in soil water diffusivity leads to a decrease in the water content near the root surface, which in turn leads to a further decrease in soil water diffusivity. The decrease in soil water diffusivity near the root surface is, however, not considered for in most models of root water uptake.

Another simplification is that roots are assumed to be distributed evenly within a certain soil layer, leading to the specification of a single “typical single root cylinder” as a geometrical model for the root system, for which calculations of water transport to roots in simple 1-D radial coordinates are sufficient. This simplification is, however, strictly speaken not correct, as roots may be distributed unevenly within the soil even within a specified depth. The uneven distribution of roots within the soil may lead to a non-uniform water content non uniform water extraction.

Despite the fact that recent computational power allows for the simulation of the water uptake from a root system in a three-dimensional explicit approach, this approach is still very time consuming and the parameterisation of the root architectural models is still a challenge.

This therefore motivates the development of simplifying approaches, which, however, take into account a) the decrease of soil water conductivity near the roots as a consequence of the decrease of the soil water content due to the water uptake itself and b) the non-uniform distribution of roots within the soil. The aim of this study is to compare different approaches to simulate root water uptake and maximum rates of water extraction. The first approach is based on the analytical solution of the radial water transport problem, which was first given by Philip (1957) and Gardner (1960). The second approach is based on the steady state and steady rate approximation for finite soil volumes, which was first given by Passioura and Cowan (1968). The third approach is based on the matrix flux potential, which was first given by van Lier et al. (2006). A numerical solution of the partial differential equation describing the water transport in radial coordinates towards a water absorbing root surface was uses as a reference ((hillel1975?), Kage and Ehlers (1996)).

1.1 Comparison of approaches to simulate root water uptake and maximum rates of water extraction

The script compares different approaches to simulate root water uptake. The first approach is based on the analytical solution of the radial water transport problem, which was first given by Philip (1957) and Gardner (1960). The second approach is based on the steady state and steady rate approximation for finite soil volumes, which was first given by Passioura and Cowan (1968). The third approach is based on the matrix flux potential, which was first given by van Lier et al. (2006). A numerical solution of the partial differential equation describing the water transport in radial coordinates towards a water absorbing root surface was uses as a reference ((hillel1975?), Kage and Ehlers (1996)). The script compares the water content distribution around roots at different soil water potentials and for three different soil textures calculated with the steady rate model.

1.2 Basics for calculation water transport towards roots

The continuity equation for soil water (\theta) in radial coordinates is: \frac{\partial \theta }{\partial t}=-\frac{\partial q}{\partial r} \tag{1}

where q is the flux density [cm^3 \cdot cm^{-2} \cdot d^{-1}] of water towards the root surface and r is the distance from the root surface. The combination of the continuity equation with the Richards equation for radial coordinates yields:

\frac{\partial \theta }{\partial t}=\frac{1}{r}\left( {{r}\cdot}{{k}_{u}}(\theta )\cdot\frac{\partial \psi }{ \partial r} \right) \tag{2}

Using the identity \frac{d\theta }{d\theta } this can be written as:

\frac{\partial \theta }{\partial t}=\frac{1}{r}\left( r\cdot {{k}_{u}}\left( \theta \right)\cdot \frac{\partial \psi }{\partial \theta }\cdot \frac{\partial \theta }{\partial r} \right) \tag{3}

The specific soil water capacity is defined as: {{C}{(\psi }}{)}=\frac{d\theta }{d\psi } \tag{4}

Combined with the unsaturated soil hydraulic conductivity this with gives the definition of soil water diffusivity:

{{D}_{w}}(\theta )={{k}_{u}}(\theta )\cdot \frac{1}{{{C}_{(\psi )}}} \tag{5}

which can be used to rewrite equ. Equation 3 as:

\frac{\partial \theta }{\partial t}=\frac{1}{r}\cdot \left( r\cdot {{D}_{w}}\left( \theta \right)\cdot \frac{\partial \theta }{\partial r} \right) \tag{6}

For this partial differential equation analytic solutions exists, but only for special cases of the boundary conditions and for specific forms of the soil specific relationship for D_w(\theta).

1.3 Analytical solutions of the radial water transport problem

Gardner (1960) and Philip (1957) first gave analytic solutions for boundary conditions given by (a) a constant rate (I_w) of water uptake per unit root length [cm^3cm^{-1}d^{-1}] where the inner boundary condition is given by:

I_w = -2\pi aK \frac{\delta\psi}{\delta r} \tag{7}

This can also be written in terms of soil water diffusivity and with a flux density [d^{-1}] on the left hand side as the term 2\pi a is the root surface area per unit root length [cm^2]:

\frac{I_w}{2\pi a} = -D_w \frac{\delta \theta}{\delta r}=q \tag{8}

at r = a, and with the assumption of an infinite soil volume available for the sink, the solution is Philip (1957), Carslaw and J.C. Jaeger (1946)):

\overline{\psi_s}-\psi_a=\frac{I_w}{4\pi K_u}\left[ln\left(\frac{4D_wt}{a^2}\right)-0.5772\right] \tag{9}

where \overline{\psi_s} and \psi_a are respectively the water potentials at the root surface and in the bulk soil, and a is the root radius. The value 0.5772 is the so called Euler-Mascheroni constant.

A similar formulation based on diffusivity and for water content differences is given by van Lier et al. (2006):

\overline{\theta_s}-\theta_a=\frac{x^2T_p}{4z D_w}\left[ln\left(\frac{4D_wt}{a^2}\right)-0.5772\right] \tag{10}

The first terms on the right hand side of the equations Equation 9 and Equation 10 contain different formulations of the dimensionless change of water content per day, which are, however, equivalent and can be transformed into terms which have the units of d^{-1}.

\frac{{{I}_{w}}}{\pi }=\frac{{{x}^{2}}{{T}_{p}}}{z} \\ \\ \frac{{{I}_{w}}}{\pi {{x}^{2}}}=\frac{{{T}_{p}}}{z} \\

The solution of the equation Equation 10 needs values for the soil water diffusivity and the unsaturated hydraulic conductivity as a function of volumetric soil water content or matrix potential, respectively.

1.4 Steady state and steady rate solutions for finite soil volumes

The above equation assumes that each root absorbs water from an infinite amount of soil, and that there is always more soil, at a steadily increasing distance, which can be drawn upon. In real root-soil systems this is not correct, since roots compete, and each root can be assumed to have a limited amount of soil for its sole exploitation. The root surface water potential will therefore fall faster than predicted above, and the mean water content and potential of the soil volume will decrease steadily. The simplest way of dealing with this situation is to assign an equal ‘equivalent cylinder’ of soil to each root, and to assume that all water entering the root originates in this cylinder. This is nearly equivalent to saying that roots are organized in equally spaced parallel array. There are no analytical solutions of the diffusion equation sufficiently simple to be of general use for this case (see Youngs and Gardner (1963)), and numerical solutions or approximations have to be employed. The latter are most useful for discussing the general features of the system. Passioura and Cowan (1968) have discussed the two usual methods, the ‘steady state’ and ‘steady rate’ approximations.

1.4.1 Steady state solution:

For a steady state, it is assumed that all the water is injected into the system at the periphery of the equivalent cylinder, and the appropriate equation is for the steady state flow in a cylinder:

\psi_x-\psi_a=\frac{I_w}{2\pi K}ln\left( \frac{x}{a}\right) \tag{11}

or with water content differences and diffusivity:

\theta_x-\theta_a=\frac{I_w}{2\pi D_w}ln\left( \frac{x}{a}\right) \tag{12} where x is the radius of the equivalent cylinder. Water is in fact being removed from the equivalent cylinder all the time, and this is allowed for by calculating the water loss after some convenient time interval at

\frac{\Delta \theta}{\Delta t}=\frac{I_w}{\pi (x^2-a^2)}\approx \frac{I_w}{\pi x^2} \tag{13}

and recalculating with the new values of \psi_b K_u or D_w to allow for the change in \theta.

The steady rate approximation (Cowan (1965)) is rather more realistic, since it embodies the assumption that water is being supplied uniformly from within the soil cylinder to maintain the constant flow (figure 1). The equation for this system is

{2\pi r{D}_{w}}\left( \theta \right)\frac{ \theta }{ r}=\overline I_w \frac{\left(x^2-r^2\right)}{\left(x^2-a^2\right)} this yields:

\theta_b-\theta_a=\frac{I_w}{2\pi D_w}\left[ \frac{x^2}{x^2-a^2}ln\left(ln\frac{x}{a} \right)-\frac{1}{2}\right] \tag{14}

The solutions are therefore closely similar, but (for constant \theta_b) but the steady rate equation predicts a value of \theta_a (for constant \theta_b) which is identical to that of e with the steady state equation at r = 1.65a.The equations are most useful in the relation to the mean water content \overline\theta, rather than \theta_b and then give (Passioura and Cowan (1968))

Steady state and with approximation for for x >> a, which is valid for most cases:

\overline{\theta}=\theta_a+\frac{I_w}{2\pi D_w}\left(\frac{x^2}{x^2-a^2}ln\frac{x}{a}-ln\frac{1}{2} \right) \approx \theta_a+\frac{I_w}{2\pi D_w}ln\frac{x}{1.65a} \tag{15}

A similar equation can be defined by using unsaturated soil hydraulic conductivity and soil matrix potential:

\overline{\psi}=\psi_a+\frac{I_w}{2\pi k_u}\left(\frac{x^2}{x^2-a^2}ln\frac{x}{a}-ln\frac{1}{2} \right) \approx \psi_a+\frac{I_w}{2\pi k_u}ln\frac{x}{1.65a} \tag{16}

Rearranging for I_w gives:

I_w=2\pi k_u \cdot \frac{\psi_a-\overline{\psi}}{ln\frac{x}{1.65a}} \tag{17}

1.4.2 Steady rate solution

The steady rate approximation (Cowan (1965)) is rather more realistic, since it embodies the assumption that water is being supplied uniformly from within the soil cylinder to maintain the constant flow. The Steady rate solution and its approximation for for x >> a is (Tinker (1976)):

\overline{\theta}=\theta_a+\frac{I_w}{2\pi D_w\left( x^2-a^2\right)}\left[ \left(\frac{x^4}{x^2-a^2}ln\frac{x}{a}\right)-\left(\frac{3x^2-a^2}{4}\right) \right]\\ \approx \theta_a+\frac{I_w}{2\pi D_w}ln\frac{x}{2.1a} \tag{18}

If one rearranges equation Equation 18 for I_w and assumes that \theta_a is equal to \theta_{PWP} (water content at the permanent wilting point) one gets the maximum possible water uptake rate under the given conditions Iw_{max}:

Iw_{max}=2\pi {{D}_{w}}\cdot \frac{{{\left( \overline{\theta }-{{\theta }_{PWP}} \right)}}}{ln\left( \frac{x}{2.1a} \right)} \tag{19}

1.5 Water content distribution around roots

The water content \theta_r in any distance r from the root center then is given by:

\theta_{r}=\theta_a\left( 1+{\frac{I_w}{2\pi Dw}}ln\left( \frac{r}{a} \right) \right) \tag{20}

The accuracy of the above shown approaches is, however, limited, as the assumption of a constant soil water diffusivity with respect to the radial distance to the root surface is obviously not correct. As the soil water content decreases near the root surface, the soil water diffusivity decreases as well, leading to a self reinforcing effect.

1.5.1 Boundary conditions

1.5.1.1 Inner boundary condition

For t > 0 and at the root surface, i.e. the inner boundary condition, the water influx rate per unis root length I_w is derived from the daily crop transpiration rate \overline{T_w} [cm^.d^{-1}] and the total root length L_r [cm^.cm^{-2}]. The most simple assumption is that the water uptake rate is equal for the whole root system:

t>0, \; r=a, \; 2\pi rD_w \left( \theta \right) \frac{ \theta }{ r} = \overline{I_w} \frac{\left( x^2-r^2 \right)}{\left( x^2-a^2 \right)} = \frac{\overline{Tr}}{2\pi a L_r} = -\frac{\overline{I_w}}{2\pi a} \tag{21}

A bit more realistic is the assumption of a diurnal variation of transpiration rate and water influx rate following the radiation intensity:

I_w\left( t \right)=\frac{Tr\left( t \right)}{L_r}=\frac{{{\overline{{{T}_{r}}}}\cdot}{{\frac{5}{3}}\cdot}\left[ 0.5+sin\left( {{\pi }^{.}}\frac{t}{12}+1.5 \right) \right]}{{{L}_{r}}} \tag{22}

It was further assumed that the water potential at the root surface can not be lower than the permanent wilting point, defined as a water potential of 1.5 MPa.

1.5.1.2 Outer boundary condition

It was assumed that there is a watershed between two neighboring roots. At this line the flow rate is zero:

t>0, \; r=x \;\;\;2\cdot{\pi }\cdot x \cdot D_w \cdot \frac{\theta }{ r}=0 \tag{23}

1.6 Simple simulation with steady rate approximation

The following code chunk reads in soil parameters for three different soils and sets up a simulation frame for water uptake simulation. The potential transpiration rate is set to 4 mm/day and a sinusoidal pattern is assumed over 10 days with a daylength of 16 hours. Three soil textures are selected, i.e. coarse sand (gS), sandy loam (Sl3) and silt loam (Ut2). The soil parameters are read in from a file and the remaining parameters are calculated. The parameters for the three soils are selected and a simulation frame is set up for water uptake simulation.

# a set of Van Genuchten Parameters for 3 soils, Parameters from "Rote Reihe"

# three example soils are selected. i.e. coarse sand (gS), sandy loam (Sl3) and silt loam (Ut2)
soils <- c("gS", "Sl3", "Ut2")

#sl3.soil <- list(v=0.5, ksat=89.9, alpha=0.047598, n=1.22044, theta_sat=0.363891, theta_res=0, theta_fk=0, theta_pwp=0)
#fSms.soil <- list(v=-0.328, ksat=507.5, alpha=0.15041, n=1.33576, theta_sat=0.4095, theta_res=0, theta_fk=0, #theta_pwp=0)
#ut2.soil <- list(v=0.5, ksat=29.3, alpha=0.008499, n=1.225240, theta_sat=0.399765, theta_res=0, theta_fk=0, theta_pwp=0)
#a.soil <- list(v=0.5, ksat=89.9, alpha=0.047598, n=1.22044, theta_sat=0.363891, theta_res=0, theta_fk=0, theta_pwp=0)

# read in the soil parameters from the collection "Rote Reihe" from a file
Params <- read.table("GenuchtenPars_roteReihe.csv", header = TRUE, comment.char = "[", sep = ";", dec = ".")
Params <- Params %>% dplyr::rename(Texture = Bodenart, theta_res=theta_r, theta_sat=theta_s)

# calculate the remaining parameters
Params$m <- 1-1/Params$n
Params$FK <- swc (psi=-10^1.8, alpha=Params$alpha,
                  n=Params$n,m=Params$m,theta_sat=Params$theta_s,
                  theta_res=Params$theta_res, psi_s = -1/alpha, lambda = m * n, type_swc = "VanGenuchten")
Params$PWP <- swc (psi=-10^4.2, alpha=Params$alpha,
                  n=Params$n,m=Params$m,theta_sat=Params$theta_s,
                  theta_res=Params$theta_res, psi_s = -1/alpha, lambda = m * n, type_swc = "VanGenuchten")
Params$nFK <- Params$FK-Params$PWP

# select the parameters for the three soils
df.soil.par <- Params[Params$Texture %in% soils,]

# set up a simulation data frame for water uptake simulation 

sinusf <- function (Stunde=12)
  
{
  output <-  pmax(0,1.64221194*(0.5+sin(pi*((Stunde+18)/12))))
  return(output)
}

SimTime <- 10 # days
dt_hour <- 0.5 # time step in hours


# set up a time sequence of 10 days with a time step of 0.5 hours
time <- seq(0, SimTime, 10/(24*2))

# set potential transpiration rate to 4 mm/day
Tp <- 4

# create a data frame with the time sequence and the potential transpiration rate
df.sim<- expand.grid(day = seq(from=0, to=SimTime, by = 1),
                     hour = seq(from=0, to=23.5, by=dt_hour))
dt <- dt_hour/24
df.sim <- df.sim %>% arrange(day, hour)
df.sim$time <- df.sim$day+df.sim$hour/24
df.sim$Tp <- Tp*sinusf( Stunde = df.sim$hour)

For an average daily transpiration rate of 4 mm/day and a sinusoidal pattern assuming a daylength of 16 hours over 10 days, the potential transpiration rate is shown in the following figure (Figure 1). It is important to notice that the maximum transpiration rate around noon is 9.85 mm/day, which is 2.5 times higher than the average transpiration rate.

Figure 1: Transpiration scenario with a sinusoidal pattern over 10 days with a potential transpiration rate of 4 mm/day.

1.7 Rooting scenario

In order to calculate a realistic water uptake scenario, the root length density for a certain soil layer has to be specified. In the following section we try to derive this scenario from some simple assumptions. a) The root length density is decreasing exponentially with depth. b) The rooting depth is increasing linearly with accumulated temperature sum.

1.7.1 Rooting depth

Rooting depth z_r (cm) is often found to increase linearly with accumulated temperature sum within certain development stages, but lag phases in rooting depth increase (Thorup-Kristensen (1998); Thorup-Kristensen and Nielsen (1998)) as well as diminishing rooting depth increases during maturity (Jaafar et al. (1993)) have been observed. Rooting depth increase [cm \cdot d^{-1}] therefore simply is:

\frac{d z_r}{dt}=b_{zr}\cdot T_{eff} {eq-dzrdt}

Where b_{zr} [cm\cdot°C^{-1}\cdot d^{-1}] is a constant. The base temperature for root growth was set to 0°C.

1.7.2 Vertical root distribution

Root length density, however, is not uniform over the soil profile. Typically it is decreasing exponentially from the highest values near the soil surface with depth (Barraclough, 1984; Gerwitz and Page, 1974; Greenwood et al., 1982):

L_{rv}=L_{rv0}\cdot {{e}^{-{{k}_{r}}\cdot z}} \tag{24}

where the constant k_r (cm^{-1}) is the fractional decrease in RLD per unit increase of soil depth and L_{rv0} is the root length density at zero soil depth. Integration of z=0 to a depth z=z_r were the root length density is very low yields the root length RL (cm \cdot cm^{-2}) (see Kage et al. (2000)):

L_r=\int\limits_{z=0}^{z=z_r}{L_{rv0}\cdot {{e}^{-{{k}_{r}}\cdot z}}\cdot dz}=\frac{L_{rv0}}{{{k}_{r}}}\cdot \left( 1-{{e}^{-{{k}_{r}}\cdot {{z}_{r}}}} \right) \tag{25} To calculate the average rooting density \overline{L_{rv}} (cm \cdot cm^{-3}) within a certain soil layer located between two soil depths z1 and z2 the above equation may be set up for both depths. The difference between RL at z2 and z1 divided by the distance z2-z1 gives the desired value of RLD_{av}:

\overline{L_{rv}}=\frac{L_{rv0}\cdot \left( {{e}^{-{{k}_{r}}\cdot {{z}_{1}}}}-{{e}^{-{{k}_{r}}\cdot {{z}_{2}}}} \right)}{{{k}_{r}}\cdot \left( {{z}_{2}}-{{z}_{1}} \right)} \tag{26} At the moment k_r and L_{rv0} remain unknown parameters.

for the depth z=z_r and introducing a new parameter, r_{RLD}, describing the ratio of RLD at z_r, RLD_{zr}, and RLD0 the following identity can be found for the parameter k_r:

{{k}_{r}}=-\frac{\ln \left( \frac{RL{{D}_{zr}}}{RL{{D}_{0}}} \right)}{{{z}_{r}}}=-\frac{\ln \left( {{r}_{RLD}} \right)}{{{z}_{r}}} (#eq:kr)

Thereby, the introduction of r_{RLD} avoids the necessity of using an iterative solution of the above set of equations.

Knowing this value the root length density at the soil surface can be calculated:

RLD_{0}=RL\cdot k_{r}\frac{1}{1-e^{-{k_r}\cdot z_r}} \tag{27}

For late sown winter wheat on a loess loam soil (Kage (2001)) estimated values for L_{rv0} of about 4 cm \cdot cm^{-3} and K_r of about 0.03.

For this depth we have an average root length density (L_{rv}) of 0.5 [cm \cdot cm^{-3}] and root length for water uptake (L_r) of 9.9 [cm \cdot cm^{-2}].

1.8 Root system as layers of single root cylinders

A horizontally homogeneous soil system with a growing crop and its root system can be simplified as consisting of a number of layers in which roots are regularly distributed in each layer in a hexagonal pattern and each hexagon can be simplified described by a so called single root cylinder. The radius of the single root cylinder is equivalent to the average half distance between roots x.

Cross section of Single root cylinder. x is the mean half distance between two roots, i.e. the radius of the single root cylinder and a is the root radius.

The root radius (a) is variable within a root system and differences between species are obvious. A value of 0.02 [cm] is assumed for this study, grain legumes for instance have higher root radii.

The half average distance between roots thereby can be calculated from the root length density L_{rv}: x=\frac{1}{\sqrt{{{\pi } \cdot}{{L}_{rv}}}} \tag{28}

Equation Equation 28 then gives the radius of the single root cylinder (x).

For the parameters assumed this gives an average half distance x of 0.8 cm.

Inversely, the root length density can be also calculated from the half distance between roots:

L_{rv}=\frac{1}{\pi \cdot x^2} \tag{29}

If we now relate the transpiration rate T_p to the root length within the layer of active water uptake we get the water uptake per unit root length (I_w [cm^3 \cdot cm^{-1}\cdot d^{-1}]). For our values of T_p and L_r this then is 0.0402 [cm^3 \cdot cm \cdot d^{-1}]. The average daily water flux density per unit root surface q is calculated from the root radius:

q=\frac{I_w}{2\cdot \pi \cdot a} \tag{30} which gives for our parameters a value of 0.3201 [cm^3\cdot cm^{-2} \cdot d^{-1}].

# The maximum transpiration rate at around noon [mm/day]
Trans <- max(df.sim$Tp)

# The maximum water uptake rate per unit root length [cm3/cm/d]
Iw_max <- Trans/(Lrv_av*20)*0.1

# three scenarios for the soil water potential [cm]
psi <- -c(200, 500, 2000)

# the half average distance between two roots, i.e. the radius of the equivalent cylinder [cm]
xl <- abstand_func(Lrv_av)

# a sequence of radii from the root surface to the edge of the equivalent cylinder
rads <- f_rad_arr(zylinderradius = xl, wurzelradius = a,  n = 20)

# a data frame consisting of the radii, the soil water potential, the soil texture and the soil parameters
df.swcr <- expand.grid(r=rads, psi=psi, Texture=soils)
df.swcr <- merge(df.swcr, df.soil.par, "Texture")
df.swcr <- df.swcr %>% group_by(Texture) %>% mutate(theta = swc(psi=psi, alpha=alpha, n=n,m=m, 
                                                      theta_sat = theta_sat, theta_res = theta_res))
df.swcr <- df.swcr %>% group_by(Texture) %>% mutate(Dw=diffusivity(psi=psi, v = l, ksat = Ks, alpha = alpha, n=n, m = m, theta_sat = theta_sat, theta_res = theta_res))

# set the influx rate to the maximum value
df.swcr$Iw <- Iw_max
#theta <- c(0.15, 0.175,0.2)
# calculate the water content at the root surface
df.swcr$theta_a <- pmax(df.swcr$PWP, baf(b=df.swcr$theta, Iw = Iw_max, Dw = df.swcr$Dw, xl = xl, a = a))

# calculate the soil matrix potential at the root surface
df.swcr$psi_a <- psi_b_f(b = df.swcr$theta_a, b_rest = df.swcr$theta_res, b_sat = df.swcr$theta_sat, n_par = df.swcr$n, alpha = df.swcr$alpha)

# calculate the water influx rate which may be limited by the maximum water uptake rate
df.swcr$Iw <- ifelse(df.swcr$theta_a<=df.swcr$PWP, iwmax(b = df.swcr$theta, bmin = df.swcr$PWP, Dw = df.swcr$Dw, xl = xl, a = a), Iw_max)

# calculate the water content distribution around the roots
df.swcr$theta_r <- bnf(ba = df.swcr$theta_a, Iw = df.swcr$Iw, Dw = df.swcr$Dw, Radius = df.swcr$r, a = a)
# 
df.swcr$psi <- as.factor(df.swcr$psi)

Given the assumed parameters for the root length density, the root radius, the soil water potential, and the soil texture, the water content distribution around the roots at different average soil water potentials and for three different soil textures is shown in the following figure (Figure 2). This calculation is based on the steady rate model (Equation 18 and Equation 20).

Figure 2: Water content distribution around roots at different soil water potentials and for three different soil textures calculated with the steady rate model. Upper dotted field apacity, lower dotted line: permanent wilting point.

1.8.1 Definition of matric flux potential

van Lier et al. (2006) with correction by Schröder et al. (2007) derived an expression to describe M as a function of the radial distance from the root surface r [cm] at the onset of the falling rate phase, based on the continuity equation for one‐dimensional axisymmetric flow Matric flux potential (M) \left[ {cm^{2} d^{−1}} \right] is defined as the integral of unsaturated hydraulic conductivity, K(h) [cm^.d^{−1}], over pressure head, or equivalently, as the integral of diffusivity, D(θ) \left[cm^2d^{−1}\right], over volumetric water content (\theta) \left[cm^{3}cm^{−3}\right]. For convenience, the permanent wilting point in terms of pressure head (h_w) \left[cm\right] or water content \theta{_w} \left[cm^{3}cm^{−3}\right] is chosen as the lower bound of the integral:

M=\int_{-h_w}^{h}K(h)dh=\int_{\theta{_w}}^{\theta}D(\theta)d(\theta) \tag{31}

The Matric flux potential can be calculated numerically by integrating k_u or D_w from the permanent wilting point. Figure 3 gives the curves for the matrix flux potential for three soil textures.

# chunk for the calculation of the matrix flux potential
# log scale for the soil water tension at permanent wilting point and field capacity
fromval <- log10(15000)
toval <- log10(0.1)

# steps and interval for the calculation of the matrix flux potential
steps <- 100
interval <- (fromval-toval)/steps 

# define a data frame for the calculation of the matrix flux potential
logpsivals <- seq(from=fromval, to = toval, by=-interval)
df.soil <- expand.grid(logpsi=logpsivals, soil=soils)
# column for back transformed matric potentials as absolut (positive) values
df.soil$psivals <- -10^logpsivals
# the changes of the matric potentials between two rows of the df
diffs <- diff(df.soil$psivals)

# merging with the soil parameters
df.soil <- merge(x = df.soil, y = df.soil.par, by.x="soil", by.y = "Texture" )

# add the soil water content values and the van Genuchten parameters
df.soil <- df.soil %>% group_by(soil) %>% mutate(swc =  swc(psi = psivals, alpha = alpha, n = n, m = 1-1/n,
                                                            theta_sat = theta_sat, theta_res = theta_res, type_swc = "VanGenuchten"),
                           kuvals = khy(psi=psivals, v = l, ksat = Ks, alpha = alpha, n=n,
              m=1-1/n, theta_sat = theta_sat, theta_res = theta_res),
              Dw = diffusivity(psi=psivals , v = l, ksat = Ks, alpha = alpha, n=n,
              m=1-1/n, theta_sat = theta_sat, theta_res = theta_res))

# calculation of summed ku values for the whole range of psi values
df.soil$sumku <- NA

for (soil in soils) {
#  soil <- soils[2]
  temp <- df.soil[as.character(df.soil$soil)==soil,]
  temp$sumku[1] <- temp$kuvals[1]
# integration with trapez-approach
  for (i in 2:length(temp$psivals)){
    temp$sumku[i] <- temp$sumku[i-1]+(temp$kuvals[i-1]+temp$kuvals[i])/2*diffs[i-1]   }

  df.soil[df.soil$soil==soil,"sumku"] <- temp$sumku
}

Figure 3: Matrix flux potential vs. soil water tension for three different soil texture types

1.9 Single root model with matrix flux potential

A solutions of eq. Equation 6 using the matrix flux potential approach have been derived by de_Willigen and Van_Noordwijk (1987) and were more recently used to derive sink term formulation of macroscopic soil water transport models by van Lier et al. (2006), Schröder et al. (2007) and de Jong van Lier, Q. de et al. (2008).

From mass conservation, and ignoring the root volume it follows that:

\frac{d \overline{\theta}}{dt}=-\frac{T_p}{z} \tag{32}

where \overline{\theta} is the average water content in the entire rhizospere and z [cm] is the rooting depth. Assuming that \frac{d\theta}{dt} is independent of r Eq. Equation 32 can be generalized:

\frac{\partial \theta}{\partial t}=-\frac{T_p}{z} \tag{33}

The right and left hand hand side can be expressed in dimensions of [d^{-1}] and it should be noted that T_p is given from the water influx rate times the root length per unit area. Assume a transpiration rate of 4 [mm^.d^{-1}] or 0.4 [cm^.d^{-1}], a root length density of average 0.5 [cm^.cm^{-2}] and the a rooting depth of 100 [cm], the total root length then is 50 [cm.cm^{-2}] and the water influx rate is 0.008 [cm\cdot d^{-1} \cdot cm \cdot cm^{-2}] which can be shortened to units [d^{-1}].

According to van Lier et al. (2006) eq. 38 a second order differential equation follows:

-\frac{T_p}{z}=\frac{I_wL}{z}=-\frac{-q}{r}-\frac{\partial q}{\partial r}=\frac{\partial M}{r\partial r}+\frac{{{\partial }^{2}}M}{\partial {{r}^{2}}} \tag{34}

For this equation the following solution is found (van Lier et al. (2006) eq. 39):

M=\frac{-T_p}{4z}r^2+C_1ln(r)+C_2 \tag{35}

where C_1 and C_2 are integration constants. This solution is known as the steady-rate solution. We assume the continuum approach required to apply We assume the continuum approach required to apply Darcy’s equation to hold for the geometry and the scales considered. van Lier et al. (2006) and Schröder et al. (2007) used the boundary condition

M=0;\;\;\;\;r=a \tag{36}

where a (cm) is the root radius, to express the conditions at the onset of the falling rate phase, i.e. when transport of water towards the roots is limiting water uptake rate.

To be representative for the entire constant rate phase, here we use:

M=M_0;\;\;\;\;\;\;r=a \tag{37}

where M_0\;[cm^2d^{−1}] is the matric flux potential at the root surface. The second inner boundary condition remains:

\frac{dM}{dr}=T_p\frac{x^2}{2za};\;\;\;\;\;r=a \tag{38}

where x (cm) is the rhizosphere radius, which is equal to the half mean distance between roots and is related to the root length density L_{rv}\;[cm\;cm^{−3}] by Eq.x.

Applying these boundary conditions (Eq. [11] and [12]) yields (Schröder et al. (2007), de Jong van Lier, Q. de et al. (2008)):

C_1=\frac{T_p}{2z}(x^2+a^2) \tag{39}

and de Jong van Lier, Q. de et al. (2008):

C_2=\frac{T_p}{2z}\left[\frac{a^2}{2}-(x^2+a^2)\:ln(a)\right]+M_0 \tag{40}

which combines to (de Jong van Lier, Q. de et al. (2008))

M-M_0=\frac{T_p}{2z}\left[\frac{a^2-r^2}{2}+\left(x^2-a^2\right)\:ln\left(\frac{r}{a}\right)\right] \tag{41}

Metselaar and de Jong van Lier, Quirijn (2007) showed by numerical analysis that M(r) under limiting hydraulic conditions has the same shape as under nonlimiting conditions and may be described with an expression equivalent to Eq. ??, with T_p replaced by the actual transpiration rate T_a and M_0 equal to the matric flux potential at the permanent wilting point (M_w), which by definition is equal to zero under limiting soil water conditions (Eq. ??). Therefore, in the falling rate phase,

M=\frac{T_a}{2z}\left[\frac{a_z^2-r^2}{2}+\left(x_z^2+a_z^2\right)\:ln\left(\frac{r}{a}\right)\right] \tag{42}

replacing T_p/z with I_w/(\pi x^2) gives an equation for the distribution of matric flux potential around the roots based on influx per unit root length:

M=\frac{I_w}{2\pi x^2}\left[\frac{a_z^2-r^2}{2}+\left(x_z^2+a_z^2\right)\:ln\left(\frac{r}{a}\right)\right] \tag{43}

To calculate water uptake per layer, de Jong van Lier, Q. de et al. (2008) substituted the T_p/z term by the root water uptake per unit of thickness for layer z, S_z (cm^3 cm^{−3} d^{−1}):

M_z-M_{0,z}=\frac{S_z}{2}\left[\frac{a_z^2-r^2}{2}+\left(x_z^2-a_z^2\right)\:ln\left(\frac{r}{a_z}\right)\right] \tag{44}

in which the index z refers to layer‐dependent parameter.

At a radial distance from the root surface \overline{r} (cm), the water content is equal to the mean (bulk) soil water content in the rhizosphere, and \overline{M} [m^2\:d^{−1}] is the corresponding matric flux potential. A coefficient d_z can be defined as

d_z=\frac{\overline{r_z}}{x_z} \tag{45}

Using analytical methods, the value of d was shown to be e^{−1/2} (≈0.607) assuming a constant diffusivity by van Lier et al. (2006) For soils with van Genuchten hydraulic functions, the cited authors showed with numerical simulations that the value of d is slightly lower (0.56 ± 0.06, with a median value of 0.53). Substituting \overline{M} and \overline{r} into Eq. [18] and incorporating Eq. [19] yields

\overline{M_z}-M_{0,z}=\frac{S_z}{2}\left[\frac{a_z^2-d_z^2x_z}{2}+\left(x_z^2+a_z^2\right)\:ln\left(\frac{d_zx_z}{a_z}\right)\right] \tag{46}

It should be noted that \overline{M_z} is the matric flux potential calculated from the average soil water tension or the average volumetric soil water content of each layer, which can be straightforward calculated by (numerical) integration of K_u(\psi) or D_w(\theta) for each soil layer.

The above equation again could be formulated to use the water influx per unit root length and setting M_{0,z} to zero for limiting soil water conditions. This then estimates the maximum water uptake per unit root length for the given soil conditions:

I_{wmax}=\frac{4\overline{M_z} \pi x^2}{a_z^2-d_z^2x_z+2\left(x_z^2+a_z^2\right)\:ln\left(\frac{d_zx_z}{a_z}\right)} \tag{47}

the value of d_z estimated to be 0.56.

2 Simulations

In order to analyze the validity of different model formulation for calculating the water uptake by roots, we simulate water uptake/transpiration, the water content distribution around roots for different soil textures, different soil water potentials and different computational approaches.

• Q-Model (steady rate) for uniform lrv distribution

• MFP-Model for uniform Lrv distribution

• Numerical solution for uniform Lrv distribution

• Q-Model for varying Lrv distribution

• MFP-Model for varying Lrv distribution

The Q-model is the simple steady rate solution without adressing a variation of root length density within a depth increment of the soil profile (Equation 19 and Equation 20). The is MFP-Model steady state model using a solution for the matrix flux potential with uniform root length density (Equation Equation 47 and ). The numerical solution is the solution for the matrix flux potential with varying root length density. The Q-Model for varying root length density is the steady rate solution with varying root length density. The MFP-Model for varying root length density is the solution for the matrix flux potential with varying root length density.

A simplifying assumption thereby is that the water uptake is either limited by the transpiration demand (see @) or is limited by the maximum soil water uptake rate, which is calculated from the maximum water uptake rate per unit root length (I_{wmax}) and the root length per root length in a soil layer (L_{r}). product of root length density, the depth of the layer and the maximum water uptake rate is the maximum transpiration rate.

T_{max} = L_{rv} \cdot \left( z_2-z_1\right) \cdot I_{wmax} \tag{48}

The actual transpiration rate then is the minimum of the potential transpiration and the transport limited water uptake rate:

T_{act} = min \left(T_{pot}, T_{max} \right) \tag{49}

sw^{t+1} = sw^{t}- I_w \cdot L_{rv} \cdot \Delta z \cdot \Delta t \tag{50}

The change of the soil water content

\overline \theta^{t+1} = \frac {sw^{t+1}}{\Delta z} \tag{51}

The numerical solution for the soil water transport towards roots is described in (Kage and Ehlers, 1996). It is based on a diffusion based approach and uses a fully implicit integration scheme for the solution of Equation 6.

2.1 Szenario parameters

The soil type is set to sandy loam, the root radius is 0.02 cm, the rooting depth is 20 cm and the root length density is 0.2 cm/cm^3. The transpiration rate is set to 0.4 cm/d. The start soil water potential is set to -100 hPa. The simulation time is 10 days with a time step of 0.1 hours. The sinus function for the transpiration rate is set to either FALSE or TRUE.

2.2 Comparison of steady rate, MFP and numerical model

Figure Figure 4 shows the soil water content distribution around roots for the different model approaches. The right header of each graph gives the time of the simulation in hours. It can bee seen that the drop in water content near the root surface is much steeper for the numerical solution and the MFP solution compared to the Q-model. The agreement between the numerical and the MFP solution is very good.

Figure 4: Simulation of the soil water content distribution around roots using 3 different model approaches. The header of each graph gives the time of the simulation in hours

Figure Figure 5 shows the average soil water content in the root zone over time for the different model approaches in an extraction scenario. The simulated extraction is faster i.e. less limited by water transport towards the roots for the Q-model compared to the MFP-model and the numerical solution. The MFP-model and the numerical solution show a similar behavior. The sinoidal scenario leads on average to a slower soil water extraction compared to the uniform extraction scenario.

Figure 5: Average soil water content in the root zone for the different model approaches in an extraction scenario

Figure Figure 6) shows the water uptake rate calculated with different model approaches in an extraction scenario. The Q-model shows a higher water uptake rate compared to the MFP-model and the numerical solution. The MFP-model and the numerical solution show a similar behavior. The sinoidal scenario leads on average to a slower soil water extraction compared to the uniform extraction scenario.

Figure 6: Water uptake rate calculated with different model approaches in an extraction scenario

In general the simulations demonstrate a very good approximation of the MFP-model to the numerical solution. The Q-model shows a higher water uptake rate compared to the MFP-model and the numerical solution and thereby overestimates the soil water availability for the root system. The sinoidal scenario leads on average to a slower soil water extraction compared to the uniform extraction scenario, because peak water influx rates are about 2.5 times higher than for the uniform water extraction scenario.

2.3 Effects of non-uniform root distribution

The next step is to analyze the effects of a non-uniform root distribution on the water uptake rate and the soil water content distribution around roots. The variation of the root length density is modeled with a log-normal distribution function. ### Calc variation parameters

The standard deviation of a log-normal distribution \sigma_l is calculated with the following formula:

\sigma_l = \sqrt{\ln\left(\frac{\sigma_{L_{rv}^2}}{L_{rv}^2}+1\right)} \tag{52}

The mean of the log-normal distribution is calculated with the following formula:

\mu = \ln(Lrv)-0.5\cdot\sigma^2 \tag{53}

The mean root length density is set to 0.1 cm/cm3 and the coefficient of variation is set to 200%.

From this distribution function 10 cumulative probability values (5, 15, … 95%) are choosen. For these 10 “moments” of the log-normal distribution function root length densities are calculated. The values are adjusted in a way that the total root length is identical to the mean value.

########## Calc  variation parameters ##### 

sd <- VK_LRV*Lrv/100

SA_ln <- sqrt(log((sd*sd)/(Lrv*Lrv)+1))
mean_ln <- log(Lrv)-0.5*SA_ln*SA_ln
# 10 cumulative probability values
cpvals <- seq(0.05, 0.95,0.1)

# calculate 10 "moments" of the log-normal distribution function for root length density
rld_moments <- qlnorm(p = cpvals, meanlog = mean_ln, sdlog = SA_ln)
# adjust values in a way that the total root length is identical to the mean value
rld_moments <- rld_moments*Lrv/sum(rld_moments)*10
#sum(rld_moments)/10

For the 10 classes the calculated mean values of root length density are shown in the table ?@tbl-RLDmoments and also as mid values in Figure 7. It becomes clear that the majority of the classes has a rooting density much lower than the mean value, on the other hand the classes with the highest root length density have a much higher value than the mean value. The log-normal distribution function is shown in Figure 7. The dashed lines show the borders of the 10 classes.

p	Lrv
0.05	0.006
0.15	0.014
0.25	0.022
0.35	0.032
0.45	0.044
0.55	0.061
0.65	0.085
0.75	0.122
0.85	0.194
0.95	0.419

cpvals <- seq(0.05, 0.95,0.1)
df.pmids <- data.frame(cpvals, rld_moments)
df.pmids$pval <- dlnorm(x = df.pmids$rld_moments, meanlog = mean_ln, sdlog = SA_ln)
tmp <- df.pmids
tmp$pval <- 0
df.pmids <- rbind(df.pmids, tmp)

rld <- seq(0,1, 0.0001)
prld <- dlnorm(x = rld, meanlog = mean_ln, sdlog = SA_ln)
df.log <- data.frame(rld, prld)

class_borders <- seq(0.1, 0.9, 0.1)
rld_borders <- (qlnorm(p = class_borders, meanlog = mean_ln, sdlog = SA_ln))
df.borders <- data.frame(class_borders, rld_borders)
df.borders$pval <- dlnorm(x = df.borders$rld_borders, meanlog = mean_ln, sdlog = SA_ln)
tmp <- df.borders
tmp$pval <- 0
df.borders <- rbind(df.borders, tmp)

Figure 7: Lognormal distribution of root length density divided into 10 classes of 10% cumulative probability, each.

### create grid for vectors of values as function of distance to root surface 

df.rad<- expand.grid(day  = seq(from=0, to=SimTime, by = 1),
                     hour = seq(from=0, to=23.5, by=dt_hour),
                     rad  = rads)
#summary(df.rad)

#dt <- dt/24
df.rad <- df.rad %>% arrange(day, hour)
df.rad$time <- df.rad$day+df.rad$hour/24
df.rad$psi <- 0
df.rad$swc <- 0


# water influx per unit root length [cm/d]
#df.sim$sinflow <- Vectorize (water_flow_func) (stunde = df.sim$hour, avg_transpi_rate=Avg_Transpi*10, L=L, sinus_f=sinus_func)*86400
#max(df.sim$sinflow)*L/1e7
#[cm3.cm-1.d-1]
#df.sim$q <- S_z

#plot(df.sim$time, df.sim$sinflow)

## choose soil type

#soil <- fSms.soil
#soil <- ut2.soil
wg_start <- swc(psi = psi_start, alpha = soil$alpha, n = soil$n, m = 1-1/soil$n, 
                theta_sat = soil$theta_sat, theta_res = soil$theta_res) 

#### setting up a data frame for MFP #### 

## select soil type
#soil <- fSms.soil
#soil <- ut2.soil
soil <- sl3.soil

###### create matrix flux potential data frame for choosen soil type ####
df.soil <- set_MFP_df(soil)

theta_fk <- swc(psi=-63, alpha=soil$alpha,
                n=soil$n,m=1-1/soil$n,theta_sat=soil$theta_sat,
                theta_res=soil$theta_res, type_swc = "VanGenuchten")

theta_pwp <- swc(psi=-10^4.2, alpha=soil$alpha,
                 n=soil$n,m=1-1/soil$n,theta_sat=soil$theta_sat,
                 theta_res=soil$theta_res, type_swc = "VanGenuchten")

##### set initial soil water tension and soil water contents ####
df.sim$psi <- psi_start
df.sim$swc <- swc(psi=psi_start, alpha=soil$alpha,
                         n=soil$n,m=1-1/soil$n,theta_sat=soil$theta_sat,
                         theta_res=soil$theta_res, type_swc = "VanGenuchten")

df.sim$sw <- df.sim$swc*Tiefe


if(sinus_func){
  df.sim$act_pTrans <- Avg_Transpi*sinusf(df.sim$hour)} else {df.sim$act_pTrans<- Avg_Transpi} 

# influx rate [cm3/cm/d] from pTrans [cm3/cm2/d] and root length [cm/cm2]
df.sim$sinflow <- df.sim$act_pTrans/(Lrv*Tiefe)  


#df.sim$sinflow <- 0
df.sim$iwmax <- 0
df.sim$Dw <- 0
df.sim$Iw <- 0
df.sim$Transpi <- 0
df.sim$swc_root <- 0
df.sim$MFP_nostress <- 0
df.sim$MFP_ <- 0
df.sim$MFP0 <- 0
df.sim$psi_root <- 0 

# copy values to data frame for matrix potential calculations
df.simMFP <- df.sim
df.radMFP <- df.rad

# temporary arrays for suction and water content profiles
psi_arr <- rads
scw_arr <- rads

#### Q-Model calculation for uniform lrv distribution ####
rads <- as.numeric(rads)
for (i in 1:(length(df.sim$time)-1))
     {
 # i <- 1
  df.sim[i, "psi"] <- -psi_b_f(df.sim[i,"swc"], b_rest=soil$theta_res, b_sat=soil$theta_sat,
                               n_par=soil$n, alpha=soil$alpha)
  df.sim[i, "Dw"] <- diffusivity(psi=df.sim[i,"psi"], v=soil$v, ksat=soil$ks,alpha=soil$alpha, n=soil$n,
                                 m=1-1/soil$n, theta_sat=soil$theta_sat, theta_res=soil$theta_res)
  df.sim[i, "iwmax"] <- iwmax(b = df.sim[i,"swc"], bmin = theta_pwp, Dw = df.sim[i,"Dw"], xl = xl,a = r_root  )
  df.sim[i, "Iw"] <- min(df.sim[i,"sinflow"], df.sim[i,"iwmax"])
  df.sim[i, "Transpi"] <- df.sim[i, "Iw"]*Lrv*Tiefe*10
  df.sim[i, "swc_root"] <- max(theta_pwp, baf(b = df.sim[i, "swc"], Iw = df.sim[i, "Iw"], Dw = df.sim[i, "Dw"], xl = xl, a = r_root))
  swc_arr <- Vectorize(bnf)(ba=df.sim[i, "swc_root"], Iw=df.sim[i, "Iw"], Dw=df.sim[i, "Dw"], Radius=rads, a=r_root)
  psi_arr <- -Vectorize(psi_b_f)(swc_arr, b_rest=soil$theta_res, b_sat=soil$theta_sat,n_par=soil$n, alpha=soil$alpha)
  df.rad[df.rad$time==df.sim[i,"time"], "swc"] <- swc_arr
  df.rad[df.rad$time==df.sim[i,"time"], "psi"] <- psi_arr
  
  df.sim[i+1, "sw"] <- df.sim[i, "sw"]-df.sim[i,"Iw"]*Lrv*Tiefe*dt
  df.sim[i+1, "swc"] <- df.sim[i+1, "sw"]/Tiefe
  df.sim[i, "MFP_"] <- NA
  df.sim[i, "MFP0"] <- NA
  df.sim[i, "psi_root"] <- psi_b_f(swc_arr[1], b_rest=soil$theta_res, b_sat=soil$theta_sat,
                                   n_par=soil$n, alpha=soil$alpha)
}  
df.sim$model <- "Q"
df.rad$model <- "Q"

df.simq <- df.sim
df.radq <- df.rad

For the model approaches considering the variation of the root length density the same approach as for the uniform distribution is used. The only difference is that the water balance is calculated for each class of the log-normal distribution function. The water balance is calculated according to Equation 50.

For the calculation of the soil water budget in each class of root length density it is necessary to calculate the water influx rate in each class. We choose a simple approach, avoiding iteration over the root classes. The water influx rate is calculated as the fraction of the maximum specific water influx rate in each class to the sum of the influx rates of all classes:

Iw_i =min \left( T_p \cdot \frac {Iw_{max,i} \cdot L_{rv,i} \cdot \Delta z} {\sum (Iw_{max,i} \cdot L_{rv,i} \cdot \Delta z)} \cdot \frac {10} {L_{rv,i}\cdot \Delta z}, Iw_{max,i} \right) \tag{54}

where Iw_i is the water influx rate in class.

# expand data grid for simulation with variation of root length density "df.simv"
df.simv<- expand.grid(lrv = rld_moments,
                      day = seq(from=0, to=SimTime, by = 1),
                      hour = seq(from=0, to=23.5, by=dt_hour)
)
df.simv <- df.simv %>% arrange(lrv, day, hour)
dt <- dt_hour/24 # time step back to  days
df.simv$time <- df.simv$day+df.simv$hour/24

# array of time steps
tsteps <- unique(df.simv$time)

# root length for each lrv-class in cm.m-2
df.simv$Wl <-  df.simv$lrv*Tiefe*1e4/10 # Wurzellänge in cm.m-2

# root length for each lrv-class in cm./ha
df.simv$L <- df.simv$Wl*1e4           # Wurzellänge in cm/ha

# single root cylinder radius
df.simv$xl  <- abstand_func(df.simv$lrv) # cm
df.simv$Zylindervolumen <- Volumen_func(df.simv$xl, r_root)

# average water influx per unit root length [cm/d]
#df.simv$sinflow <- Vectorize (water_flow_func) (stunde = df.simv$hour, avg_transpi_rate=Avg_Transpi*10, L=L, sinus_f=sinus_func)*86400
if(sinus_func){
  df.simv$act_pTrans <- Avg_Transpi*sinusf(df.simv$hour)} else {df.simv$act_pTrans<- Avg_Transpi} 

# influx rate [cm3/cm/d] from pTrans [cm3/cm2/d] and root length [cm/cm2]
df.simv$sinflow <- df.simv$act_pTrans/(Lrv*Tiefe)  

df.simv$psi <- psi_start
df.simv$swc <- swc(psi=psi_start, alpha=soil$alpha,
                   n=soil$n,m=1-1/soil$n,theta_sat=soil$theta_sat,
                   theta_res=soil$theta_res, type_swc = "VanGenuchten")

df.simv$sw <- df.simv$swc*Tiefe
#df.sim$sinflow <- 0
df.simv$iwmax <- NA
df.simv$Dw <- NA
df.simv$Iw <- NA
df.simv$swc_root <- NA
df.simv$MFP_nostress <- NA
df.simv$MFP_ <- NA
df.simv$MFP0 <- NA
df.simv$psi_root <- NA 
df.simv$ku <- NA
df.simv$Transpi <- NA

# copy values to data frame for matrix potential calculations
df.simvMFP <- df.simv

#df.simv$ftime <- as.factor(df.simv$time)
#tsteps <- unique(df.simv$ftime)

#### Q model with variation #####

for (i in 1:(length(tsteps)-1))
{
  # i <- 1
  # create temporary df for one time step
  tmp <-df.simv[df.simv$time==tsteps[i],]   
  
  tmp$psi <- -Vectorize (psi_b_f)(tmp$swc, b_rest=soil$theta_res, b_sat=soil$theta_sat,
                                  n_par=soil$n, alpha=soil$alpha)
  tmp$Dw <- Vectorize (diffusivity)(psi=tmp$psi, v=soil$v, ksat=soil$ks,alpha=soil$alpha, n=soil$n,
                                    m=1-1/soil$n, theta_sat=soil$theta_sat, theta_res=soil$theta_res)
  # calculate the maximum water uptake rate per unit root length [cm³/cm/d]
  tmp$iwmax <- Vectorize (iwmax)(b = tmp$swc, bmin = theta_pwp, Dw = tmp$Dw, xl = tmp$xl, a = r_root  )
  #  tmp$Iw <- pmin(tmp$sinflow, tmp$iwmax)
  #  tmp$Iw <-  tmp$sinflow*tmp$iwmax * tmp$lrv*sum(tmp$iwmax * tmp$lrv)
  
  # distribute influx according to ratio of maximum influx to sum of influxes
  tmp$Iw <-  pmin(tmp$act_pTrans*(tmp$iwmax*tmp$lrv*Tiefe) /sum(tmp$iwmax*tmp$lrv*Tiefe)*length(tmp$lrv)/(tmp$lrv*Tiefe), tmp$iwmax)
  #  tmp$Iw  <- pmin(tmp$sinflow, tmp$iwmax )
  # soil water content at root surface  
  tmp$Transpi <- tmp$Iw*tmp$lrv*Tiefe
  tmp$swc_root <- pmax(theta_pwp, Vectorize(baf)(b = tmp$swc, Iw = tmp$Iw, Dw = tmp$Dw, xl = tmp$xl, a = r_root))
  tmp$psi_root <- pmin(10^4.2, -psi_b_f(b = tmp$swc_r, b_rest = soil$theta_res, b_sat = soil$theta_sat,
                                        n_par = soil$n, alpha = soil$alpha))
  df.simv[df.simv$time==tsteps[i],] <- tmp  
  df.simv[df.simv$time==tsteps[i+1],"sw"] <- tmp$sw-tmp$Iw*dt*Tiefe*tmp$lrv
  #  df.simv[df.simv$time==tsteps[i+1],"swc"] <- df.simv[df.simv$time==tsteps[i+1],"sw"]/(tmp$Zylindervolumen*Tiefe)
  df.simv[df.simv$time==tsteps[i+1],"swc"] <- df.simv[df.simv$time==tsteps[i+1],"sw"]/(Tiefe)
}  
#df.simv <- df.simv[-nrow(df.simv),]

df.simv$model <- "QVar"

## MFP model with variation


for (i in 1:(length(tsteps)-1))
{
  #  i <- 2
  tmp <-df.simvMFP[df.simv$time==tsteps[i],]   
  
  tmp$psi <- -Vectorize (psi_b_f)(tmp$swc, b_rest=soil$theta_res, b_sat=soil$theta_sat,
                                  n_par=soil$n, alpha=soil$alpha)
  tmp$Dw <- Vectorize (diffusivity)(psi=tmp$psi, v=soil$v, ksat=soil$ks,alpha=soil$alpha, n=soil$n,
                                    m=1-1/soil$n, theta_sat=soil$theta_sat, theta_res=soil$theta_res)
  tmp$MFP_ <- MFP_f(tmp$psi, df.soil)
  
  #maximum water uptake per unit root length for the given soil conditions
  #  tmp$iwmax <- (4*pi*xl^2*tmp$MFP_)/((r_root^2-0.56^2*tmp$xl^2)+2*(tmp$xl^2+r_root^2)*log((0.56*tmp$xl)/r_root))
  #tmp$iwmax <- (4*pi*xl^2*tmp$MFP_)/((r_root^2-0.56^2*xl^2)+2*(xl^2+r_root^2)*log((0.56*xl)/r_root))
  tmp$iwmax <- MFP_Iwmax_f(MFP_ = tmp$MFP_, xl = tmp$xl, r_root = r_root)
  
  tmp$Iw <-  pmin(tmp$act_pTrans*(tmp$iwmax*tmp$lrv*Tiefe) /sum(tmp$iwmax*tmp$lrv*Tiefe)*length(tmp$lrv)/(tmp$lrv*Tiefe), tmp$iwmax)
  tmp$Transpi <- tmp$Iw*tmp$lrv*Tiefe
  
  
  #matric flux potential an the root surface
  tmp$MFP0 <- MFP0_f(MFP_ = tmp$MFP_, xl = tmp$xl, r_root = r_root, Iw = tmp$Iw)
  #    tmp$MFP0 <- pmax(0, tmp$MFP_- tmp$Iw/(2*pi*tmp$xl^2)*((r_root^2-(0.56*tmp$xl)^2)/2+(tmp$xl^2+r_root^2)*log((0.56*tmp$xl/r_root))))
  
  #  tmp$psi_root <- pmax(-10^4.2, Vectorize(psi_f)(MFP = tmp$MFP0, df = df.soil))
  tmp <-  tmp %>% group_by (lrv) %>% mutate(psi_root= psi_f( MFP0, df.soil))
  #  tmp$psi_root
  tmp <-  tmp %>% group_by (lrv) %>% mutate(swc_root = swc(psi=psi_root, alpha=soil$alpha,
                                                           n=soil$n,m=1-1/soil$n,theta_sat=soil$theta_sat,
                                                           theta_res=soil$theta_res, type_swc = "VanGenuchten"))
  # update water content for next time step
  
  df.simvMFP[df.simvMFP$time==tsteps[i],] <- tmp  
  
  df.simvMFP[df.simvMFP$time==tsteps[i+1],"sw"] <- tmp$sw-tmp$Iw*tmp$lrv*Tiefe*dt
  df.simvMFP[df.simv$time==tsteps[i+1],"swc"] <- df.simvMFP[df.simv$time==tsteps[i+1],"sw"]/(Tiefe)
  
  #  df.simvMFP[df.simvMFP$time==tsteps[i+1],"swc"] <- df.simvMFP[df.simvMFP$time==tsteps[i+1],"sw"]/(tmp$Zylindervolumen*Tiefe)
}  

df.simvMFP$model <- "MFPVar"

The simulations with the different models are compared in the following plots (Figure 8). The first plot shows the extraction rates of the different models. The second plot shows the soil water contents at the root surface and at the soil matrix.

From the comparison of the simulated water extraction rates it becomes clear, that the approach considering a variation of root length density leads to much lower water extraction rates. The Q-Model with uniform root length density distribution on the other hand has the highest extraction rates.

Figure 8: Comparison of extraction rates for Q-Model and MFP-Model with uniform and variable root length distribution

Figure 9: Comparison of soil water contents for Q-Model and MFP-Model with uniform and variable root length distribution

The different root length density leads to a non uniform extraction of soil water from the different classes. Wheras in the class with the highest root length density the soil water content is reduced to the permanent wilting point, the soil water content in the class with the lowest root length density is only slightly reduced (Figure 9).

Figure 10

The extraction rates of the different models are shown in the following plot (Figure 11). The red line the average extraction rate of the root length density classes calculated according to the MFP model, and the blue lines are the extraction rates for the different root length density classes calculated according to the MFP model.

Figure 11: Comparison of extraction rates for MFP-Model with root length distribution

p <- ggplot(data=df.simv, aes(x=time, y=psi, color=lrv))
p <- p + theme_bw(base_size = 16)
p <- p + geom_point()
p <- p + ylab("Matric potential [cm]")
#p <- p + scale_y_continuous(limits=c(-1,1))
p

Figure 12: Comparison of matric potentials at the root surface for the differenct root length density classes

p <- ggplot(data=df.simv)
p <- p + theme_bw(base_size = 16)
p <- p + geom_point(aes(x=time,y=Iw/iwmax, color=lrv))
p

Warning: Removed 10 rows containing missing values or values outside the scale range
(`geom_point()`).

Figure 13: Comparison of the ratio of actual water extraction to maximum water extraction rates

References

Barraclough, P.B., 1984. The growth and activity of winter wheat roots in the field: Root growth of high-yielding crops in relation to shoot growth. Journal of agricultural Science (Camb.) 103, 439–442.

Carslaw, H.W., J.C. Jaeger, 1946. Conduction of heat in solids. Clarendon Press, Oxford.

Cowan, I.R., 1965. Transport of water in the soil-plant atmosphere system. J. appl. Ecol. 2, 221–239.

de Jong van Lier, Q. de, van Dam, J.C., Metselaar, K., Jong, R.D., Duijnisveld, W.H.M., 2008. Macroscopic root water uptake distribution using a matric flux potential approach. Vadose Zone Journal 7, 1065–1078. https://doi.org/10.2136/vzj2007.0083

de_Willigen, P., Van_Noordwijk, M., 1987. Roots, plant production and nutrient use efficiency (Ph.D. thesis). Agricultural University Wageningen, The Netherlands.

Gardner, W.R., 1960. Dynamic aspects of water availability to plants. Soil Sci. 89, 63–73.

Gerwitz, A., Page, E.R., 1974. An empirical mathematical model to describe plant root systems. J.appl.Ecol. 11, 773–781.

Greenwood, D.J., Gerwitz, A., Stone, D.A., Barnes, A., 1982. Root development of vegetable crops. Plant and Soil 68, 75–96.

Jaafar, M.N., Stone, L.R., Goodrum, D.E., 1993. Rooting depth and dry matter development of sunflower. Agronomy Journal 85, 281–286.

Kage, H., 2001. Improving nitrogen use efficiency of intensive cropping systems (Habilitationsschrift). Hannover, University of; University of Hannover, Hannover.

Kage, H., Ehlers, W., 1996. Does transport of water to roots limit water uptake of field crops? Zeitschrift für Pflanzenernährung und Bodenkunde 159, 583–590. https://doi.org/10.1002/jpln.1996.3581590609

Kage, H., Kochler, M., Stutzel, H., 2000. Root growth of cauliflower (brassica oleracea l. Botrytis) under unstressed conditions: Measurement and modelling. Plant and Soil 223, 131–145.

Metselaar, K., de Jong van Lier, Quirijn, 2007. The shape of the transpiration reduction function under plant water stress. Vadose Zone J 6, 124–139. https://doi.org/10.2136/vzj2006.0086

Passioura, J.B., Cowan, I.R., 1968. ON SOLVING NON-LINEAR DIFFUSION EQUATION FOR RADIAL FLOW OF WATER TO ROOTS. Agricultural Meteorology 5, 129+. https://doi.org/10.1016/0002-1571(68)90044-7

Philip, J.R., 1957. The physical principles of soil water movement during the irrigation cycle, in: ICID (Ed.), Proc. Int. Congr. Comm. Irrig. Drain., Proc. Int. Congr. Comm. Irrig. ICID, New Dehli, pp. 125–154.

Schröder, T., Javaux, M., Vanderborght, J., Vereecken, H., 2007. Comment on “root water extraction and limiting soil hydraulic conditions estimated by numerical simulation.” Vadose Zone Journal 6, 524–526. https://doi.org/10.2136/vzj2007.0036

Thorup-Kristensen, K., 1998. Root growth of green pea (pisum sativum l.) genotypes. Crop Science 38, 1445–1451.

Thorup-Kristensen, K., Nielsen, N.E., 1998. Modelling and measuring the effect of nitrogen catch crops on the nitrogen supply for succeeding crops. Plant and Soil 203, 79–89.

Tinker, P.B., 1976. Roots and water. Phil. Trans R. Soc. Lond. B 273, 445–461.

van Lier, Q. de J., Metselaar, K., van Dam, J.C., 2006. Root water extraction and limiting soil hydraulic conditions estimated by numerical simulation. Vadose Zone Journal 5, 1264–1277. https://doi.org/10.2136/vzj2006.0056

Youngs, E.G., Gardner, W.R., 1963. A problem of diffusion in the infinite hollow cylinder. Soil Science Society of America Journal 27, 475–476. https://doi.org/10.2136/sssaj1963.03615995002700040033x