Spatial Heterogeneity Enables Higher Root Water Uptake in Dry Soil but Protracts Water Stress After Transpiration Decline: A Numerical Study

A common assumption in models of water flow from soil to root is that the soil can be described in terms of its representative or effective behavior. Microscale heterogeneity and structure are thereby replaced by effective descriptions, and their role in flow processes at the root‐soil interface is neglected. Here the aim was to explore whether a detailed characterization of the microscale heterogeneity at the scale of a single root impacts the relation between flow rate and pressure gradient. Numerical simulations of water flow toward a root surface were carried out in a two‐dimensional domain with a randomized configuration of spatially variable unsaturated hydraulic conductivities and varying boundary conditions, that is, increasing and decreasing root water uptake rates. By employing Matheron's method, the soil hydraulic properties were varied, while the effective hydraulic conductivity (corresponding to the geometric mean) remained unchanged. Results show that domains with a uniform conductivity could not capture important features of water flow and pressure distribution in spatially variable domains. Specifically, increasing heterogeneity at the root‐soil interface allowed to sustain higher root water uptake rates but caused a slower recovery in xylem suction after transpiration ceased. The significance of this is that, under critical conditions, when pressure gradients and flow rates are high, microscale heterogeneity may become an important determinant and should not be neglected in adequate descriptions of water flow from soil to root in dry soil.


