Enhancing the Noah‐MP Ecosystem Response to Droughts With an Explicit Representation of Plant Water Storage Supplied by Dynamic Root Water Uptake

Plants are able to adapt to changing environments and thus survive droughts. However, most land surface models produce unrealistically low ecosystem resiliency to droughts, degrading the credibility of the model‐predicted ecohydrological responses to climate change. We aim to enhance the Noah‐MP modeled ecosystem resilience to droughts with an explicit representation of plant water storage supplied by dynamic root water uptake through hydrotropic root growth to meet the transpiration demand. The new model represents plant stomatal water stress factor as a function of the plant water storage and relates the rate of root water uptake to the profile of model‐predicted root surface area. Through optimization of major leaf, root, and soil parameters, the new model improves the prediction of leaf area index, ecosystem productivity, evapotranspiration, and terrestrial water storage variations over most basins in the contiguous United States. Sensitivity experiments suggest that both dynamic root water uptake and groundwater capillary rise extend the plants' “memory” of antecedent rainfall. The modeled plants enhance their efficiency to use antecedent rain water stored in shallow soils mainly through more efficient root water uptake over the U.S. Southwest drylands while use that stored in deep soils and aquifers with the aid of groundwater capillary rise in the Central United States. Future plant hydraulic models should not ignore soil water retention model uncertainties and the soil macropore effects on soil water potential and infiltration.


