A Big Root Approximation of Site‐Scale Vegetation Water Uptake

Land surface model (LSM) predictions of soil moisture and transpiration under water‐limited conditions suffer from biases due to a lack of mechanistic process description of vegetation water uptake. Here, I derive a “big root” approach from the porous pipe equation for root water uptake and compare its predictions of soil moistures during the 2010 summer drought at the Wind River Crane site to two previously used Ohm's law analog plant hydraulic models. Due to a fuller representation of pressure gradients and flows within a complex root system architecture, the new formulation achieves somewhat improved fit and significantly lower bias compared to the Ohm's law analog models. A key advantage of the improved physical representation is the increased robustness of fits and predictions, making it less liable to overfitting. This new mechanistic model advances our understanding of vegetation water limitation at site scale with potential to improve LSM predictions of soil moisture, temperature and surface heat, water, and carbon fluxes.


Introduction
Vegetation access to soil water can play a key role in limiting terrestrial fluxes of water (Schlesinger & Jasechko, 2014), heat (Wild et al., 2015) and carbon dioxide (Fatichi et al., 2016;Nobel, 2009), yet it remains poorly understood at relevant scales. A lack of accurate descriptions of plant water limitation gives rise to significant prediction errors in terrestrial or land surface models (LSM), which aim to predict these surface fluxes as a bottom boundary condition on atmospheric circulation in earth system simulations (Dahlin et al., 2015;Dietze et al., 2011;De Kauwe et al., 2015). Results from Coupled Model Intercomparison Project Phase 5 (CMIP5) simulations of the Intergovernmental Panel on Climate Change (IPCC) show overestimates of evapotranspiration (ET) and associated overland cooling in most terrestrial regions, especially under soil moisture limited ET regimes (Mueller & Seneviratne, 2014). CMIP5 simulations underline the critical role of plant access to variable soil moisture in mediating surface-atmosphere fluxes and vegetation dynamics (Green et al., 2019). This study therefore focuses on mechanistic description of vegetation access to soil moisture as a means of improving predictions of soil moisture, vegetation water uptake, and attendant fluxes in terrestrial models.
A common approach in large-scale land surface, dynamic vegetation or ecosystem dynamics models is to calculate water uptake from each soil layer in proportion to an assumed root length density fraction reduced by a stress factor, which is often a linear function of local soil moisture status (Best et al., 2011;De Kauwe et al., 2015;Drewniak, 2019;Feddes et al., 2001;Fisher et al., 2015;Hidy et al., 2016;Kala et al., 2015;Krinner et al., 2005;Le Moigne et al., 2012;McMurtrie et al., 1990;Medvigy et al., 2009;Oleson et al., 2013). It has been shown that changing this root water uptake formulation can give rise to significant differences in predicted transpiration across models (Jing et al., 2014;Ferguson et al., 2016). As awareness of the importance of plant hydraulics for ecosystem ecology and earth system science has grown (Sack et al., 2016), Figure 1. Comparison of conceptual models: Resistance diagrams for the Ohm's law models with subsequent soil layers connected to the root base in parallel (a) and in series (c). Conceptual diagram of big root model (e) as a solution to the continuous, nonlinear porous pipe equation. Circles represent water potentials at the root base (green, top left), in the soil (rightmost column) and within the root (inner nodes in resistance networks in (a) and (c), and along the solid root water potential curve in (e)). Sparsity patterns and composition of linear systems for layer root water potential (̄x) in the Ohm's law parallel (b), series (d), and big root (f) models, color coded to correspond to elements in conceptual diagrams. Coefficient matrices show nonzero diagonals present in each model, with yellow on the main diagonal (layer being solved for), orange and red for layers above and dark and pale blue for layers below. Right-hand side vector entries are linear combinations of the soil water potentials of individual layers, identically color coded. Green denotes a term for the root base boundary condition. Colored nodes in conceptual diagrams show elements involved in the solution for root water potential in yellow layer.
it has become apparent that more mechanistic descriptions of plant hydraulics are needed in LSMs to capture vegetation behavior during drought (Sperry and Love, 2015), especially since the frequency of such water-limiting events is likely to increase in the future (Allen et al., 2010).
This has given an impetus to replace the classical plant water stress formulation with a more process-based model (Bonan et al., 2014;Kennedy et al., 2019;Xu et al., 2016), based on an Ohm's law analog for water flow in plants (Sperry et al., 2016). While such an approach works well for stems, its weakness in representing root water uptake is that it linearizes all plant resistances and places them between the root base and each soil layer either exclusively in parallel (Sperry et al., 2016) or in series (Amenu & Kumar, 2008). Each of these configurations is an end-member of the spectrum of possible root system architectures (RSA) in that all other RSAs combine resistances to within-plant flow between layers in series with resistances in parallel. The actual arrangement of resistances in RSA is crucial to controlling uptake from variably saturated soils Doussan et al., 2006;Javaux et al., 2013;Lobet et al., 2014;Sulis et al., 2019), as it determines the rate at which water uptake dissipates water potential gradients (Bouda & Saiers, 2017;Bouda et al., 2018), which drive flow along the soil-plant-atmosphere continuum (Sperry et al., 2016). As water uptake in a given layer dissipates the potential gradient between adjacent layers, it reduces how much water can be taken up in deeper layers. The strength of this effect depends on the soil water potentials in adjacent layers and the full arrangement of resistances. Whereas the parallel resistance structure (Figure 1a) allows uptake in one layer to affect all others directly via the potential at the root base, the model that places all layers in series ( Figure 1c) assumes that the plant water potential value effective in uptake and in transport are the same in a given layer, oversimplifying the effects of adjacent water potential values on each other. Neither is a suitable conceptual model for the complex cross-layer effects known to emerge in spatially explicit, 3-D RSA simulations (Bouda & Saiers, 2017).
A conceptual alternative is provided by the "porous pipe" model of root water uptake (Landsberg & Fowkes, 1978), which describes the continuous variation of water potential ( ) along the length (s) of roots taking up water with the differential equation where the left-hand side term represents the dissipation of the gradient of root water potential ( x ), which is proportional to the water potential difference between the root and the soil ( s ), with the ratio of conductance into (k r ) and along (K x ) the root as the constant of proportionality. Analytical solutions to this equation have recently been extended from the case of single roots with homogeneous properties to arbitrary RSAs (Bouda et al., 2018;. This paves the way for solutions to this continuous model that predict the root water potential effective in water uptake in each layer (̄x) as a function of varying soil water potentials and canopy demand ( Figure 1e). Here, I show that solving the porous pipe equation analytically for values of̄x, under the assumption that all roots at a site can be represented with a single "big root," yields a new model that better represents the effects of root system architecture on belowground pressure gradients than Ohm's law analogs. This big root model yields improved vegetation water uptake predictions during drought at a forested site, compared to the Ohm's law models.

Model Assumptions and Problem Setup
Terrestrial models divide the subsurface into n depth layers, which are assumed to have homogeneous hydraulic properties. Soil hydrological equations, describing the movement of water in porous media, are t a n h ( i ∕2)∕ i -Linear coefficient discretized over these layers to obtain solutions for volumetric water content and water potential s in each layer. Water uptake is calculated in line with the resulting hydration status of each soil layer and the evaporative demand given by atmospheric conditions and a model for photosynthesis and stomatal transpiration. So as to fit neatly with such a description of subsurface hydrology, let us assume that both the plant and the soil are subdivided into depth layers with homogeneous properties. Let us also assume a single value of the conductances K i x and k i r for the root, though such assumptions may be overcome in the future, as in .
The point of departure for the present model is not to linearize equation (1) a priori, but rather to solve the nonlinear problem for a mean value of water potential,̄i x , in each layer, i ∈ 1, .., n, which represents the value effecitve in root water uptake. This is the value that correctly yields water uptake Q r , when used in the equation where k r is the soil-root conductance per unit length and S is the total length of a root segment. If we assume that each layer is traversed by a single root segment, then this mean water potential is given bȳ

Assembly of System of Equations
Solutions to equation (1) can be used to evaluate the integral in this equation, yielding six possible expressions for the mean root water potential, Step 2: Substitute expressions for boundary conditions in layer i in boundary conditions at layers i ± 1.
Step 3: Replace boundary conditions at layers i ± 1 with expressions in mean root and soil water potentials in layers i ± 1 and i ± 2. The result is an equation for mean root water potential in a given layer in terms of those in the four surrounding layers and the corresponding soil water potentials.
These expressions are linear in the layer boundary conditions: soil water potential ( i s ) and the root water potential ( i ) or its gradient (G i ) at layer boundaries (subscript 0 at layer bottom and 1 at top); linear coefficients c i , ∈ 1, .., 5 are determined by root properties (k i r , K i x , S i ), as shown in Table 1. Systems of linear equations assembled from the relations (4a)-(4f) can be used to find simultaneous solutions for̄i x in all layers i ∈ 1, .., n. The assembly of such systems can be illustrated in three steps ( Figure 2). The assembly of a system using layer boundary water potential gradients (G i 0 , G i 1 ) as boundary conditions begins (Step 1) from equation (4a), where the boundary conditions are as yet unknown. To eliminate the unknown G i 1 ( Step 2), substitute from equation (4a) evaluated in layer i − 1 (note that K i from conservation of mass). The equation now contains the unknown G i−1 1 . An expression for this boundary condition can be found (Step 3) by combining equation (4b) evaluated in layer i − 2 with equation (4d) in layer i − 1. As i−1 0 = i 1 (from continuity), we can solve both (4b) and (4d) for the water potential at the boundary between the two layers and set the resulting expressions equal to each other. Rearranging yields an expression for G i−1 1 , which can then be then substituted into the overall equation. Analogous steps are used to eliminate G i 0 , only using expressions evaluated in layers i + 1 and i + 2. The resulting equation can be rearranged algebraically to the form which shows that the mean xylem water potential in a given layer is a linear function of the mean xylem and soil water potentials spanning the two surrounding layers. Expressions for coefficients and for each layer are given in Table S1 in the supporting information.
A total of 16 separate systems can be derived, by choosing different implicit boundary conditions. The complement to the system described above is that, which uses i 0 and i 1 as boundary conditions in all layers instead of gradients. To derive this, we take equation (4f) as the basis in Step 1 and again for the substitution in Step 2. In Step 3, we use equations (4b) and (4d) again, but eliminate the gradient terms to solve for values of x at the boundary between layers i − 2 and i − 1 (i + 1 and i + 2). After rearrangement, we retrieve the form of equation (5), albeit with new expressions for coefficients and (again, given in Table S1). Each of these separate parameterizations thus yields a distinct system of linear equations (with unique coefficients), each of which is derived from the underlying root parameters and all of which yield the same solutions for a given set of boundary conditions.

Domain Boundary Conditions
The root is assumed to have no flow at the lower boundary, consistent with the "porous pipe" model, simplifying the resulting equations for layers n − 1 and n. This assumption follows from the idea that the rooting depth is located somewhere in the last layer and thus no water uptake takes place deeper down.
Boundary conditions at the top of the domain may either specify the water potential ( B ) or gradient (G B ) at the stem base (or root collar). In either case, the fixed value replaces the respective expression for the boundary condition in the equations for Layers 1 and 2, simplifying the resulting coefficients. For transient simulations of root system water uptake from a numerically represented soil, these values can vary at every time step to simulate water demand from the above-ground portion of the vegetation. In the present study, the root conceptualization is considered in isolation from aboveground processes, but these boundary conditions should ultimately result from coupling with an above-ground plant and atmospheric model.

Simultaneous Solution
Equation (5) can be used to construct a system of linear equations for successive layers with boundary conditions appropriately incorporated, which can be rearranged into the form where A is the matrix of coefficients , x is the solution vector of values̄i x for i ∈ 1, .., n and b is the vector of values resulting from the second term on the right-hand side of equation (5) and boundary conditions in Layers 1 and 2. Thus, if we know the root parameters k r , K x , and S for each layer, then for any set of boundary conditions B or G B and i s in each layer, this pentadiagonal system provides the analytical, exact solution for the resultinḡi x in each layer. This is termed here the "big root" model.

Big Root Model
The difference between the present big root model and the Ohm's law analog formulations can be demonstrated on the sparsity of the coefficient matrix in equation (6), and how it arises. This sparsity pattern determines how uptake rates in each model interact across layers. It corresponds to each model's representation of the effects of RSA on root water potential gradients.
As the Ohm's law models begin by linearizing equation (1), each necessarily yields a linear system of equations for̄x, where the sparsity patterns of matrix A follows directly from the arrangement of resistances assumed, as does the position of the stem base boundary condition in the right-hand-side matrix. The parallel Ohm model assumes each layer's root water potential value depends only on the same layer's soil water potential and the root base boundary condition (Figure 1a), resulting in a simple diagonal matrix and the stem base condition affects each entry of the right hand side vector, propagating this information to each layer (Figure 1b). The series Ohm model instead connects each layer by a linear resistance to the adjacent layers (Figure 1c), resulting in a tridiagonal matrix, and restricts the stem base information to the top layer (Figure 1d). The present big root approach, by contrast, never assumes linearity. When the continuous form of equation (1) is solved analytically for the values̄x in each layer, the derivation yields a linear system with a pentadiagonal coefficient matrix and the stem base condition only affects the top two layers (Figure 1f). The entries of the matrix are related to the root parameters by exact expressions. In this case, the linear system is retrieved from their derivation, rather than assumed up front and its precise form results deductively from the problem itself.
Remarkably, a five-banded linear dependence has been found between values of̄x calculated post hoc from analytical solutions to equation (1) on explicit 3-D root networks in layered soils (Bouda & Saiers, 2017). That is, a version of equation (5) with an intercept and no constraints on or produced almost perfect fits of the results, regardless of the 3-D RSA used to generate the values. One might conjecture that equation (5) is one of a class of models, which can represent analytical solutions to equation (1) for other RSA with a pentadiagonal linear system. In these models, the specific RSA does not determine the sparsity of the matrix, as in the Ohm's law models, but rather the analytical expressions for the coefficients and . The expressions derived for the coefficients here may thus be understood as a special case of, or "big root constraints" on, this general class of models, previously called an RSA Stencil.
The denser matrix sparsity pattern of the present model represents more complex interactions across layers than in the Ohm's law models-it allows for more cross-layer effects on water uptake. Bouda and Saiers (2017) found that models with simpler sparsity patterns were insufficient to fully represent the cross-layer effects of water uptake in more complex RSAs. The greater the uptake in each layer, the more the water potential gradient inside the root system is dissipated there, leaving less "sucking force" in deeper layers. In the parallel model, this effect is assumed to be absent. In the serial model, it is overly constrained by a first-order approximation, in which a single value of water potential is assumed effective both in vertical transport through the root and soil-root transport in a layer (Bouda & Saiers, 2017). The five-diagonal scheme allows the two values of to be decoupled, an important effect of the nonlinearity of equation (1), enhanced by RSA complexity. Thus, even though the present model imposes constraints on the coefficient matrix entries that follow from the assumption of a single root, its matrix sparsity pattern may be expected to provide a better representation of how RSA affects water uptake than the corresponding (series) Ohm model. How it performs compared to the parallel Ohm model depends on two things: whether RSA can be ignored under the parallel assumption and whether constraining the RSA Stencil by single root assumptions limits its versatility too much.

Model Inversion
To employ the big root idealization for a single plant, a species, or a land cover type in practical applications, model parameters have to be estimated by inversion. To establish whether, and to what extent, inversion could yield theoretically true values, I first compared them to the exact values obtained analytically for synthetically constructed single roots. Twenty cases of a single root with randomly varying properties (k r , K x , S) across seven layers were generated as follows. Central parameter values (L = 5 cm, k r = 2.83e−8, K x = 5.03e−10) were multiplied by a factor distributed uniformly between 0.5 and 1.5, independent in each layer. Values for̄x and Q R in all layers were then found analytically for each case for all combinations of the following boundary conditions: B varied from −0.4 to −1.2 MPa and s from −0.2 to −1 MPa, each by increments of 0.2 MPa. Big root stencils were then fit to these solutions using the soil water potentials as boundary conditions, as well as either the water potential at the stem base or its gradient.
It has already been shown that 16 distinct formulae for parameters and can be obtained analytically, which yield distinct values. Each parameter set (linear system) is separately capable of producing exact values of̄i x . Preliminary work showed, however, that convergence rates of model inversions on correct parameter values are inadequate when using only a single set of parameter constraints (e.g., one or the other expressions for and in Table S1), leading to stalled optimization runs or very high numbers of iterations (results not shown). By contrast, simultaneously using complementary sets of expressions (e.g., both sets in Table S1) in the inversion leads to high rates of convergence (Table 2), from an uninformed initial guess (all values start at 1). Either of the two complementary sets of coefficients can be used in predictive runs, as they yield identical results.
The results also show that regardless of the data set employed in the inversion, the resulting linear coefficients ( , ) yield high accuracy predictions of root uptake (Q R ). Specific values of the linear coefficients ( , ) do not, however, uniquely determine the underlying root parameters (k r , K x , S), as only ratios or products of the latter occur in the expressions for the former. Inversion of the model thus cannot yield these three parameters directly.
Depending on the information used, the inversion can, however, uniquely constrain the physically meaningful quantities k r S (total layer soil-root conductance) and K x ∕S (integrated layer vertical conductance). When both̄i x and Q i R are used in each layer to constrain the model, a unique value of k r S is implied, yielding also the correct K x ∕S (Table 2, rows 1 and 2). This is true regardless of the boundary condition set at the stem base (water potential B or its gradient G B ).
As field data are almost certain not to contain target values for̄i x , a second set of inversions used soil-root flows alone as the data. In this case, with a boundary condition on B , the physically meaningful values for total soil-root and total vertical conductance are still obtained with high accuracy in all but the nth layer. In the nth layer, a different set of linear coefficients is obtained that correctly predicts flows (Q i R ). This shows that for the degenerate form of the equation that occurs in the nth layer (where G n 0 = 0), the values of k r S and K x ∕S are not uniquely determined by the flow solution alone. Instead, the inversion converges on one of a set of equivalent stencils ( , ). With a boundary condition on G B , the values of k r S and K x ∕S are only approximated, with errors on the order of 10%, though again a stencil (set of linear coefficients , ) yielding exact solutions for Q R is found.
Depending on the data used in model inversion, different interpretive value thus can (or cannot) be ascribed to the effective big root parameters underlying the resulting stencil. These results suggest that measuring stem base water potentials is a valuable component of field site design if physically meaningful big root parameters are sought. In their absence, the resulting stencil ( , ) values only yield approximations of the underlying conductances but are still constrained by the physical behavior of a single root with some real conductance values that yield matching uptake rates at the calibration points.

Wind River Crane Case Study 3.2.1. Site and Data Description
To evaluate the big root model as a tool for predicting water uptakes at field scale, its predictions were compared to soil moisture observations, along with the two Ohm's law models. Each of the three models was used to predict vegetation water uptake and soil moisture at the Wind River Crane AmeriFlux site, a well-described old growth Douglas fir (Pseudotsuga menziesii) site with deep soils of volcanic origin (Shaw et al., 2004;Warren et al., 2005;Wharton, 1998). The site has a high mean annual precipitation (>2,000 mm/yr) and also a pronounced summer drought period when less than 5% of the annual total falls. The site's valley topological position leads to a relatively shallow water table, which at times intrudes on the root zone.
Precipitation and soil moisture data at depths of 20, 30, 40, 50, 60, 100, and 150 cm, at half-hourly resolution for the year 2010 (Wharton, 1998) were used in conjunction with published soil textures to define soil characteristic curve parameters (Brooks & Corey, 1964;Clapp & Hornberger, 1978;Warren et al., 2005). The summer drought period was subselected from the data (Figure 3a) as the period with no precipitation events that noticeably affect shallow-layer soil moisture. This was done in order to minimize soil water flows in the resulting data set. During the selected drought period (Figure 3b), there is a notable diurnal pattern in the soil moisture data down to the 60 cm depth (Figure 3d, also visible in Figure S2), suggestive of a dominant influence of vegetation water uptake.

Assumptions and Scope
To construct a data set of boundary conditions from these observations that would be sufficient for the model evaluations, several assumptions had to be made. One is on the attribution of soil moisture differences between time steps to soil-root transfer versus vertical soil water movement. In the absence of any data on this, the simplest assumption was made: that all water movement during this period is due to plant water uptake or loss. This assumption is least unrealistic during the drought period, when soil moisture conductivities are lower, there are no periods of infiltration and rapid drainage, and there is a notable dominance of the diurnal transpiration-refilling signal in the soil moisture data. While rather restrictive, this assumption greatly simplifies the problem as it enables root water uptake to be calculated for each time step using centered differences of volumetric water content.
Another assumption is required on how best to determine the total transpiration rate of the vegetation. In order to use measured evapotranspiration, an assumption would have to be made on its partitioning into components. Using the resulting transpiration rate as the stem base boundary condition on water uptake would assume no time-lag between the two as well as a full correspondence between total transpired amounts and uptake in the layers measured. The first assumption is unlikely to hold in a 60-m-tall old growth Douglas fir stand, while the latter is unlikely since no soil moisture data shallower than 20 cm were available and plants may also have access to the phreatic zone. Differences between total uptake and transpiration due to time-lag as well as soil-root flows outside the observed domain would be difficult to resolve with any degree of certainty. Thus, total flow at the stem base was instead assumed to simply be the sum across the soil column of the layer uptakes, as calculated above.
These assumptions, as well as the use of soil moisture characteristic curve parameters derived from texture, clearly limit the scope of the simulation exercise. The different conceptualizations are being compared here as models of root water uptake, specifically with respect to how water potential gradients propagate on the assumed RSA, given this interpretation of the data. In this respect, their performance cannot be expected to differ based on choosing a less simple set of assumptions, while any resulting increases in overall simulation realism would remain disputable.

Inverse and Forward Simulation Methods
To ensure model comparisons under like conditions, and since big root parameters must be obtained by fitting the model, these parameters as well as Ohm's law resistances were found by model inversion. For consistency, a constraint was imposed on the total soil-root conductance in each model, so that at a water potential difference of 0.5 MPa, the uptake is 0.2 mm/hr, close to the maximum value found from the data. All inversions were performed using an interior-point algorithm (Byrd et al., 1999(Byrd et al., , 2000Waltz et al., 2006), with initial values of soil-root conductance in each layer following the relative distribution of root area (Figure 3c) reported for the site (Warren et al., 2005).
Inversions were performed on limited subsets of data, but always evaluated in comparison to data from the entire drought period. Model fit was estimated from inversions using a larger, evenly distributed calibration data set: three 72-hr periods were subselected from the drought period, at its start, middle, and end ( Figure 3b). The robustness of the fits and prediction error was then estimated by performing 1,000 separate model inversions on shorter, 72-hr subsets of the data. In this case, each such subset consisted of three separate 24-hr blocks, chosen at random from the overall data set without replacement. Calibration data included soil water potentials and calculated root water uptakes in each layer, as well as the flow boundary condition. The time course of soil volumetric water content was not included in calibration data. The resulting 1,000 separate inversions produced a distribution of prediction error for each model.
Forward simulations with each model used soil moisture in each layer at the start of the drought period as the initial condition and solved for soil water contents at each subsequent time step, with water potentials calculated from predicted soil moistures. The simulation employed a Crank-Nicolson time-stepping scheme and represented the and s relationship implicitly, using an iterative solver to account for its nonlinearity (Smith et al., 2002;Nocedal & Wright, 2006). Aside from bias, and root-mean-square error (RMSE), model fit was quantified using Nash-Sutcliffe Efficiency (NSE), the metric most commonly used to evaluate hydrological models (Ritter & Muñoz-Carpena, 2013). The metric has support from − inf to 1, with 0 indicating that the model performs no better than the mean of the data in predicting individual values and 1 representing perfect fit. Values deemed acceptable vary by author and range between 0.6 and 0.8 in the literature, though such judgments clearly depend on model application. The question can be quantitatively reformulated in terms of the number of times n t that the standard deviation of the observations exceeds RMSE, yielding one possible acceptability threshold of NSE = 0.65, that corresponds to n t ≈ 0.7 (Ritter & Muñoz-Carpena, 2013). Best practices (Ritter & Muñoz-Carpena, 2013) involve using the metric in conjunction with others (especially bias) to provide a more complete description of error and evaluating the robustness of model fit by constructing probability distributions of NSE through repeated fits to data subsets, as was done here.

Model Comparison
Overall performance in predicting the constructed instantaneous uptake rates (Figure 4), as measured by NSE, was somewhat better for the big root model (0.70) than either the parallel (0.67) or series (0.58) Ohm models. The parallel Ohm model performed nearly as well as the big root model in all layers except at the 50-cm depth, where it performed as poorly as the series model. All three models had trouble capturing the dynamics of the two deepest layers, where observations did not exhibit diurnal cycling of moisture loss.
The assumption of the calibration procedure that all water movement is through the vegetation was likely significantly off in these deeper layers, making it hard for all models to square physically consistent uptake in the overlying layers with these. The greatest differences among the models were apparent in the biases  Note. Distribution of fit statistics and prediction error in (n = 1,000) inversions to 72-hr subsets of data. NSE reflects the mean across all layers, bias is measured as the total absolute bias across all layers, and RMSE as the total across layers. Reported values of fit statistics are median and observed 95% confidence interval limits in parentheses. Computational cost shown for a single run of each model, showing cpu time in seconds of a single matrix inversion (flow calculation), the number of these evaluations in the entire run, and the resulting total time spent calculating flow by each model.
of these predictions, measured as the mean difference between prediction and data. On this metric, the big root model performed substantially better than the other two. Absolute bias summed across layers was 0.14 mm/day for the big root model, compared to 0.31 mm/day for the parallel and 0.55 mm/day for the series models, with the latter incurring an especially heavy penalty (0.20) in the deepest layer.
There were also significant differences between models in the robustness of these model fits (Table 3). The parallel model's fit to instantaneous uptake was particularly affected, as indicated by the lower end of the 95% confidence interval of its NSE values at 0.62. More than for the other two models, its performance was increasingly degraded for data sets using more clustered times, as these were increasingly less representative of the overall time series. This suggests a tendency of this model toward overfitting, which is a disadvantage if predictions are to extrapolate beyond the range of calibration data. 93.6% of all big root NSE values were higher than for the parallel model (99.6% than the series model). The parallel and big root models showed similar robustness in total absolute bias, with the latter again showing lower overall values. The series model showed slightly lower bias overall than with the larger data set, indicating very different fits with the limited data sets than the larger one.
These biases propagated to predicted soil volumetric moisture content ( Figure 5) and the resulting prediction errors ( Figure 6). The parallel model had total RMSE higher than the big root model in 96.6% of cases with a mean difference of 1.2% of volumetric soil moisture. The mean RMSE was separately lower for the big root model in each layer except at 50 cm. The series model not only had higher total RMSE than the big root in 99.4% of cases, with a mean difference of 2.7% of volumetric soil moisture, it also had an entirely different distribution of biases (and resulting errors) among layer depths, as compared to the baseline inversion. Its distributions of errors show a consistent positive bias in shallow and negative bias in deep uptake, resulting in a systematic dry bias in shallow layer soil moisture predictions and wet bias in deep layers. These biases are inherent to the model, which propagates negative water potentials from the root base linearly through shallow layers on to deeper ones. The absence of such biases in the fit to the longer, evenly distributed data set is notable, but clearly not robust, as the biases emerge in the prediction error distributions. A sensitivity analysis shows that these results do not change significantly in response to perturbations of the soil characteristic curve parameters (Text S1 and Figure S1).
The big root model and the series Ohm model assume the same RSA and differ only in its mathematical treatment-based either on an Ohm's law linearization of resistances or on seeking the mean root water potential while taking into account the nonlinearity of equation (1). Comparing the quality of resulting sets of predictions suggests that, for a given RSA, a porous pipe conceptualization is strictly better than Ohm's law, which leads to lower fits, greater biases and predictions that suffer from a rigid imposition of the linearized propagation of water potential gradients through the RSA.
The differences between the parallel and big root models are more complex. The parallel model makes a powerful simplification by assuming that water potential gradients propagate through woody roots with negligible hinderance, as this allows it to introduce the stem base boundary condition directly to all layers. The big root model replaces this parallel flow assumption by propagating the water potential gradient as it would in a single root, taking into account within-layer gradients and cross-layer effects. The present case study indicates that the latter is capable of somewhat better and more robust fits, consistently lower biases, and slightly but consistently better predictions, when both are fit by inversion. Taken together, these results suggest that its physical consistency makes the big root model less liable to overfitting and a better guide to predictions under nonstationary conditions than the parallel Ohm model.
Without fitting, as in LSM implementations, the parallel model has a known tendency to overpredict hydraulic redistribution, sometimes requiring drastic constraints to be placed on it (Kennedy et al., 2019). This follows from the ease with which water in this model simultaneously travels down the gradient from wet soil layers to the root base and from the base out to dry layers. It is a direct consequence of the assumed parallel resistance structure, which places the stem base boundary condition in direct contact with each soil layer. In a big root formulation, this issue is again strongly reduced: during uptake, water may be lost from the root to drier soil layers, but the root-soil gradient is always locally-not globally-in competition with the vertical root gradient, so that major root-soil losses can only occur once transpiration ceases.
The parallel model also produced substantial error in a comparison with analytical solutions for uptake on explicit 3-D RSA, when using actual per-layer resistance values (Bouda & Saiers, 2017). This suggests the model assumption of no effect of architecture on uptake was broken, at least for the root systems used in that exercise. The present big root approximation is the most constrained instance of the pentadiagonal model class, whose most general form led to substantially lower error in that case. Fitting unconstrained pentadiagonal stencils to the Wind River Crane 2010 drought data (Text S2 and Figures S2 and S3) resulted either in better fits to instantaneous uptake or better soil moisture predictions, but not both. Unconstrained, the model is thus more versatile and can represent a wider range of RSA, but is liable to overfitting and its predictions are not robust. The big root constraints on the RSA Stencil derived here sacrifice some of the flexibility of the general model class to achieve greater prediction robustness by enforcing physical consistency with the behavior of an actual root. Such robustness, based on accurate representations of mechanisms, is especially important in the context of future climate projections, where vegetation water uptake potentially takes place in no-analog scenarios.

Computational Cost
The computational cost of the prediction run with each of the models (Ohm's law series, parallel, and big root) is summarized in Table 3. With matrix inversion optimized out of the flow calculations, the cost per evaluation and the total cost of the big root model was within the range of values for the two Ohm's law models, showing broadly that the models involve comparable levels of computational cost. The computational cost in any LSM implementation is likely to be more sensitive to code optimization than model formulation. While there were 2,280 time steps in the simulation, the number of evaluations of each model are higher, as the nonlinears relationship was implicit in the equations solved, requiring an iterative solver. The number of evaluations of each model thus reflects the influence that the choice of root uptake model has on the convergence rate of the solution of the soil moisture characteristic curve. The greater number of evaluations for the parallel model reflects poorer convergence, which may translate into significant numerical errors when used with a direct solver that linearizes the soil hydraulics equations, as LSMs commonly do.

Limitations and Further Work
The main limitation of the big root model follows from the assumption that soil moisture is uniform in each layer. Lateral soil moisture heterogeneity is a key feature of unsaturated zone hydrology and neglecting it necessarily leads to model prediction error . This assumption is not novel in the big root model, however, but is instead a feature of all LSMs. In fact, it was assumed at the outset of model development to ensure compatibility of the result with relevant LSMs and maintain sufficiently low computational cost. The model might, however, be extended to account for lateral soil moisture heterogeneity, for example, by separately modeling a "dry" and "wet" soil column, with a two-root model simulating a split root situation. Microscale heterogeneity, between the rhizosphere and "bulk soil" can also be the key limiting factor to soil-root flows (Carminati et al., 2016). While it was not represented in this study, this resistance can be modeled separately, in series with the RSA model (Williams et al., 2001). Also, it is notable that not representing this mechanism in the big root model seemingly did little to degrade the quality of its predictions for the Wind River Crane data.
A limitation of the present implementation of all three models is that no roots are considered to occur below the unsaturated zone. This is especially problematic as groundwater likely influences vegetation water balance at sites with shallow water tables (Loheide et al., 2005;Nachabe et al., 2005) and seasonal droughts (Gou & Miller, 2014), as is the case here. Unfortunately, there were insufficient data to estimate the magnitude of uptake from the saturated zone, and the modeling exercise was thus constrained to the vadose zone alone, although that is also where errors due to the root model formulation are likely to be highest. Groundwater uptake can, nevertheless, be incorporated into vegetation water uptake models using either the Ohm's law (Gou & Miller, 2014) or the present formulation, and has been productively implemented in land surface models (Wang et al., 2018) and used to study plant-water feedbacks in a climatic context (Fan et al., 2017).
To allow for broader application of a big root model in LSMs, its parameter values need to be established for a number of sites with distinct plant functional types. The calibration procedure is data intensive and subject to some restrictive assumptions. It is particularly important to arrive at an adequate partitioning of root water uptake from soil water movement prior to calibration. Model calibration is also sensitive to the soil moisture characteristic (s ) curve, as the model is linear in soil water potential, which is itself very nonlinear in soil moisture. Ideally, the calibration data set would contain both sets of values, directly observed. It should be noted that calibrating Ohm's law models requires the same set of inputs as for the big root model. The latter requires the estimation of at most one more parameter per layer, as it decomposes root conductances into conductivities and length.
Given its physical basis, it is possible to use the model to find meaningful hydraulic parameters of the "effective root" (i.e., total within-layer soil-to-root conductance and the cross-layer root conductance) rather than just constrained RSA Stencil coefficients that predict water uptake, as was done for the test case. The necessary condition is to have data on xylem water potential at the base of the root system for each time step, which can be obtained through stem psychrometry (McBurney & Costigan, 1982). Obtaining time series of how such effective big root parameters change at different time scales in real vegetation may lead to insights concerning the importance of underlying biological mechanisms (e.g., root growth and senescence, aquaporin activity) at larger scales. As it is based on an improved conceptualization of root water uptake, the big root model can be expected to yield more accurate estimates of per-layer root conductances than a comparable Ohm's law model, using the same calibration data.
Finally, following from the observation that pentadiagonal models can be fit to the uptake of complex RSA and the conjecture that the big root model may be one of a class of such models related analytically to specific assumed RSA configurations, further work in this direction is warranted. The aim would be loosening the big root constraints by establishing analytical expressions for coefficients of pentadiagonal models for more complex RSA. This should yield greater model versatility without sacrificing physical consistency.

Conclusions
A big root model for site-scale vegetation water uptake can be derived from first principles describing water movement through plant roots. The novel formalism provides a simple and mechanistically based framework to describe a key process limiting land surface exchanges of matter and energy. Considering the nonlinearity of potential gradients on this effective root provides more robust predictions of water uptake than commonly used Ohm's law analog models. Robust, physically based models are especially valuable in forecasting applications with nonstationary climate, where empirical or heuristic formulations are likely to fail. This makes the present conceptual advance a promising line of inquiry to address the persistent knowledge gap concerning root water uptake and its limiting influence on transpiration. As it is derived under the very restrictive assumption of a single, big root, one way to further improve the model may be to generalize this approach to other instances of root system architecture. Potential applications of the theory beyond terrestrial modeling include improving drought stress formulations in existing crop models, which often rely on empirical or heuristic formalisms to predict crop yields, and thus food production, in future climate scenarios. I am deeply grateful to Sonia Wharton and the team at Wind River Crane who collected and kindly provided the data underlying the test case. I gratefully acknowledge my colleagues in the Herben and Brodersen research groups, as well as two anonymous reviewers, all of whose comments improved this manuscript. This work was supported by the Ministry of Education, Youth, and Sports of the Czech Republic through the institutional project MŠMT CZ.02.1.01/0.0/0.0/16 019/0000803 and the Large Infrastructures for Research, Experimental Development and Innovations project "IT4Innovations National Supercomputing Center LM2015070". All original results reported here and original code implementing the models and their inversions are available at the public repository (https://github.com/ mbouda/bigRootManuscriptCode).

Raw observations for the Wind River
Creek site can be obtained from the present PI of the project.