Introduction
Spatial variability, on all scales of observation, is the most pertinent feature of water flow and solute transport in porous media. A sizeable number of studies therefore focused on developing different approaches to characterize the spatial variability of the saturated and unsaturated zones in order to improve predictions of water flow and solute transport for the sustainable management and the protection of subsurface water resources (Gerke, 2006;Gharari & Razavi, 2018;Jarvis et al., 2016;Kuhlmann et al., 2012;Renard & de Marsily, 1997;Sanchez-Vila et al., 2006;Šimůnek et al., 2003;Vereecken et al., 2007;Zhang, 2010). It is well known that heterogeneities in the unsaturated zone can induce flow instabilities, which lead to preferential and nonequilibrium flow during infiltration, and drastically decrease the mean residence time of water and solutes, such as manure, fertilizers, or pesticides, in the vadose zone, thereby threatening groundwater resources (Allaire et al., 2009;Jarvis et al., 2016;Šimůnek et al., 2003;Vogel et al., 2010). This key challenge in vadose zone research has led to a great emphasis in addressing the problem of wetting of spatially variable, structured soils. The problem of drying in heterogeneous and structured soils, a crucial process in the exchange of water between soil, vegetation, and atmosphere, however, has received less attention.
Two basic ecosystem processes, which are an integral part of the hydrological cycle, primarily drive the drying of soil. Evaporation from the soil surface and, driven by transpiration demand, root water uptake. In contrast to wetting, which often involves a wide range of pores and where flow instabilities are caused by heterogeneities such as fissures, cracks, or root channels (Gerke, 2006), water flow in dry soils usually only takes place in the fine pores of the soil matrix. Water movement is therefore generally assumed to be largely uniform and to occur, for a given boundary condition, while pressure head and water content are locally in equilibrium. On these assumptions, the soil can be described as a homogeneous domain with a representative or effective behavior. Spatial heterogeneity of soil hydraulic properties within the matrix Figure 1. Illustration of the emergent relationship between absolute xylem pressure and increasing and decreasing rates during a diurnal cycle of root water uptake (after Carminati et al., 2017). When root water uptake decreases after state B (orange line), xylem pressure during recovery (blue line) remains higher than the expected trajectory (dotted line), hence inducing an apparent hysteresis in the relation (stage C). For a more detailed explanation, see text.
are thereby commonly replaced by effective descriptions of the soil (Gerke, 2006;Vogel & Roth, 2003). The reason often is that models need to be easy to handle, or the lack of information or computational resources make it impossible to explicitly consider the heterogeneity in hydraulic properties (Zhang, 2010). This has made effective parameters efficient in a wide range of studies of water movement in dry soil across scales. For many of these applications it is reasonable to assume that the water flow is uniform, because the flow rate of water in dry soil is comparably low. However, it is not clear whether this assumption remains true for all soil hydraulic conditions during drying and whether heterogeneity within the soil matrix can affect this assumption under certain conditions.
A typical example to explore such conditions is water uptake by plant roots. Root water uptake occurs when, due to water loss by transpiration, the reduction in xylem pressure and the ensuing tension are propagated to the root following the cohesion-tension principle (Steudle, 2001). When the xylem pressure is more negative than in the surrounding soil, water flows from soil into the root along the paths with the lowest resistance.
Figure 1 provides a stylized example of the relationship between root water uptake (L/T) and absolute xylem pressure (L) for a diurnal cycle of increasing and decreasing root water uptake. At the beginning of the diurnal cycle (stage A) the xylem pressure increases linearly with root water uptake. The increase during stage A is proportional to the resistance of the plant (Deery et al., 2013;Passioura, 1980). As the rate of root water uptake further rises (stage B), significant pressure gradients develop in the soil, while xylem pressure increases nonlinearly in order to sustain root water uptake. In stage B xylem pressure and the rate of root water uptake are largely dependent on the hydraulic condition and properties of the soil. When after stage B root water uptake rates decrease (because stomata close) absolute xylem pressure remains markedly above the anticipated recovery curve (dotted line; Carminati et al., 2017;Deery et al., 2013;Passioura, 1980). The result is an apparent hysteresis in the relationship (stage C). The pressure difference between the expected and the observed trajectory is denoted by ΔP, which corresponds to the pressure dissipation in the soil. Carminati et al. (2017) showed that by using measured effective descriptions of the soil hydraulic properties, a single root model of water uptake could only provide an adequate account of the observed pattern during the phase of increasing rates of water uptake but failed to reproduce the hysteretic pattern observed during declining rates of water uptake. Several possible explanations exist for the hysteresis in the relation between root water uptake and xylem pressure. These include, among others, the buildup of a significant osmotic pressure along the root surface, due to the accumulation of salts; an alteration of the hydraulic behavior of the rhizosphere during drying; a local depletion of water along the entire root system in the regions where root water uptake is relatively high; and an increase of the flow resistance at the root-soil interface, due to root shrinkage and loss of contact. But these tentative explanations alone may not be enough to describe the observed hysteresis to its full extent.
So far, however, the usage of effective descriptions of the soil in models of plant-water relations has not been challenged. Most studies of water flow from soil to root rely on uniform descriptions of the soil properties around roots, thereby neglecting the role soil heterogeneity and structure may have on water flow from soil to root. It is therefore not known whether effective descriptions undermine the attempts to properly characterize the hysteresis in the relation between root water uptake and xylem pressure.
Within the sizeable body of literature of saturated and unsaturated water movement, it is well established that spatial heterogeneity within a flow domain can produce hysteretic flow patterns. (Davies & Beven, 2015;Gharari & Razavi, 2018;Mantoglou & Gelhar, 1987b). By investigating three-dimensional transient variably saturated flow, Binley et al. (1989) aimed to define effective hydraulic properties that could reproduce different time-dependent deterministic events but did not find a single effective parameter for low-permeability soils. Mantoglou and Gelhar (1987a) applied their stochastic theory of three-dimensional transient unsaturated flow (Mantoglou & Gelhar, 1987b) to evaluate effective hydraulic conductivities in stratified soils. They found that the effective hydraulic conductivities showed a significant hysteresis, which resulted from the spatial variability of the hydraulic soil properties. Davies and Beven (2015) use the multiple interacting pathways model to show that the storage-discharge relationship is dependent of the antecedent condition, because of storage-related differences between flow velocities and celerities.
Most hydrological models investigating the role of spatial heterogeneity on water flow place great focus on the pedon up to catchment scales, while only little work has focused on the effects of microscale heterogeneity . In particular, little attention was dedicated to how small-scale heterogeneity (scale < 1 cm) may affect soil-plant water relations. The objective of this study was to explore whether a more detailed account of the spatial variability of soil hydraulic properties at the single root scale could improve modeling the plant physiological response to varying soil hydraulic conditions-specifically the observed delay in the recovery of the xylem pressure after stomatal closure. It is hypothesized that at this scale the role heterogeneity becomes important when significant pressure gradients develop around the root (stage B, Figure 1). At this critical point, microscale heterogeneity in hydraulic properties causes water to flow along nonuniform flow patterns, leading to disequilibrium conditions and a slow recovery of xylem pressure after stomatal closure (stage C, Figure 1).
In order to test this hypothesis, numerical simulations of water flow toward a root surface were carried out in a spatially heterogeneous two-dimensional flow field. The numerical experiments were repeated for different degrees of spatial variability of the hydraulic conductivity function. The effective conductivity, however, remained unchanged during all simulations by applying Matheron's conjecture (Matheron, 1967) for effective conductivities in two-dimensional flow fields (Zhang, 2010). The rationale for choosing a two-dimensional numerical approach was based on the assumption that essential features of water flow in three-dimensional coordinates could be described in 2-D with a moderate model complexity (Dagan, 1989;Roth, 1995).