Introduction
Arid and semiarid (or dryland) ecosystems cover~40% of the global terrestrial surface (Reynolds et al., 2007) and are expected to expand under projected warming climates (Huang et al., 2015). Also, droughts may become more frequent, intense, and last longer (Trenberth et al., 2014). Due to increases in the atmospheric demand in terms of potential evapotranspiration (PET), the global aridity index (defined as the ratio of PET to precipitation) will increase, and the global drylands are expected to expand (Dai, 2012;Feng & Fu, 2013;Huang et al., 2015) despite uncertainties caused by the way how PET is computed  or whether plants' stomatal closure as a result of CO 2 fertilization is considered (Milly & Dunne, 2016). Therefore, terrestrial ecosystems over drylands under persistent water stress and those in dry-to-wet transition regions under episodic water stress will suffer from more frequent and intense droughts.
In dryland ecosystems or ecosystems during droughts, water becomes a dominant resource limiting biomass productivity. In general, ecosystem gross primary productivity (GPP) and net primary productivity (NPP) increase as the mean annual precipitation (P) increases. For mesic ecosystems, the rain use efficiency (RUE), that is, biomass productivity per unit rainfall (GPP/P or NPP/P; g C/kg H 2 O), decreases with increasing P, as other resources such as nutrient and light become limiting factors (with more water loss through runoff and evaporation). However, the ecosystems' RUE tends to increase and converge across biomes to a common maximum value during the driest conditions (Huxman et al., 2004;Ponce-Campos et al., 2013). Dryland ecosystems show much stronger resiliency or higher capacity to tolerate water stress as indicated by the widespread woody encroachment (Andela et al., 2013), dryland greening (Fensholt et al., 2012), and increasing trends in the net carbon sink in semiarid regions (Ahlstrom et al., 2015) despite warming associated droughts.
Most state-of-the-art land surface models (LSMs), however, produced unrealistically low ecosystem resiliency in terms of RUE under water stress. For instance, the DOE ELMv1 produced much lower RUE under water stress, and so do versions of the NCAR Community Land Model (CLM) (Zhu et al., 2019). An earlier version of CLM (CLM4.5) produced a decreasing trend in leaf area index (LAI) larger than remote sensing observations in the drying drylands (Mao et al., 2013). Also, Noah-MP  produced lower-than-observed GPP, LAI, and evapotranspiration (ET) during droughts in the Central United States (N. . One major limitation to most LSMs is the use of prescribed, static root profiles as a uniform, exponential, or asymptotic profile (Hao et al., 2005;Smithwick et al., 2014) through type-dependent parameters controlling maximum rooting depth and vertical root distribution (Jackson et al., 2000;Zeng, 2001). The static root profile disconnects the interactions between changes in belowground water and nutrient resources and aboveground plant carbon assimilation.
Plant roots are known to grow unevenly in both lateral and vertical directions toward subsurface water through hydrotropism (Dietrich et al., 2017) or hydropatterning, enabling roots to perceive microscale water potential gradient (Bao et al., 2014). Soil water profile becomes a dominant factor regulating plants' rooting depth over many other factors (Fan et al., 2017). While many roots penetrate soil concretions and bedrock fractures and terminate in the capillary fringe to avoid oxygen stress (Fan et al., 2017), phreatophytes, widely distributed in global drylands, extend their aerenchymatous roots deep into the saturated zone (or phreatic zone) (Gary, 1963;Naumburg et al., 2005;Scott et al., 2014;Stromberg, 2013), extracting water directly from groundwater and capillary fringe (Ehleringer & Dawson, 1992;Le Maitre et al., 1999;Scott et al., 2014). In response to declines in the groundwater level, rapid root growth (up to 15 mm/day) toward the water table has been observed for desert phreatophytes (Naumburg et al., 2005;Orellana et al., 2012;Vonlanthen et al., 2010). A recent modeling study of phreatophytes at a riparian site under hyper-arid climates (with P of only~35 mm/year) indicates that direct root water uptake from the capillary fringe and groundwater (~6 m deep), which is supplied by lateral saturated flow, accounts for~84% of the total transpiration (~472 mm/year), about 10 times P (Wang et al., 2018). In addition, dryland ecosystems have evolved thinner roots to efficiently enhance root length and surface area, improving the efficiency of carbon invested in roots (Z. . There are two general approaches to modeling root water uptake at microscopic and macroscopic scales (Feddes & Raats, 2004;Šimůnek & Hopmans, 2009;Warren et al., 2015). Microscopic models describe water extraction by individual roots in cylindrical (or radial) coordinates with three-dimensional (3-D) root architecture models that predict evolution of the 3-D branching architecture, for example, ROOTMAP (Diggle, 1988), SimRoot (Lynch et al., 1997), and SPACSYS (Wu et al., 2007) and others as reviewed by Dunbabin et al. (2013). The 3-D coupled root-soil water model, R-SWMS (Draye et al., 2010;Javaux et al., 2008), which couples a 3-D flow and transport model (Somma et al., 1998) and a 3-D model of root architecture dynamics and root hydraulics (Doussan et al., 1998), represents one of the most complex microscopic models. Macroscopic models describe the lumped effects of vertical distribution of root density in length (Hartmann et al., 2018) or surface area (Schymanski et al., 2008) on root water uptake from different soil layers. To save computational costs, most recent developments of dynamic root water uptake models in LSMs (e.g., Drewniak, 2019;Gayler et al., 2014;Wang et al., 2018;Zhu et al., 2017) followed the macroscopic approach for use at landscape, regional, and global scales (Feddes et al., 2001).
Like many other LSMs, Noah-MP predicts bulk root biomass as the residual of photosynthate after allocation to leaves and wood for balancing the carbon budgets but does not link the root biomass, structure, and function to soil water distributions. Following the macroscopic approach, Wang et al. (2018) implemented the dynamic root water uptake model of the vegetation optimality model (VOM) developed by Schymanski et al. (2008Schymanski et al. ( , 2009 into Noah-MP  and applied the coupled model at a field site. In this study, we develop a layered root biomass model that it allows more carbon translocation to roots when the plants are under water stress and represents root hydrotropism through more carbon translocation to roots in wetter soil layers but with less turnover following Parton et al. (1978) (Figure 1). Also, we modified the carbon allocation scheme to allow roots to first receive photosynthate followed by leaves, stems, and wood (Wu et al., 2007). We then used the vegetation-type-dependent specific root area data of Jackson et al. (1997) to convert the layered root biomass to layered root surface area for use in the root water uptake model. Following VOM, the new model explicitly predicts plant water storage, which is demonstrated to buffer water stress during daytime through nocturnal water uptake by roots (Huang et al., 2017) and possibly at seasonal scales in some tropical woodlands (Tian et al., 2018). We then used the plant water storage to parameterize the water stress factor for plant photosynthesis and stomatal conductance. We tested the model performance over the contiguous United States (CONUS) river basins and investigated, through sensitivity experiments with the new model, the role of root water uptake in buffering drought stress with a focus on the U.S. Southwest drylands under persistent droughts and the Central United States that experiences more frequent and intense episodic droughts at river basin scales. We use monthly RUE and its response to aridity index as a critical measure of the modeled ecosystem resiliency to droughts, because monthly RUE measures how efficiently plants use rain water stored in the subsurface during previous months. Throughout the paper, we use soil matrix water potential, soil matrix suction head, and soil water pressure interchangeably.

Model Description 2.1.1. Plant Water Storage and Root Water Uptake
The dynamic root model explicitly predicts plant water storage and root water uptake, which is related to the root surface area (Figure 1). The governing equation for the total amount of water stored in the living plant tissues, M q (kg m −2 or mm), is where Δt is the time step size (= 3,600 s in this study), Q R is the rate of water uptake by roots (kg m −2 s −1 or mm s −1 ), and Q T is transpiration rate (mm s −1 ). Q R is the sum of root water uptake of each root layer, Q R,i , which is related to the root surface area in the layer, A R,i , and the hydraulic gradient between the matric suction head, h s,i (mm), and the root suction head, h R (mm), following Schymanski et al. (2008) (Figure 1): Figure 1. Diagram of the plant water storage (M q ) and root water uptake (Q R ) model: (1) The model explicitly predicts changes in M q as a residual of Q R and transpiration (Q T ); (2) Q R , controlled by root surface area (A R,i ), is converted from live biomass of fine roots in a soil layer (C R,i ); (3) the stomatal water stress factor (β) is parameterized as the relative excess of M q over the plant wilting water storage; (4) it allows more carbon translocation to roots (F R ) when the plants are under water stress [F R = 0.3(1 − β)]; and (5) it represents root hydrotropism through more carbon translocation to roots weighted by a wetness factor (w i ) but less turnover (D w,i ) in wetter soil layers.
where N R is the total number of soil layers that contain roots, A R,i is the root surface area per unit ground area (m 2 m −2 ) or root area index (RAI), which is converted from the model-predicted root carbon mass of the ith layer, and C R,i (g m −2 ): A R,i = A R,s C R,i /1,000, where A R,s is specific root surface area per unit carbon mass [m 2 /(kg C)] for different vegetation types (Jackson et al., 1997; see Table 1). Ω R (s) is the root radial resistance to water flow through the roots (Table 1), and Ω s,i (s) is the resistance to water flow through soil matrix from ambient soils toward the root surface in the ith layer (Appendix A). h R is the root suction head, and h s,i = h sat (θ i /θ sat ) −b (m) is the soil matrix suction head, where θ i and θ sat are liquid soil moisture (m 3 m −3 ) in the ith layer and that at saturation, respectively (Clapp & Hornberger, 1978). Note that, in Equation 2, the elevation head is neglected for both h s,i and h R .
We assume that all parts of the plants are in an equilibrium hydraulic state, neglecting the elevation head gradient between leaves, stems, xylem, and layers of roots. Therefore, the root suction head h R approximates the leaf water head, and h R = −c bm P b , where c bm = 10.2 m bar −1 , and P b (bar) is the leaf balance pressure, that is, the pressure needed to squeeze the water out of leaves (Roderick & Canny, 2005): where M dry is dry biomass (kg m −2 ) and M q,max is the maximum water storage in the plant living tissues when the plants are at full hydration. Supported by measurements of 156 leaf samples of 10 plant species with the maximum leaf water fraction [M q,max /(M q,max + M dry )] ranging from 52% to 86%, Roderick and Canny (2005) suggest that P b shows a linear relationship with the relative leaf water loss, (M q,max − M q )/ (M q,max + M dry ), at a slope of κ [bar/(kg m −2 )], which is generally decreasing with M q,max . We take the form of Schymanski et al. (2008) for the slope, κ, which was obtained by fitting the data of Roderick and Canny (2005): where r w = M q,max /M dry , a vegetation-type-dependent parameter (Table 1), and the two constants: c 1 = 750.0 bar, and c 2 = 1.0 bar. Neglecting the smaller second term, Equation 4 suggests that dryland species with a smaller r w would produce a greater root suction head (h R ) when losing the same percentage of water.  ); δ T = the decay rate of D T,max (-); L TO = leaf turnover rate (10 −6 s −1 ); A leaf,s = specific leaf area A leaf,s (m 2 kg −1 ). a A R,s values are adopted from Jackson et al. (1997) and slightly revised, and the rest are calibrated against MODIS leaf area index (LAI) and FLUXNET MTE GPP and ET data sets (see section 2.3).
The dry biomass, M dry = 2.5 × (C leaf + C stem + C root + 0.02 × C wood ), where C leaf , C stem , C root , and C wood are the carbon mass stored in leaf, stem, root, and wood (assuming 2% sapwood), which are predicted by the dynamic vegetation module of Noah-MP. The factor 2.5 is the ratio of dry biomass to carbon mass (1 kg C = 2.5 kg biomass; through the formulae of sugar, C 12 H 22 O 11 , while considering other heavier tracer elements; Parton et al., 1978).

Root Carbon Mass
The governing equation for the carbon mass stored in the ith root layer, C R,i (g m −2 ), is where R GPP,i is photosynthate (or GPP) allocated to the ith soil layer roots (g m −2 s −1 ), R EX,i root exudates (g m −2 s −1 ), which is simply assumed 30% of R GPP,i , R MR,i root maintenance respiration (g m −2 s −1 ) based on a Q 10 function of soil temperature (Q 10 = 2.0), R GR,i root growth respiration (g m −2 s −1 ), and R TO,i root turnover due to temperature and water stresses (g m −2 s −1 ). R GPP,i is proportional to the layer thickness, Δz i (m), and a weighting factor, w i : where A is the total carbon assimilation rate of the sunlit and shaded leaves (g m −2 s −1 ) through photosynthesis . F R is the fraction of A partitioned to roots, where β is the plant stomatal water stress factor (section 2.1.3), representing that plants tend to partition more carbon to roots under water stress with a maximum of 30%. The weighting factor, w i , is parameterized as a function of layer depth and θ i (m 3 m −3 ): where z i is the layer node depth, θ wilt the wilting soil moisture (m 3 m −3 ), and θ ref a reference soil moisture close to field capacity (m 3 m −3 ). Equation 7 represents two principles of root carbon allocation: (1) more carbon to shallower layers and (2) more to wetter layers (Parton et al., 1978).
The root turnover rate R TO,i (g m −2 s −1 ) of soil layer i is computed as a function of root death factors due to soil temperature stress, D T,i , and moisture stress, D W,i , in layer i: where C R,min = 0.01 g m −2 , and R TO,max is the maximum root death rate in extremely dry and unfavorable temperature conditions (Table 1): We use a constant temperature factor, D T,i = 0.005, not following the functional form of Parton et al. (1978), which may produce too little root biomass in too cold or too warm conditions.

Plant Stomatal Water Stress Factor
We further develop a new scheme of plant water stress factor for plant carboxylation and stomatal conductance as a function of plant water storage (or plant water potential): where M q,wilt is the plant wilting point in water depth (mm or kg m −2 ). Equation 10 means that β is actually the plant water available for transpiration (M q − M q,wilt ) relative to its maximum (or the maximum water that a plant can lose through transpiration until its wilting point; M q,max − M q,wilt ). M q,wilt can be then converted from the plant wilting point in pressure through Equations 3 and 4: For a plant with r w = 9.0, M q,wilt = 0.562 M q,max , meaning that a plant at full hydration can lose 43.8% of its water until the wilting point. Also, Equation 11 indicates that a plant with a greater r w can lose a greater percentage of its water until the wilting point. In this study we simply take P b,wilt = 30 bar or 3.0 MPa, representing dryland phreatophytes (Wang et al., 2018), whereas it may vary widely with plant types (Bonan et al., 2014).
Because of the linear relationship between leaf water storage and leaf water potential (Equation 3), Equation 10 can be easily transformed into a function of leaf water potential, P b , We assume Q R,i = 0, when h s,i at a soil layer reaches the wilting point, reducing M q to M q,wilt and thus resulting in a complete stomatal closure.

Data 2.2.1. Climatic Forcing, Vegetation, and Soil Parameters
We used the hourly, 0.125°North American Data Assimilation Phase 2 (NLDAS-2) near-surface meteorological forcing data of air temperature, specific humidity, wind speed, surface pressure, downward shortwave radiation, downward longwave radiation, and precipitation (Xia et al., 2012) to drive the model. We used the global 1-km hybrid State Soil Geographic Database and the United States Geological Survey (USGS) 24-category vegetation data to determine the soil and vegetation parameters, respectively. Both the 1-km data sets are aggregated to 0.125°with the dominant soil and vegetation types to match the spatial resolution of the NLDAS-2 forcing data ( Figure 2). The soil and vegetation parameters are then determined for each soil type and vegetation type through the look-up tables of Noah-MP.

FLUXNET MTE GPP and LE
To calibrate and validate the modeled ET and GPP, we used the FLUXNET model tree ensembles (MTEs) GPP and latent heat (LE) from the Max Planck Institute for Biogeochemistry. The MTE was first trained with measured flux data of 198 FLUXNET towers across a wide range of biomes worldwide. With inputs of fraction of absorbed photosynthetic active radiation derived from remote sensing, climate, and land cover data, this approach generates monthly, 0.5°gridded GPP and LE data during 1982-2011 over global continents (Jung et al., 2010(Jung et al., , 2011. While the uncertainty in the MTE products is not negligible due to the density of flux towers selected for training the model tree, we have the greatest confidence over the "data-rich" CONUS, where most AmeriFlux sites were incorporated (Jung et al., 2010).

MODIS LAI
We calibrated and validated the modeled LAI against a global MODIS LAI product that were reprocessed considering the presence of cloud and seasonal snow cover, instrument problems, and uncertainties of retrieval algorithm in the original product (Yuan et al., 2011). We have collected 15-year (2002-2016) global 0.05°reprocessed MODIS LAI product. These improved MODIS LAI data are closer to LAI reference maps in magnitude and also more continuous and consistent in both time series and spatial domains.

GRACE Terrestrial Water Storage Anomaly
Changes in terrestrial water storage (TWS) of soil water, groundwater, and surface reservoirs are greatly affected by terrestrial surface ET. As an indirect validation of the modeled ET, we compared the modeled TWS anomaly (TWSA) to the 1°, monthly Gravity Recovery and Climate Experiment (GRACE) TWSA products released by three different processing centers. We used the mean of the three products, because the noise of different solutions can be effectively reduced by arithmetic average of the three products (Sakumura et al., 2014). Also, we used the gain factor to compensate for the signal loss during postprocessing of the original surface mass variations (Landerer & Swenson, 2012). According to Kumar et al. (2016), the total measurement error of the GRACE TWSA ranges from 0 to 40 mm over most parts of CONUS, but 10.1029/2020MS002062

Journal of Advances in Modeling Earth Systems
higher errors may occur on the West Coast, the lower Mississippi River basin, and Florida due to the spatial leakage error (Seo et al., 2006).

USGS Water Watch Hydrological Unit Runoff
We used the monthly USGS Water Watch hydrological unit runoff data. USGS divides CONUS into successively smaller hydrological units with a unique "hydrological unit code" (HUC; Seaber et al., 1987). The USGS runoff of each hydrological unit was generated by combining the historical flow data collected at stream gauges, drainage basins of the stream gauges, and boundaries of the hydrological units (Brakebill et al., 2011). This data set contains monthly runoff from 1901 to 2016 at various basin scales from 18 first-level two-digit basins (HUC2; Figure 2) to eight-digit basins (HUC8) and has been regarded as a close surrogate of natural runoff for hydroclimate studies (Ashfaq et al., 2013;Ma & Szilagyi, 2019;Oubeidillah et al., 2014;Velpuri et al., 2013). The HUC2 basins are further divided into 334 third-level six-digit basins (HUC6) (Seaber et al., 1987). In this study, we used the monthly HUC6 runoff from 1982 to 2008 to derive an ET climatology for use to evaluate the modeled and FLUXNET ET climatology.

Model Experiments
We first conducted two model experiments: one with a previous version of Noah-MP (hereafter OLD) and one with the new developments described in section 2.1 (NEW) to test the models' performance. To discern the exact mechanisms leading to the differences between the two models, we then conducted sensitivity experiments using the new model. All the model experiments were conducted at hourly time step from 1 January 1980 to 31 December 2015 driven by the NLDAS-2 forcing under the CO 2 level of 360 ppm and spun up for two loops (totally 72 years) to reduce the effects of unknown initial conditions, and the modeling results from the second loop were analyzed.
Experiment OLD, which includes all the augmentations in Niu et al. (2011), uses the same model physical options as the EXP6 in Yang et al. (2011) and later evaluated by N.  (see options listed in Table 2). This version together with the manually optimized parameters against various data sets of runoff, TWSA, LAI, and greenness fraction  has been released by NCAR and widely used in many other applications, for example,   . In this study (also in EXP6 of Yang et al., 2011), we choose the Noah water stress factor as a function of soil moisture: where z R is the rooting depth. The plant transpiration controlled by β is evenly partitioned into different root layers and thus proportional to layer thickness, being treated as root water uptake. For lack of subsurface data, some key hydrological parameters are assumed spatially constant over the domain, for instance, the maximum groundwater discharge rate, the saturated hydraulic conductivity decay factor, and the micropore volume fraction (f mic ) (see more details in Niu et al., 2011). Noah-MP adopted the simple bucket-type groundwater model of Niu et al. (2007) to represent recharge into the "bucket" in wet periods and capillary rise of groundwater from the "bucket" during dry periods. It also introduced f mic to reduce the capillary rise of groundwater to account for the presence of soil macropores and thus helped improve the modeled soil moisture variability in the State of Illinois . f mic ranges from 0.0 to 1.0, and f mic = 0.0 represents a free drainage lower boundary condition at the 2 m deep bottom of the model soil column, while f mic = 1.0 represents a full effect of capillary rise of groundwater. In general, a larger f mic produces a wetter soil with smaller soil moisture variability. In Experiment OLD, f mic = 0.2, a constant over the whole CONUS domain.
On top of OLD, Experiment NEW includes all the new developments described in section 2.1 (also the Appendixes A and B) and thus differs only in the "dynamic vegetation" and the "water stress factor for stomatal conductance" as listed in Table 2. We adjusted the spatially constant f mic from 0.2 in OLD to 0.6 in NEW over the whole domain. We calibrated the vegetation parameters against the data sets (section 2.2) averaged over each type and assessed the modeling results using the same data sets but averaged over each HUC2 and HUC6 basins. Through manual optimization (incomplete search of optimum parameter values), we first calibrated the vegetation parameters to minimize the model errors in the modeled monthly GPP, LE, and LAI, which are averaged over the same vegetation types in CONUS (Supporting Information Figures S1-S3) as measured by the Nash-Sutcliffe efficiency (NSE). We also calibrated root turnover parameter (R TO,max in Table 1) against total root biomass data of Jackson et al. (1997) for major vegetation types. The optimized model performs fairly well for most vegetation types except "irrigated crop" and "shrubland." We then calibrated the soil parameters of only "bedrock" for the whole soil profile, which are adjusted closer to sand parameters (to hold more water for plant use) against the monthly GPP, LE, and LAI data sets (section 2.2). The "bedrock" exists mainly over the U.S. Southwest, for example, the Rio Grande and the Lower Colorado, where there are some shrubs and grasses (Figure 2). The calibration is challenging because a same vegetation type (e.g., shrubland, grass, and evergreen needle) covers a large area across various climate zones. For instance, the evergreen needle covers a large area of both the Southeast and Northwest (Figure 2), which experience totally different seasonality in surface air temperature. For this reason, we developed a temperature-sensitive function for leaf turnover (D T,leaf in Equation B3), which is critical for improving the simulation of LAI seasonality in different climate zones. The resulting optimum vegetation parameters for major vegetation types are listed in Table 1.
Due to the many differences in process representations and associated parameters between OLD and NEW, it becomes difficult to interpret the differences between NEW and OLD. Therefore, we conducted sensitivity experiments using the new model by changing a single process representation to mainly isolate the combined effects of two major "pumping" mechanisms for plants to survive droughts: plant root water uptake and groundwater capillary rise. The first experiment (STATIC) reverts the explicit representation of plant water storage and dynamic root water uptake in NEW back to the Noah water stress factor for stomatal conductance as a function of soil moisture (Equation 13) and a static, evenly distributed vertical root profile used in OLD, and f mic remains the same value as in NEW (f mic = 0.6). The second experiment (FD) changes f mic = 0.0 to represent free drainage or zero groundwater capillary rise while retains all other schemes as in NEW. In all the experiments, N R = 4 (rooting depth = 2.0 m), while we will 10.1029/2020MS002062

Journal of Advances in Modeling Earth Systems
discuss the role of rooting depth in section 4 as inspired by an additional sensitivity experiment with N R = 3 (rooting depth = 1.0 m).

Modeled LAI, GPP, and ET
We evaluated the model performance over CONUS at a grid scale in terms of absolute and relative biases in the climatological mean (Figures 3 and 4) and over the HUC6 basins with NSE of the monthly outputs ( Figure 5). The same as in N. , the modeled ET biases by OLD show a similar pattern with those of GPP, that is, more ET and GPP in the Eastern United States and less ET and GPP over the Southeast coastal regions and the California (Figures 3e and 3h). OLD overestimates GPP by 40% and ET by 22% over the whole CONUS domain, while with a prescribed monthly climatology of LAI, the relative bias in ET drops to 4%, indicating that the dynamic vegetation module in Noah-MP is problematic (N. . OLD underestimates LAI in the Eastern United States and the coastal West and slightly overestimates LAI in the Central United States, a pattern quite different from that of GPP. This suggests that LAI and GPP are not necessarily linked in the model due to different major controls on leaf dynamics and photosynthesis. In the Northeast United States, OLD overestimates GPP due mainly to overestimation of maximum rate of carboxylation at 25°C (V c,max in Table 1), while it underestimates LAI due mainly to overestimation of the leaf turnover rate.
NEW largely reduces the model biases in LAI, GPP, and ET by OLD (Figures 3c, 3f, and 3i). However, NEW still produces positive biases in GPP and ET over the High Plains and the Lower Mississippi basins. Also, NEW produces negative biases in GPP, ET, and LAI over the northern energy-limited basins, for example, the Great Lakes, Rockies, and the coastal Northwest as well as southern water-limited regions, for example, the southern tip of the Texas-Gulf basin and California.
To more precisely evaluate the modeled ET climatology (30-year mean from 1982 to 2011), we then compared the modeled and FLUXNET ET climatology to an estimate derived from the water balance, ET wb = P − R, where P and R are NLDAS precipitation averaged over the HUC6 basins and USGS HUC6 runoff, respectively. This estimate, ET wb , assumes that (1) the water storage change in a basin is negligible compared to cumulative fluxes over the 30-year long-term period and (2) the inter-basin water transfers (e.g., through irrigation channels or lateral groundwater flow) are negligible. Due to the higher accuracy in measuring streamflow (~5%) than measuring ET with eddy-covariance sensors (~5-20% at landscape scale), ET wb may represent the "best" estimate of ET climatology and was used to evaluate other estimates including the FLUXNET (Ma & Szilagyi, 2019). FLUXNET underestimates ET over most of the basins in the southern United States, more apparently over basins in the High Plains and the southwest United States (the Lower Colorado and Rio Grande basins) by up to 20% (Figure 4d). It overestimates ET in the Great Lakes and southern tip of Texas by 5-10% and more significantly over the Pacific Northwest and southern California by more than 20% (Figure 4d). Both the absolute and relative ET biases produced by NEW appear smaller than those of the FLUXNET in the southern United States (Figures 4b and 4e). NEW ET is better than the FLUXNET ET by about 10-20% over the water-limited basins of the High Plains and Southwest (Figure 4f). However, it is apparent that NEW underestimates ET and overestimates runoff ( Figure S7) over basins along the coastal East, Great Lakes, Rockies, and California. Except California, where ET is underestimated due most likely to neglect of irrigation by the model, the basins along the coastal East, Great Lakes, and Rockies are relatively wetter, and the soil types are mostly sandy soils (including sand, loamy sand, and sandy loam; see Figure 2b), which may facilitate drainage of water, resulting in more runoff ( Figure S7).
We also evaluated the model temporal variability in LAI, GPP, and ET using NSE over the 334 HUC6 river basins ( Figure 5). OLD predicts ET (Figure 5g) better than GPP ( Figure 5a) and LAI (Figure 5d) with NSE > 0.4 over most of the HUC6 basins. But GPP and LAI modeled by OLD are unacceptable with NSE < 0 over most of southern CONUS basins including water-limited basins (e.g., Arizona) or basins experiencing episodic water stress (e.g., Texas). NEW largely improves the simulation of GPP and LAI over these basins (Figures 5c and 5f) and slightly improves the simulation of ET over most of the CONUS basins ( Figure 5i). However, NEW still produces very low NSE values of LAI, GPP, and ET (NSE < 0.0) at HUC6 basin scale over the U.S. Southwest (Figures 5b, 5e, and 5h), because it is more challenging for a model to capture the large interannual variability in GPP, LAI, and ET, in response to the large temporal variability in the climatic forcing (e.g., the monsoon rainfall). However, NEW performs fairly well at HUC2 basin scale over the Lower Colorado and Rio Grande basins in the U.S. Southwest ( Figure S8).

Modeled TWS in the Central United States
Seasonal variations in TWS are greatly controlled by seasonal variations in ET; a low bias in the modeled ET due to low transpiration would result in a shallower trough in the modeled TWS. Therefore, we evaluated the modeled TWSA against the GRACE product as an indirect evaluation of the modeled ET. The modeled TWS is the sum of the modeled groundwater storage, snow water equivalent, soil moisture, and canopy-intercepted water. NEW also includes the plant water storage (M q in Equation 1). Consistent with postprocessing of the GRACE TWSA, we also removed the static field averaged over all months from 2004 to 2009 from the TWS time series. OLD produces shallower troughs during droughts, more apparently during the widespread 2012 Central U.S. drought and the 2011 Texas drought, and higher ridges in relatively wetter periods (left panels of Figure 6). NEW largely improves the simulation of TWSA during the droughts, resulting in greater NSE values except the Arkansas-White-Red basin, where the NSE value drops to 0.45 from 0.75 (see section 3.3 for the causes). The difference in the modeled TWSA between OLD and NEW (shaded areas in Figure 6) can be mainly attributed to the cumulative effect of the modeled ET; NEW produces more ET during drier seasons (blue shaded area in the right panels of Figure 6) but less during wetter seasons (green shaded area). Also consistent with the modeled TWSA, NEW improves the simulation of ET with greater NSE values except the Arkansas-White-Red basin, where the ET NSE value drops to 0.81 from 0.85. This suggests that NEW produces more persistent ET during droughts than OLD due presumably to a longer "memory" of antecedent rain.
The mechanisms resulting in the more persistent ET by NEW can be further attributed to enhanced transpiration (left panels of Figure 7) and GPP (right panels of Figure 7) during dry seasons. NEW produces more transpiration, and GPP does OLD during dry seasons, more pronounced in the Lower Mississippi, Arkansas-White-Red, and Texas-Gulf basins (blue shaded areas in Figure 7) and over all the basins during the 2012 summer drought in the Central United States. NEW largely improves the simulation of GPP as indicated by its NSE values compared to those of OLD. However, over the Texas-Gulf basin, GPP produced by NEW is still lower than the FLUXNET GPP over extremely dry summers of 2006, 2009, and 2011 (see section 4 for the causes). All the GPP, transpiration, and ET fluxes show the same pattern of differences between

Journal of Advances in Modeling Earth Systems
NEW and OLD (Figures 6 and 7). This suggests that, compared to OLD, NEW provides a "memory" mechanism for the plants to sustain their productivity (GPP) with a greater transpiration (Figure 7), resulting in deeper troughs in TWS during droughts ( Figure 6).

Mechanisms of the Modeled Ecosystem Responses to Droughts
NEW largely improves the modeling of GPP and LAI over almost the whole domain ( Figure 5) and ecosystem responses to the varying water stress in the Central United States (Figures 6 and 7). Through analyses of the sensitivity experiments, we learned that these improvements are through two different "pumping" mechanisms that sustain the plant water storage during droughts: root water uptake and groundwater capillary rise in different regions.
We present the results from the sensitivity experiments with a focus on the dryland ecosystems over the U.S. Southwest HUC2 basins including the Rio Grande and Lower Colorado under persistent water stress ( Figure S9a) and ecosystems in the Central U.S. HUC2 basins including the Arkansas-White-Red and Texas-Gulf under episodic water stress ( Figure S9d). We use the response of monthly RUE (= monthly GPP/monthly P) to monthly aridity index (AI = monthly PET/monthly P), of which PET is computed with the Penman-Monteith equation assuming that surface resistance is zero, as a critical measure of ecosystem resiliency. We set up a minimum value of 0.01 mm/month for a month when the monthly P is zero. Both the modeled and the FLUXNET monthly RUE closely follow AI, Figure 4. ET biases relative to ET wb of (a) FLUXNET ET (mm/day) and (b) NEW ET (mm/day) and (c) bias reduction by NEW (i.e., the difference between the absolute value of FLUXNET ET bias and that of NEW ET bias; NEW is better than FLUXNET over "blue" basins while worse over "red" basins) as well as their corresponding relative biases (d and e; %) and relative bias reduction (f; %) over the HUC6 basins.
increasing almost linearly with AI (Figures 8 and 9), suggesting that, during a drier month with little rain (contributing to a higher AI), plants tend to use antecedent rain water stored in the soil under favorable light and temperature conditions (also contributing to a higher AI) for photosynthesis. Therefore, the slope of RUE against AI (Figures 8c, 8d, 9c, and 9d) reflects the strength of a plant's "memory" of antecedent rain water stored in the soil or the strength of ecosystem resiliency to the present water stress.
Over the two representative semiarid U.S. Southwest HUC2 basins (Figure 8), the RUE difference between NEW and STATIC is much larger than that between NEW and FD ( Figure S10), suggesting that the root water uptake plays a dominant role over groundwater capillary rise. The ecosystem productivity becomes highly dependent on how efficiently the plants use antecedent rain water stored in shallower soils (2 m in this study) in the semiarid basins with a very deep groundwater table (>100 m), which may be decoupled from the soil water. It is apparent that the plants represented in NEW are more efficient to use antecedent rain water, resulting in a higher RUE under the same rainfall than do STATIC (Figures 8a and 8b); and the larger AI, the larger difference in RUE between NEW and STATIC (Figures 8c and 8d).
Over the Central U.S. HUC2 basins (Figure 9), the RUE difference between NEW and FD is much larger than that between NEW and STATIC ( Figure S11), suggesting that the groundwater capillary rise plays a dominant role than does the root water uptake. NEW, with a greater groundwater capillary rise (f mic = 0.6), produces a greater RUE, suggesting more efficient use of antecedent rain water stored in deeper soils or aquifers in antecedent wetter seasons, than does FD (with f mic = 0.0). However, the spatially constant f mic parameter in NEW may be overestimated for the Arkansas-White-Red, resulting in deeper-than-GRACE TWSA troughs during the 2011 and 2012 droughts, while it improves the TWSA simulation over other the Central U.S. basins ( Figure 6). It is apparent that, over the Texas-Gulf, RUE produced by NEW is still lower than that derived from the FLUXNET data during the more severe droughts of 2011 (Figure 9b), consistent with the low model biases in GPP and ET produced by NEW (Figures 6 and 7).

Discussion
In NEW, extracted by transpiration (Q T ) while replenished by root water uptake (Q R ), the plant water storage (M q ) plays an important role in modulating transpiration through its control on the opening of plants'  Figure 10f). We attribute this model deficiency to insufficient root water uptake to sustain the plant water storage and thus to meet the transpiration demand during the more severe droughts. In NEW, the root water uptake is driven by the hydraulic pressure gradient between roots (h R ) and the ambient soil (h s,i ). Like many other LSMs, Noah-MP adopts the Campbell soil water retention model (Campbell, 1974;Clapp & Hornberger, 1978;Cosby et al., 1984) to convert the model-predicted soil moisture to soil matrix water potential ( Figure 11; used in all experiments in this study). However, the Campbell model generally results large errors in the matric potential over the dry and wet ends of the soil moisture spectrum. It produces a much larger soil matric suction than does the van Genuchten (1980) model by a few orders of magnitude for almost all soil types, becoming more apparent for clay and when the soil becomes drier (Figure 11). The resulting wilting point for clay at 30 bar corresponds to 0.27 m 3 /m 3 in soil moisture from the Campbell model, whereas it is equivalent to 0.11 m 3 /m 3 from the van Genuchten model, suggesting that the van Genuchten model would allow more water to be extracted by roots (or push more water to the plants through roots). In the Texas-Gulf basin, with a large portion of clay soil (45.4%; Figure 2b), the model plants compete against the strong soil matric suction of water produced by the Campbell model used in NEW, resulting in insufficient root water uptake. Whereas in OLD, the wilting point is 0.138 m 3 /m 3 (prescribed through the look-up table of Noah-MP; used in Equation 13), which allows the plants to withdraw more water than NEW. Figure 11 also indicates that the uncertainties in the total soil water available for plant The soil water retention characteristics may be further affected by the presence of macropores in the soil matrix. Most soils contain macropores, which are formed by burrowing animals, plant roots, fissures due to freezing/thawing cycles, aggregates due to shrinking/swelling of clay soils, and soil pipes due to subsurface erosion and deposition of smaller soil particles (Beven & Germann, 1982). These macropores tend to reduce the soil matrix suction head, allow the plants to withdraw more water, and enhance infiltration at the soil surface. One may argue that they tend to enhance subsurface water flow and gravitational drainage, reducing the total amount of soil water available for plants. However, the soil macropores may be more abundant in the topsoil due to more intense mechanical, chemical, and biological weathering in the topsoil with more organic matter inputs at the surface in addition to the effects of soil erosion, deposition, and compaction. As a result, the soil permeability decays with depth, which is supported by the data of Beven and Germann (1982), Beven (1984), and Elsenbeer et al. (1992), thereby slowing down the bottom drainage, retaining the water longer for plant use, and enhancing the plants drought resilience. At present, there is still not an adequate physical theory or observational techniques to support parameterizations for use in large-scale models (Beven & Germann, 2013). The dual-domain modeling approaches with one domain representing rapid preferential flow through macrochannels and the other representing slow capillary flow through matrix (Šimůnek & van Genuchten, 2008) may represent a promising step forward for use in large-scale LSMs.
The current root water uptake model does not have a capability of predicting rooting depth, and all the above experiments assume a fixed rooting depth with N R = 4 (rooting depth of 2.0 m). An additional sensitivity experiment with N R = 3 (rooting depth of 1.0 m) indicates that the modeled GPP, transpiration, and ET averaged over CONUS drop by 17.64%, 19.73%, and 7.99%, respectively, more pronounced over water-limited basins ( Figure S12) and during droughts ( Figure S13). This indicates that rooting depth is also critical for

Journal of Advances in Modeling Earth Systems
surface water and carbon exchanges with the atmosphere. The macroscopic root growth model predicting rooting depth in HYDRUS (Hartmann et al., 2018) seems suitable for developing a dynamic rooting depth model in LSMs considering various stress factors by linking the photosynthate translocation to roots to the growth rate of rooting depth. Our layered root biomass (controlling root surface area) model should also be revised to consider other stress factors such as soil strength and aeration in addition to the soil moisture and temperature stress factors.
The new model provides a capability of representing "hydraulic redistribution" (HR; Dawson, 1993;Prieto et al., 2012) through the exchange of water between roots and the ambient soil matrix as described by Equation 2, which results in positive (water uptake from soils to roots) or negative (releasing water from roots to soils) fluxes at different depths. The resulting negative water fluxes can be regarded as HR. During dry-down periods, the roots moistened by the deep soil release water to the ambient soil in the topsoil driven by the difference between the root water head and soil matrix head, resulting in hydraulic lift ( Figure S14c). During rainfall, the roots moistened by the topsoil release water to the ambient deep soils, resulting in hydraulic descent ( Figure S14d). The modeled HR is highly dependent on the modeled root biomass ( Figure S14a), surface area density ( Figure S14b), and the soil hydraulic conductivity, which is too small to have an effect for the dry soils over the Western drylands. The present model may overestimate hydraulic lift and underestimate hydraulic descent due to the assumption of an equilibrium plant hydraulic state (neglecting elevation head). By representing the root hydrotropic growth, the present model may produce a relatively smaller amount of HR due to lesser root biomass and root surface area in the drier topsoil than a static, evenly distributed root profile. The realism of the present model in modeling HR and the effects of HR on plant water use strategy are subject to further studies (see Schymanski et al., 2008, for more discussions).   water stress factor (β) by NEW (Equation 10) and OLD over four HUC2 basins (where NEW > OLD is filled with blue, while NEW < OLD is filled with green). M q = plant water storage; Q R = root water uptake; Q T = transpiration.

10.1029/2020MS002062
Journal of Advances in Modeling Earth Systems

Summary
In a warming climate, ecosystems over drylands under persistent water stress, for example, the U.S. Southwest, and those in the dry-to-wet transitional climates under episodic water stress, for example, the Central United States, may experience longer, more frequent, and intense droughts. We developed an explicit representation of plant water storage supplied by dynamic root water uptake to improve the Noah-MP ecosystem responses to droughts. The new model explicitly represents root water uptake as a function of root surface area density at different soil layers and relates the plant stomatal water stress factor to the plant water storage available for transpiration relative to its maximum. It also allows more carbon translocation to roots when the plants are under water stress and represents root hydrotropism through more carbon translocation to roots but less turnover in wetter soil layers.
Through optimization of major leaf, root, and soil parameters against observation-based data sets averaged over the same vegetation and soil types, the model improves the prediction of LAI, GPP, and ET, and TWSA averaged over most of the 334 HUC2 basins in CONUS, except the U.S. Southwest basins facing the challenge of the strong interannual variability in the climatic forcing (e.g., the monsoon rainfall). However, the model performance is acceptable at the HUC2 basin scale over the Lower Colorado and Rio Grande basins in the U.S. Southwest. Compared to the "best" estimates of ET climatology derived from the water balance, ET wb , using the USGS HUC6 runoff climatology, the new model produces better ET climatology than FLUXNET over most of the Southern HUC6 basins including the water-limited basins in the Southwest drylands and the Central United States. However, it produced ET worse than ET wb over the northern energy-limited basins with more sandy soils, which may facilitate drainage of water. The improvement in the modeled ET during the Central U.S. droughts is further confirmed by the closer agreement of the modeled TWSA with that derived from the GRACE gravity fields.
The response of monthly RUE to monthly AI reflects how efficiently plants use antecedent rain water stored in the subsurface and thus can be regarded as a critical measure of ecosystem resiliency. The sensitivity experiments using the new model suggest that both root water uptake and groundwater capillary rise provide a "memory" for plants to use antecedent rain water and thus enhance ecosystem RUE during droughts. The explicit representation of plant water storage supplied by dynamic root water uptake enhances the plants' efficiency to use antecedent rain water stored in shallow soils over the U.S. Southwest drylands. Whereas in the Central U.S. basins, the modeled plants tend to buffer droughts through more efficient use of antecedent rain stored in deeper soils or aquifers with the aid of groundwater capillary rise.
In this study, we tested the new model's performance through data products of the aboveground carbon and water fluxes and further confirmed through the improved simulation of TWSA. Also, the new model provides a useful tool for improving our understanding of the mechanisms for plants to survive droughts. However, the dynamic root model developed in this study is still tentative and subject to further

10.1029/2020MS002062
Journal of Advances in Modeling Earth Systems developments in rooting depth dynamics by taking into consideration of more root traits and the effects of soil hydraulic properties and subject to validation against various belowground root biomass, length, surface area, and depth data. Future plant hydraulic models should not ignore the uncertainties in soil water retention models and soil macropore effects on soil water retention and infiltration.

Appendix A: Soil Matrix Resistance
Following Schymanski et al. (2008), the soil matric resistance to water flows to the root surface in the ith soil layer, Ω s,i (s), is where r R is the root radius (m), which is 2 × 10 −4 m for all plant types, and K i is hydraulic conductivity K i = K sat (θ i /θ sat ) 2b+3 (Clapp & Hornberger, 1978), where K sat is the saturated hydraulic conductivity (m s −1 ).

Appendix B: Leaf Dynamics
We have revised the schemes of carbon translocation and leaf death of the dynamic leaf model of Noah-MP (Dickinson et al., 1998;Niu et al., 2011). The leaf carbon, C leaf (g m −2 ), is where F leaf is the fraction of photosynthate translocated to leaves, L TO the rate of leaf turnover (Table 1), and R leaf the leaf respiration rate including maintenance and growth respiration.
where T death is a type-dependent constant of leaf temperature at which the death rate reaches its maximum, D T,max (Table 1), and δ T (Table 1) is a decay factor for D T,max in warmer conditions (>T death ). In this study, the water stress for leaf death is removed but considered through its effects on photosynthesis through Equation 10 or 12. LAI is then converted from C leaf using specific leaf area A leaf,s (m 2 g −1 ) ( Table 1).
The maintenance respiration rate, R leaf,m (g m −2 s −1 ), is a Q 10 (= 2.0) function of leaf temperature, T v (in°C), and leaf water content: R leaf ; m ¼ R leaf ; 25 Q 10 e Tv − 25 10 LAI R g β; where R leaf,25 (g m −2 s −1 ) is the leaf respiration rate at 25°C, a type-dependent parameter, R g is a growing season factor, and β is the plant stomatal water stress factor (Equation 10 or 12).