Water Movement in Dry Soil and Water Uptake by Plant Roots
Root water uptake is driven by pressure gradients along the pathway of water flow from the soil toward the root surface and into the root (Carminati et al., 2016;Passioura, 1988;Raats, 2007). The pressure gradients are proportional to the flow rate and to the hydraulic conductivity of the root and the soil, which in unsaturated soil is a function of the pressure head. Early accounts of models of root water uptake (Gardner, 1960;Philip, 1957) assumed that root water uptake was driven by the pressure gradient between the bulk soil and the root surface. But it could later be shown that water flow along the root-soil interface was, under most conditions, controlled by the difference between the xylem pressure and the pressure head at the root surface, and the hydraulic conductivity along the root-soil interface (Taylor & Klepper, 1975;Newman, 1969aNewman, , 1969bMolz, 1976). Only when in dry soil the soil hydraulic conductivity is much lower than the rate of root water uptake do significant pressure gradients develop in the soil near the root. At this stage root water uptake is mainly controlled by the soil hydraulic condition (Carminati et al., 2016;Passioura, 1988).
Water movement across the root tissue from the root surface to the xylem takes place along apoplastic, symplastic, and transcellular pathways (Steudle, 2000). There is still lack of clarity which of the pathways is predominant, and available evidence indicates that the relative importance of each pathway can depend on plant species, development stage, and the hydraulic condition of the root (Kim et al., 2018;Ranathunge et al., 2017). However, root water uptake can be well approximated by a single pathway and by using an average root hydraulic conductivity (Landsberg & Fowkes, 1978;Zarebanadkouki et al., 2014). Water flow in the root can thus be described by the following equation (Carminati et al., 2016;Couvreur et al., 2014;Javaux et al., 2013) where q root (cm/s) is the root water uptake rate, K root (cm/s) is the overall hydraulic conductivity of the root tissue between the root surface and the xylem, h root (cm) is the pressure at the root surface, and h xylem (cm) is the xylem pressure. Since the fundamental work of Gardner (1960), an extensive range of studies simulated the water flow from the bulk soil to the root in terms of a cylindrical geometry. In this case the root is conceptualized as a linear sink in the center of a soil cylinder, and the outer radius of the soil cylinder is often defined as half the mean distance between roots, which is related to the root length density.
Water flow in unsaturated soils is described by the Darcy-Buckingham law. For soils with two-dimensional spatial variability and assuming no influence of gravity the equation can be modified to is the vector of the spatial coordinate, h (cm) is the pressure head and K (cm/s) is the hydraulic conductivity. Angle brackets indicate the expected value of the probability distribution function.
Combining the Darcy-Buckingham law with the mass conservation equation yields the Richardson-Richards equation (Raats & Knight, 2018), which in its h-based form (Celia et al., 1990) is expressed by where C(h) = (h) h is the soil hydraulic capacity (cm −1 ), is the water content (cm 3 /cm 3 ), and t is time (s). The Richardson-Richards equation is widely used to describe the redistribution of water and the emerging gradient in pressure along the soil-root pathway in unsaturated soils. Here an important extension for the case of water uptake by plant roots in heterogeneous soil is the introduction of the spatially variable unsaturated hydraulic conductivity K(x, h). The characterization of K(x, h) is discussed in more detail in the following section.

Characterization of Spatial Heterogeneity and the Effective Hydraulic Conductivity
In heterogeneous porous media the spatially variable hydraulic conductivity function K(x) can be considered a spatial random field with a given probability density function (PDF). The PDF of K(x) is usually assumed to follow a log-normal distribution. This is supported by many field and numerical studies where it was found that spatially variable hydraulic conductivity fields are often positively skewed. (Jury et al., 1987;Matheron, 1967;Neuman & Orr, 1993;Nielsen et al., 1973;Vereecken et al., 2007). It is therefore often preferred by many authors to use the general relation where Y (x) is the log-transform of the spatially variable hydraulic conductivity with a Gaussian PDF. For a Gaussian PDF the spatial variability of Y (x) can be entirely described by its first-and second-order statistical moments, that is, its expected value ⟨Y (x)⟩ and its variance 2 Y . Y (x) can also be further expanded into a deterministic and a random component (Reynolds' decomposition) so that where ⟨Y (x)⟩ is the expected value and Y (x) ′ is a zero mean perturbation with the variance 2 Y . In most practical applications it is not convenient to continuously define K(x) within a given flow domain because there is considerable uncertainty about the exact spatial distribution of the hydraulic properties Table 1 Mualem- van  and because the exact implementation of the spatial variability would require very large computational resources for most problems dealing with water flow in heterogeneous media (Vereecken et al., 2007;Zhang, 2010). In the past, an enormous amount of research was therefore devoted to identifying proper ways to describe heterogeneous media in terms of their effective hydraulic behavior and to replace the spatially variable parameter K(x) by a single effective parameter (e.g., effective retention or effective hydraulic conductivity). By this method, the spatially variable medium is essentially replaced by a homogeneous medium of the same size that shows the same effective hydraulic behavior as the heterogeneous example (Zhang, 2010). It is worth emphasizing at this point that the same notion of an effective behavior also forms the basis of a number of experimental approaches, ranging from the scale of a single soil core to the field scale, to determine the hydraulic properties of porous media. Therefore in experimental results as well as in models, spatial variability at scales smaller than the representative elementary volume is commonly neglected and incorporated in the medium properties Gerke, 2006).
There exist a large number of approaches to define the effective hydraulic conductivity in heterogeneous media. However, most efforts focus on saturated flow problems (Sanchez-Vila et al., 2006). One of the early findings was that, regardless of the dimensionality and spatial arrangement of K(x), the lower boundary of the effective conductivity K eff is the harmonic mean K H and the upper boundary is the arithmetic mean K A so that K H ≤ K eff ≤ K A (Renard & de Marsily, 1997).
For log-normal statistically isotropic porous media, Matheron (1967) conjectured that the effective saturated hydraulic conductivity corresponds to the harmonic mean K H in media with one-dimensional heterogeneity and to the geometric mean K G in media with two-dimensional heterogeneity. Later it was shown that this conjecture was valid under saturated conditions (Neuman & Orr, 1993). Zhang (2010) extended Matheron's method to unsaturated soils and used a numerical exercise with Miller-similar synthetic soils to show that for different spatial dimensions of heterogeneity it can well predict the effective unsaturated hydraulic conductivity in isotropic porous media. Given the relationship Y h = ln K(x, h) and an average pressure headh, Zhang (2010) described the effective unsaturated hydraulic conductivity K eff (h) by where is the variance of Y h for a givenh, and d expresses the dimension of the heterogeneity within the flow field. Note that according to equation (7) K eff (h) also approaches K H for d = 1 and becomes K G (h) for d = 2 (Zhang, 2010). The following will refer to this framework, since it always assumed here that K(x) is log-normally distributed and the case of statistical anisotropy is beyond the scope of this study.

Soil Hydraulic Properties
A sandy loam potting mix (Carminati et al., 2017) was utilized for the soil hydraulic base properties of the numerical experiments. The soil was highly aggregated at the mesoscale, had a relatively high organic matter content, and its pore range was such that it remained well aerated even at high water contents. The water retention and the hydraulic conductivity of the soil were measured with the HYPROP ® device, METER Group, Munich, Germany, following the evaporation method (Peters & Durner, 2008;Schindler, 1980). The Mualem- van Genuchten (MvG) model (van Genuchten, 1980) was inversely parametrized by fitting the Richardson-Richards equation to the measured pressure heads (supporting information Figure S1). A more detailed description of the soil hydraulic parameters is given in Text S1. The results of the parametrization are presented in Table 1.

Spatial Variability of Hydraulic Properties
In order to test how different degrees of heterogeneity affect the outcome of the numerical experiments, the MvG parameters (x) and K s (x) were assumed to be log-normal spatial random fields, while the other y of the log-normally randomized flow parameters (x) and K s (x). The parameters were randomly varied around the geometric mean K G (red line) in order to maintain the same effective conductance for different values of 2 y ( Figure S2). The arithmetic mean K A is shown by the green line. Gray lines illustrate the K(h) relationship corresponding to 10 (5th, 10th, 25th, 35th, 45th, 55th, 65th, 75th, 85th, and 95th) quantiles of the log-normal distribution of (x) and K s (x), which were employed in the random flow domain. parameters remained constant. A simple linear positive correlation was assumed between and K s ( Figure  S3). After being log-transformed to y (x) = ln (x) and y K s (x) = ln K s (x) the parameters were expanded by their zero mean perturbations y ′ (x) and y K s ′ (x), respectively (Vereecken et al., 2007). By applying the Reynolds' decomposition the parameters can be rewritten as The variance 2 y of the normally distributed zero mean perturbations y ′ and y K s ′ was changed between 0 and 6 ( Figure 2) to generate different degrees of heterogeneity of K(x, h) within the flow domain. For this modeling exercise the perturbations were simplified to 10 classes for each 2 y corresponding to the 5th, 10th, 25th, 35th, 45th, 55th, 65th, 75th, 85th, and 95th quantiles of the generated frequency distributions of y ′ and y K s ′ (see Figure 2). All calculations presented in this section were performed by using the statistical programming language R (R Core Team, 2018). The resulting frequency distributions of (x) and K s (x) were used to generate statistically isotropic random flow fields with different degrees of spatial variability, which is described in the following section.

Characteristics of the Flow Domain
A two-dimensional rectangular flow domain was used for this numerical study. The domain was set up to give a simplified representation of a longitudinal cross section of the root-soil interface. It had a height of 4 cm and a total width L dom of 5.2 cm. Along L dom the domain was subdivided into a soil domain and a root domain with distinct hydraulic properties. The root domain had a width (root radius) of 0.1cm and was located at the left boundary of the flow domain Γ left . It consisted of two 0.05 cm wide longitudinal compartments. The first compartment had a low resistance to water flow to represent the xylem vessel. The second compartment described the root tissue, separating the xylem vessel from the soil domain. Based on data provided in Carminati et al. (2017) and assuming an active root length L r,active of 10% of the total root length, an overall conductivity K root of 4.5819 × 10 −11 cm/s was calculated for the membrane compartment. The soil domain had a width of 5.1 cm and was composed of symmetric hexagonal subdomains with a width L hex of 0.1 cm. L hex here is equal to the correlation length within the flow domain. The soil hydraulic properties of the soil domain were defined using the MvG model with the parameters given in Table 1. The parameters s , r , n and were assigned equally to all hexagons. In order to allocate the spatially variable parameters (x) and K s (x), the hexagons were separated into 10 classes. The classes were randomly assigned to the hexagons and no spatial autocorrelation was assumed between the subdomains at this scale (Figure 3a). For a given spatial variance 2 y (Figure 2) each of the 10 classes represented one of the 5th, 10th, 25th, 35th, 45th, 55th, 65th, 75th, 85th, and 95th quantiles of the frequency distribution of (x) and K s (x) (see above). For 2 y = 0, which represented the case of a "homogeneous" control, (x) = and K s (x) = K s (Table 1). To test whether K eff of the generated soil domains was equivalent to K G and did not vary for different 2 y , steady state numerical simulations were carried out ( Figure S2). Results showed that the differences of K eff between the domains were only about 0.01% (Table S1) and therefore regarded as insignificant. At initial condition the ratio of hydraulic conductivities above K G (h 0 ) was 0.44 corresponding to a ratio of 0.56 of hydraulic conductivities below K G (h 0 ) (Figure 3b).
Contrary to many studies that have used various forms of the one-dimensional radial model (see section 2) it was chosen here to apply equation (1) in a two-dimensional setting and along the longitudinal direction of the root instead of the radial direction. This had several advantages. First, it remained relatively simple to conceptualize the spatial variability of the unsaturated hydraulic conductivity in the numerical approach. Second, due to the larger surface a wider range of the spatially variable unsaturated hydraulic conductivity could be obtained along the longitudinal cross section of the root as compared to the radial cross section. Third, in two dimensions the effective hydraulic conductivity K eff is directly described by the geometric mean K G (see above). Fourth, the results of the numerical simulations could be easier generalized to other two-dimensional problems of water flow in heterogeneous porous media.

Variation of the Correlation Length L hex
It was also tested how the width of the hexagonal subdomains (correlation length) L hex changed the outcome of the model experiments. L hex was both increased to 0.2 cm and reduced to 0.05 cm, while the spatial variance was kept constant at 2 y = 2. For both variations of L hex six realizations of the flow domain were generated in order to account for changes in the random configuration of the hexagonal subdomains. Inverse

10.1029/2019WR025501
simulations of steady state water flow were used to assure that K eff did not significantly change between the realizations.

Numerical Experiments
Numerical experiments of root water uptake were conducted in the flow domain by solving the Richardson-Richards equation (equation (4)) in two dimensions. A constant head (Dirichlet) boundary condition of h = −2, 000 cm was set at the right side Γ right of the flow domain. At the left domain boundary Γ left , where the xylem compartment was located, varying rates of root water uptake q root (Neumann boundary) were applied during the experiments. Both top Γ top and bottom Γ bottom boundary conditions were fixed to a zero flux boundary. The initial condition within the flow domain was h 0 = −2, 000 cm. This initial condition was chosen because a number of experiments have shown that close to this value the soil resistance becomes a limiting factor for the flow of water from the soil into the root (Passioura, 2010(Passioura, , 1980). An adaptive triangular computational grid was utilized to discretize the flow domain. The lengths of the elements were between 3.6 × 10 −4 and 1.8 × 10 −2 cm so that 144,256 elements were obtained. The flow problem was solved with the finite element method using the PDE solver provided in the Partial Differential Equation Toolbox™ in MATLAB.
The numerical experiments were repeated for different degrees of spatial variability by varying 2 y (Figure 2), while the "homogeneous" case 2 y = 0 was used as a control. An experimental cycle consisted of a stepwise increase of the root water uptake rate q root and a subsequent drop in water uptake. The objective was to mimic a diurnal increase in transpiration demand followed by stomata closure when the transpiration demand can no longer be sustained by soil water flow. At the beginning of each experiment q root was −1 × 10 −7 cm/s. Afterward q root was raised by q root = 1 × 10 −7 cm/s every 20 min. To ensure a stable solution q root was not raised abruptly, but gradually within a time step of 3 min. When the relation between q root and the pressure at the root surface h root showed first signs of nonlinearity q root was only raised until an additional pressure drop of max. Δh root = 10,000 cm was recorded. At the maximum pressure h max = h nlin +Δh root , where h nlin is the pressure at the beginning of the nonlinear deviation, root water uptake was ceased and the system was allowed to reequilibrate. The time step between two consecutive model solutions was 5 min in the linear phase and 1.5 min during the nonlinear phase.
The total pressure lost in the soil within the first 3 hr of the reequilibration phase 3h was calculated by integrating h root over time where t max (min) is the time at h max and t 3h (min) is the time 3 hr after h max was reached.

Results
In Figure 4 the results of the numerical experiments over time are shown for different 2 y . The simulation cycles started in dry soil with an initial pressure of h 0 = −2,000 cm. There were considerable differences in soil water content depending on the degree of the spatial heterogeneity. For 2 y = 0 the average volumetric water content was ⟨ v ⟩ = 0.12 and for 2 y = 6 it was ⟨ v ⟩ = 0.15. q root was increased stepwise to find at which point the transition between steady state root water uptake and nonsteady state occurred, or more precisely when the relation between q root and h root began to be nonlinear. For all 2 y this was already the case within the first hour of the numerical cycle. The rate of increase, however, was smaller as 2 y increased. When q root was below 3 × 10 −7 cm/s the xylem pressure h xylem increased faster than h root , because the resistance of the root membrane was still higher than the soil resistance. Above q root = 3 × 10 the soil resistance increased rapidly, and the increase of h root caught up with h xylem showing a highly nonlinear increase until the maximum pressure h max . At h max the trajectory was almost vertical for all 2 y and the system was quickly failing. Values of h max varied between −16,100 cm ( 2 y = 0.5) and −17,500 cm ( 2 y = 6), which is not far above the permanent wilting point, often defined at −15,000 cm. With a higher 2 y , h max was reached at bigger maximum root water uptake rates q max . This is also shown in Figure 5 where the relation between pressure at the root surface and the rate of root water uptake is compared for different 2 y . During the nonlinear phase, when the soil resistance became important, bigger 2 y also enabled higher root water uptake rates at similar pressure heads. This is mirrored by the less rapid deviation from the linear relationship between q root and h root at low root water uptake rates. When the maximum pressure h max was reached and root water uptake was halted, h root recovered markedly slower with increasing 2 y . While the first record during the recovery (5 min after reaching h max ) was approximately −6,600 cm for 2 y = 0 (homogenous flow domain), h root stayed at approximately −9,800 cm in the flow domain with 2 y = 6. The total pressure lost within the first three hours of the recovery phase was calculated by equation (8). The results are plotted in Figure 6b for varying degrees of soil heterogeneity. For low degrees of heterogeneity the pressure lost in the soil 3h after 3 h already showed a rapid increase and continuously rose with higher 2 y . Figure 7 illustrates the two-dimensional pressure profile within the soil domain at h max for different 2 y . The profiles show that at h max pressure gradients extended further into the soil with higher 2 y . While for 2 y = 0 significant gradients in pressure developed in the first 1.5 cm near the root, in the most heterogeneous case ( 2 y = 6) significant pressure gradients extended up to nearly 3 cm into the soil. Notable differences could also be observed regarding the location of the steepest pressure gradients. In the flow domain with 2 y = 0  the steepest gradient in pressure developed in the first millimeter near the root surface. With increasing 2 y the location of the steepest pressure gradient moved further away from the root surface. This was accompanied by the development of low pressure zones near the root that extended up to 0.5 cm into the soil. In increasingly heterogeneous soil domains these low pressure zones developed along highly conductive areas, while the steepest gradients occurred at discontinuities of zones with a high conductivity (Figures 7 and S4). y at the maximum pressure h max and for the maximum root water uptake rate q max , which was applied at the left-hand side of the domain. The homogeneous flow domain ( 2 y = 0) clearly shows a one-dimensional flow velocity profile. When flow properties were varied, on the other hand, two-dimensional flow profiles with preferential flow paths emerged. With an increasing 2 y these preferential flow paths became more pronounced and higher flow velocities were reached along these pathways, leading to a higher overall difference in flow velocity within the flow domain. With a large heterogeneity in K(x, h) water predominantly moved along a few preferential flow paths. In many areas of the domain flow velocities remained extremely low, even in some areas close to the  root surface, thus, contributing only little to the overall flow. With increasing 2 y a network of preferential flow paths evolved, spanning between the right and left domain boundaries.
Overall the variation of the correlation length L hex did not lead to significant differences in the model outcome ( Figure 9). The median of q max as well as of 3h were equal and larger than the homogeneous case for both L hex ( 2 y = 0, Figure 6). However, for the larger L hex the model results showed a higher variability among realizations. Figure 10 provides an exemplified spatial representation of the two-dimensional pressure head profile and flow velocities for varied L hex . For L hex = 0.05 cm zones with the steepest pressure gradient were more uniformly shifted away from the root surface, while for L hex = 0.2 cm such shifts were much more localized. This is also reflected in the spatial representation of flow velocities. In domain realizations of L hex = 0.2 cm, zones with a high flow velocities were confined to only a few, but large, preferential pathways. While the length of these high velocity zones did change for different 2 y , no such effect on the extent of the high velocity pathways was observed for varied L hex . For all realizations of L hex these pathways had a length of 3 cm (see also Figure 8).

Discussion
The findings of this study support the hypothesis that local variations in soil properties affect the emergent relation between flow rate and pressure profile in heterogeneous soil. The simulations showed that there was a clear divergence between the behavior of flow fields with an explicit account of spatial heterogeneity and the uniform control, which became more pronounced as heterogeneity increased.
The numerical simulations showed that, even when the effective conductivity across flow domains is equal in steady state, differences in the microscale heterogeneity of the hydraulic conductivity impacts the emergent relation between water flow and the pressure profile. In more heterogeneous domains higher root water uptake rates could be sustained (Figures 5 and 6a). The presence of highly conductive pathways possibly caused the notable reduction of the flow resistance. The steepest gradients emerged further away from the root surface (Figure 7), with water flowing along a network of preferential flow paths and bypassing a substantial fraction of the domain. This channeling became more pronounced as heterogeneity increased (Figure 8). The results also showed a notable overall deceleration in the reequilibration of pressure in domains with a higher spatial heterogeneity, after root water uptake was ceased (Figures 5 and 6b). The correlation length did not impact the maximum sustainable flux and the hysteresis. However, for higher correlation length, the simulations showed a greater variability among realizations.
The impact of heterogeneity was most apparent during two phases: first, when water uptake rates were high and large gradients in pressure emerged in the flow domain; and second, in the reequilibration phase after root water uptake was ceased. In the phase when the rate of root water uptake was high, a big proportion of the water flow was constrained to conductive regions near the root surface. These relatively conductive pathways led to the observed channeling within the flow domain (Figure 8), which caused the notable increase in the effective conductivity of the domain. When root water uptake rates were low or water flow was at steady state ( Figure S2), however, no such divergence was observed and the effective behavior of all flow domains was the same. The overall deceleration in the recovery of pressure, observed during the follow-up phase of no transpiration, may have been caused by a multitude of local flow imbalances and nonuniform pressure profiles. The imbalances were the result of the juxtaposition of zones with hydraulic conductivities that differed by several orders of magnitude. The presence of relatively conductive zones near the root surface induced local alleviations of the pressure gradient (Figure 7) so that the steepest pressure gradients shifted further away from the root surface and to zones with a low hydraulic conductivity ( Figure S4). As the gradients extended further from the root surface, their relaxation was slower, causing the observed recovery delay of the root pressure after water uptake was ceased. This effect was maintained also when the correlation length was reduced, although zones with steep pressure gradients were shifted away more uniformly from the root surface. A higher correlation length led to a more localized presence of relatively conductive zones. The length of these preferential flow paths increased with spatial variability, but did not depend on the correlation length. Even in domains with a low degree of spatial variability, the observed lengths of these high velocity zones well exceeded the mean radius of exploitation typically assumed in root system models (Metselaar & De Jong van Lier, 2011).
The fact that the length of the preferential flow paths did not depend on the hexagon size Figure 10 is surprising, as one could expect the preferential pathways to scale with the correlation length. A possible explanation is that the length of these pathways depends on the probability of high conductive regions to be connected, which is scale independent, and on the dissipation of pressure along these pathways. The dissipation of pressure decreases with increasing heterogeneity-i.e. the conductivity differences become more pronounced with increasing variance (Figure 8). The problem of juxtaposed conductive and nonconductive zones within a flow domain has been addressed in the percolation theory. In the percolation theory the critical (percolation) probability threshold at which a network becomes connected is scale independent (Berkowitz & Balberg, 1993). That such scale independence is similar to the outcome of the present analysis and the invariant length of the connected pathways is an intriguing hypothesis that deserves further studies.
The above considerations indicate that during specific drying conditions, which are critical for water uptake by plant roots (Carminati et al., 2016), the relationship between water flow and pressure gradient cannot be captured, when heterogeneity in soil hydraulic properties is replaced by uniform descriptions of the soil. It may be conceivable to modify the effective soil hydraulic conductivity in such a way that it captured the relative attenuation in the effective resistance caused by the channeling at high root water uptake rates. However, as the results presented here have also shown, it is not clear whether first, such an adjusted conductivity can be easily found (Binley et al., 1989;Fiori & Russo, 2007;Vogel et al., 2010); and secondly, such attempts would fail to describe the observed delay in the pressure redistribution during the recovery phase. An adjustment of the effective hydraulic conductivity alone will therefore probably not suffice to properly predict the observed behavior in spatially variable soil matrices. Models decoupling the relation between water pressure and water content Ross & Smettem, 2000) or including hysteresis in the hydraulic properties might enable to achieve a better fit. To the authors' knowledge, the problem of nonuniform flow under the critical drying conditions tested here has not been rigorously approached before. However, the literature on unsaturated nonuniform flow during wetting conditions is vast and may offer potential avenues to approach this issue Gharari & Razavi, 2018). Vogel et al. (2010), for instance, offered a suitable framework, which included 1) conceiving the soil hydraulic properties in terms of a probability space rather than a fixed functional relation, 2) the usage of stochastic moments to characterize the spatial variability and 3) collecting information on the spatial configuration and correlation lengths of the soil properties.
Such a framework could be particularly promising with regard to the characterization of the hydraulic properties of the rhizosphere, which is a key component in the water flow from soil to root. Important advances in imaging techniques make it possible to obtain high resolution images for the much necessary quantification of the physical soil structure around the roots (Carminati et al., 2016;Hinsinger et al., 2009). These images have shown that compared to the bulk soil, the rhizosphere exerts a considerable heterogeneity driven by the propagation of mucilage within the pore space, changes in soil structure around roots (e.g., by soil swelling and shrinkage), the emergence of water repellent zones, the formation of rhizosheath and physical interactions with microorganisms (Carminati et al., 2016). Spatial variability in root-soil contact (partial contact) may further add to this complexity (Carminati et al., 2013;de Willigen et al., 2018;Herkelrath et al., 1977;North & Nobel, 1997;Veen et al., 1992). Some of these features (mucilage, partial contact) have already been used in attempts to explain the hysteresis in the relation between xylem pressure and rates of root water uptake (Carminati et al., 2017). Yet it has not been thoroughly studied whether a more general understanding of how these features contribute to local heterogeneities in water flow and the delay in the redistribution of pressure can provide new insights into the puzzling hysteresis.
The results of this study are based on numerical simulations of water flow, which were carried out in a two-dimensional flow field with a random configuration of hexagonal subdomains. It was assumed that the PDF of the spatially variably hydraulic conductivity function was positively skewed and could be characterized by a log-normal distribution. A range of field and numerical studies have found positively skewed hydraulic conductivity functions, however, most analyses were carried out at a much larger scale. A more detailed characterization of the hydraulic properties at the rhizosphere scale is therefore needed to test the plausibility of this assumption, since it has been shown that conductivity fields often deviate from this ideal case (Neuweiler & Vogel, 2007). It remains unknown whether the variances chosen here properly represent the heterogeneity encountered at this root-soil interface. Moreover, in our models the conductivities were distributed randomly. However, in the rhizosphere this might not be the case. For instance, air-filled gaps with low hydraulic conductivity might develop along and around all or part of the root (Carminati et al., 2013). Despite these simplifications, our study shows how small-scale heterogeneity at the root-soil interface impacts the relation between transpiration and xylem pressure head. The consequences are manifold: First, our results challenge the validity of models of water flow in unsaturated conditions during variable boundary conditions using effective properties; second, they show how microheterogeneity at the root-soil interface impact plant scale plant water relation when the soil is dry and transpiration is high. Interestingly, heterogeneity partly explains the observed hysteresis in the relationship between transpiration and xylem pressure head, although other reasons exist. Furthermore, it is known that plants have mechanisms to respond and adapt to soil heterogeneity, for instance, by hydropatterning (Bao et al., 2014) and the growth of root hairs and mycorrhiza (Carminati et al., 2017;Keyes et al., 2013). Explicitly considering heterogeneity in the rhizosphere might help to understand soil and root water relations and plant adaptation to dry soil.