Disentangling the Effects of Vapor Pressure Deficit and Soil Water Availability on Canopy Conductance in a Seasonal Tropical Forest During the 2015 El Niño Drought

Water deficit in the atmosphere and soil are two key interactive factors that constrain transpiration and vegetation productivity. It is not clear which of these two factors is more important for the water and carbon flux response to drought stress in ecosystems. In this study, field data and numerical modeling were used to isolate their impact on evapotranspiration (ET) and gross primary productivity (GPP) at a tropical forest site in Barro Colorado Island (BCI), Panama, focusing on their response to the drought induced by the El Niño event of 2015–2016. Numerical simulations were performed using a plant hydrodynamic scheme (HYDRO) and a heuristic approach that ignores stomatal sensitivity to leaf water potential in the Energy Exascale Earth System Model (E3SM) Land Model (ELM). The sensitivity of canopy conductance (Gs) to vapor pressure deficit (VPD) obtained from eddy‐covariance fluxes and measured sap flux shows that, at both ecosystem and plant scale, soil water stress is more important in limiting Gs than VPD at BCI during the El Niño event. The model simulations confirmed the importance of water stress limitation on Gs, but overestimated the VPD impact on Gs compared to that estimated from the observations. We also found that the predicted soil moisture is less sensitive to the diversity of plant hydraulic traits than ET and GPP. During the dry season at BCI, seasonal ET, especially soil evaporation at VPD > 0.42 kPa, simulated using HYDRO and ELM, were too strong and will require alternative parameterizations.

Inclusion of a mechanistic representation of plant hydraulics can provide an important link between soil and evaporative demand and improve the prediction of plant response to drought (Christoffersen et al., 2016;Grossiord et al., 2020;Kennedy et al., 2019;Liu et al., 2020;Sperry & Love, 2015). Motivated by notable improvements in models that link water transport to stomatal conductance, plant hydraulics have started to be physically represented in land surface models and biosphere models, such as Ecosystem Demography Model (ED2) (Xu et al., 2016), CLM5 (Kennedy et al., 2019), a tractable land surface model (Sabot et al., 2020), the Joint UK Land Environment Simulator (JULES) (Eller et al., 2020), LM3PPA-TV (Cano et al., 2020), and Noah-MP (Li et al., 2021). Recently, a plant hydrodynamics (HYDRO) model (Christoffersen et al., 2016) has been implemented in the Functionally Assembled Terrestrial Ecosystem Simulator (FATES) (Fisher et al., 2018), which is a dynamic vegetation module that is currently coupled to CLM and the land model (ELM) of Energy Exascale Earth System Model (E3SM) (Golaz et al., 2019;Leung et al., 2020).
It has become an emerging research topic whether atmospheric moisture stress is a dominant limitation for vegetation productivity and water use during periods of hydrologic stress (Liu et al., 2020;Massmann et al., 2019;Novick et al., 2016;Sulman et al., 2016;Yuan et al., 2019) as concurrent drought (soil water stress) and high atmospheric water deficit (VPD) often co-occur . These two interactive factors both constrain transpiration and vegetation productivity and they should be accounted for together to explain changes in ecosystem fluxes (Stocker et al., 2019;and references therein). Yet few studies using land-surface models have separately resolved the response of plant functioning to high VPD and soil water stress. This can limit model predictability of future climate impacts on terrestrial ecosystems . In this study, we use a numerical model that represents plant hydraulics and observation data analysis to disentangle the importance of VPD and soil water stress on canopy conductance (G s ) at a tropical forest site in Barro Colorado Island (BCI), Panama during the 2015-2016 El Niño-Southern Oscillation (ENSO) event. The goals of this study are to improve understanding of tropical forest responses to drought and to inform experiments and observational studies for new data collection, as well as new model development needs.

Model Description
The E3SM land model version 1 (ELMv1) for this study is identical to the Community Land Model version 4.5 (Oleson et al., 2013) in most aspects except for some biogeochemistry components (Burrows et al., 2020; FANG ET AL.

10.1029/2021JD035004
2 of 20 Ricciuto et al. 2018), which are not turned on in this study. As in CLM4.5, stomatal conductance in the standard ELM version is heuristically related to soil water stress (Oleson et al., 2013): where β is the stress factor, ψ S,i is the soil water matric potential of layer i (m), ψ C is the soil water potential (m) when stomata are fully closed, ψ O is the soil water potential (m) when stomata are fully open, and r i is the root fraction in soil layer i. β ranges from 0 to 1. β = 1 when vegetation is not water stressed. This stress factor only relies on soil water potential and it does not represent plant adaptation to available water (Mencuccini et al., 2019;Sabot et al., 2020;Xu et al., 2016; and references therein).
Evapotranspiration (ET) is calculated as water extraction by roots from each soil layer, and the fraction of root uptake of water from soil layer i follows Equation 2: where r e,i is the effective root fraction. Transpiration multiplied by r e,i for each soil layer represents the local sink term (extracted from the soil) in the soil hydrology model. In this model, water transport within the roots is not simulated.
In order to overcome some of the limitations of the heuristic approach and explore mechanistically the role of plant hydraulics, we used a plant hydrodynamics scheme (HYDRO), which is an extension of the model described in Christoffersen et al. (2016). HYDRO is a standalone module that can be coupled with any host land models via a user interface. HYDRO simulates the water transport through different organs in the plants, from roots to leaves, taking into account the plant internal water storage that can buffer the imbalance of root water uptake and transpiration demand (Bartlett et al., 2019;Huang et al., 2017;N. G. Phillips et al., 2003;Pineda-Garcia et al., 2013;Tian et al., 2018). The model schematic is shown in Figure 1. An individual tree is divided into four porous media: leaf, stem, transporting root (troot), and absorbing root (aroot). Except for the absorbing root, which is divided into compartments for each soil layer, the other types are each represented by single compartments. The soil in each layer is radially discretized into cylindrical shells representing the rhizosphere around an absorbing root ( Figure 1).
The transient water mass balance equation in one dimension along the hydraulic path for each compartment i can be written as: where k r is the relative hydraulic conductance (dimensionless), K max,i is the maximum conductance (kg MPa −1 s −1 ) at the boundary, K max,aroot is the maximum conductance (kg MPa −1 s −1 ) of the absorbing root, K max,shell1 (kg MPa −1 s −1 ) is the maximum conductance of the soil shell (shell1) next to the absorbing root, k r,aroot and k r,shell1 are relative conductance (dimensionless) for the absorbing root and shell1 (nearest shell to the root), respectively. k r i 's at the compartment boundaries shared by the absorbing root and soil shell are estimated by the harmonic mean method, k r i 's at the other boundaries are estimated using the upstream weighting method, i.e., they are chosen in relation with the direction of flow (Haverkamp & Vauclin, 1979). K max,i incorporates not only the inherent plant tissue or soil conductivity per unit volume (k max,i ; kg MPa −1 m −1 s −1 ), but also sapwood area, the path length, and a correction factor for xylem conduit taper as described in Christoffersen et al. (2016).
The relationship between water potential ψ i and the relative water content s l,i (water content θ i relative to the maximum water content or porosity θ sat on a per volume basis) comprises three successive regions, i.e., a capillary region, an elastic drainage region, and embolism region. In each plant compartment, the water potential is described using pressure-volume (PV) curves given by the following equations: FANG ET AL.

10.1029/2021JD035004
4 of 20 From the top to the bottom, the three equations in Equation 6 describe three successive dehydration phases: a capillary drainage region, followed by an elastic drainage region, and then a final embolism region. These regions distinguish the locations that are responsible for the release of plant water storage (Christoffersen et al., 2016;Steppe, 2018). Capillary region represents the release of sapwood capillary water storage. In the elastic region, elastic water stored within living cells is released, and water is released in the embolism region when air replaces the water-filled volume. m cap is the slope for the linear drainage of capillary water, ψ 0 is the saturated water potential (MPa). ψ sol,i is the solute potential (or osmotic potential) (MPa) (<0) caused by the presence of solutes in living cells, ψ p,i is the pressure potential or cell wall turgor pressure (MPa), which is normally positive or zero, s ft,i is the relative water content at which elastic drainage begins, s tlp,i is the relative water content corresponding to the turgor loss point (TLP), s r,i is the residual relative water content. ψ sol and ψ p are given by: where π o (MPa) is osmotic potential at full turgor, and ε (MPa) is the bulk elastic modulus which describes the flexibility of plant tissues, which is defined as the ratio of the change in cell turgor to change in the relative leaf cell volume (Saito et al., 2006).
The xylem vulnerability curve in every compartment follows: The soil water characteristic curve is described using the Clapp & Hornberger empirical equations (Clapp & Hornberger, 1978) in Equations 10 and 11: where P 50,i is the water potential (MPa) at which 50% of maximum conductance is lost and a i is a shape index (dimensionless), λ is pore size distribution index (dimensionless), ψ b is the air entry pressure head (MPa).
where g s,max is the maximum stomatal conductance (µmol H 2 O m −2 s −1 ) in the absence of water supply limitation as given by the Ball-Berry stomatal conductance model (Collatz et al., 1991) in ELM, and β is a water stress fraction that depends on leaf water potential according to Equation 13 which replaces Equation 1, where P 50,gs is the leaf water potential at 50% stomatal closure, and a gs is shape parameter for stomatal control. Equation 13 has been implemented also in Xu et al. (2016) to calculate the stress factor, but TLP was used instead of P 50,gs . These two parameterizations (P 50,gs and TLP) of Equation 13 will be included in the following sections.
The original plant hydrodynamic model developed by Christoffersen et al. (2016) for individual trees of different sizes was modified to conform to the big leaf version of ELM, using average tree height (m) to calculate the total individual tree biomass, canopy biomass, belowground biomass for each plant functional type (PFT). The number of trees n i of the ith PFT is calculated using the following equation: where A g is the ground area (m 2 ), sla i is the specific leaf area of the ith PFT (m 2 /gC), LAI is the leaf area index (m 2 /m 2 ), wt i is the fraction of leaf area for the ith PFT, and B leaf,i is leaf biomass per individual in carbon units (gC/individual). With the exceptions of the stem and absorbing root compartments, the volume of each compartment is derived from the biomass pools calculated using the allometric equations specified in FATES as function of average height. The volume of the stem is the sapwood area times the tree height. The total volume of the absorbing root is calculated using a given specific root length (m g −1 ), root biomass (g), and absorbing root radius (m). The absorbing root volume in each soil layer is a fraction of the total root length in that layer times the total volume. When HYDRO is activated, fluxes calculated based on ET fractions in Equation 2 are replaced with the water fluxes calculated from HYDRO. It can be negative or positive depending on whether soil shells lose or gain water.
For this study, the soil in each layer is radially discretized into 5 cylindrical shells around an absorbing root ( Figure 1). The mass balance equations for the plant-soil system in each soil layer and compartment are solved simultaneously using the Newton-Raphson approach with water potentials as prognostic variables.
The coupling of ELM-HYDRO can be easily implemented through an application interface because HY-DRO is designed such that it works as a library that can be plugged into a selection of driver land models such as CLM and ELM. ELM and HYDRO share the same soil volume. Via ELM's soil hydrology scheme, runoff, drainage, and soil water fluxes between layers are calculated given infiltration as top-boundary conditions and root water extraction calculated in HYDRO as sources/sinks. The change in water content in each soil layer is then passed back to the plant hydraulics model through the application interface in ELM. A soil layer in ELM is represented by one grid cell, while it is represented by five grid cells in HYDRO. The change in ELM soil water content is downscaled to water contents of the rhizosphere shells of the same layer following a heuristic approach. That is, in the case of a net gain in soil moisture calculated in ELM relative to the previous time step, shells in HYDRO are filled with water in order from the driest to the wettest shell, and for a net loss of soil moisture, water in soil shells are drained in order from the wettest to the driest shell. For comparison, the model that ignores plant hydrodynamics is simply referred to as ELM in the following sections, and the model that considers hydrodynamics by coupling ELM with HYDRO is simply referred to as HYDRO.
FANG ET AL.

Study Sites
Our study site is located at BCI (9°10'N, 79°51'W), Panama. The 50-ha forest plot at BCI is one of the best-studied tropical forest sites in the world (Kupers et al., 2019). The wet season at BCI is roughly from May to December and the dry season is from late December to April. Meteorological data from 2003 to 2016 is available from a meteorological tower near the Lutz catchment at BCI (Faybishenko et al., 2018). BCI had a long dry season in 2016 during the strong El Niño event of 2015-16 . Annual mean precipitation during the simulation period is 2382.7 mm, while mean precipitation in the dry season is 219 mm. Average annual and dry season relative humidity at BCI are 88.5% and 83.5%. Mean annual and dry season temperatures are respectively 26.0°C and 26.4°C at BCI.

BCI Data set
Field efforts at BCI in the last 35 years have generated a rich set of data, such as continuous meteorological forcing, streamflow, ET, soil water potential, soil moisture content, as well as leaf water potential, sapflow, and plant hydraulic traits for trees at the higher elevation plateau of the island, more prone to water stress in drought conditions. This data set provides a unique opportunity to evaluate the plant hydraulics model and tropical forest response to drought conditions at this site.
At BCI, ET was measured by an eddy-covariance system installed on a tower, the AVA tower (∼1.5 km from the Lutz catchment) located 41 m above the ground on the top plateau. The eddy-covariance system comprised a sonic anemometer (CSAT3, Campbell Scientific, Logan, UT) and an open-path infrared CO 2 /H 2 O gas analyzer (LI7500, LiCOR. Lincoln, NE). High-frequency (10 Hz) measurements were acquired by a datalogger (CR1000, Campbell Scientific) and stored on a local computer. Data were processed with a custom program using a standard routine described in Detto et al. (2010). Hydrometeorological variables collected on the same tower and used in the postprocessing included: air temperate, relative humidity, air pressure, radiation components, rainfall, and soil moisture. These variables were stored at 5-min interval.
Eddy-covariance fluxes were computed on a 30-minute averaging window, a de-spiking algorithm was applied to detect and remove spikes (values greater than six standard deviations in a one-minute window), values exceeding a reasonable physical range, or low diagnostic instrument values. Fluxes were corrected for air density fluctuations (Detto & Katul, 2007). A QA/QC was applied to the 30-min turbulent fluxes to remove measurement during and immediately after rain events, out of reasonable physical limits or not satisfying other turbulent criteria (Pau et al., 2018). Gaps were filled using Artificial Neural Network and hydrometeorological inputs.
Soil moisture was monitored using three Time Domain Reflectometers (TDR, CS616, Campbell Scientific) installed vertically in the vicinity of the tower in July 2012 at 0 to 15 cm depths. The TDR probes measured the apparent permittivity of soil that can be related to the soil water content using an ad hoc calibration curve (Kelleners et al., 2005). Seven in situ gravimetric soil water content samples (0-15 cm) were collected near the probes during different soil moisture regimes (30 campaigns). TDR soil water content compares well with gravimetric sampling.
Leaf water potential was measured on four canopy trees (Inga pezizifera, Gustavia superba., Miconia argentea, and Simarouba amara) during the 2016 dry season, on March 12, 2016. These trees are relatively abundant in the footprint of the eddy covariance tower. Between predawn (0600 h) and late afternoon (1,800 h), six measurements of water potential were taken . At each time, three sun exposed leaves were collected from each tree with a pruning pole or slingshot. The leaves were immediately sealed in humidified plastic bags, placed in a cooler with ice, and, within two hours, water potential was measured with a pressure chamber (PMS Instruments Albany, OR, USA). Sap flux density was measured on each of those trees with custom-built thermal dissipation probes (Granier, 1985). On each tree, two probes were placed on opposite sides of the main bole at 1.3 m height, extending 1 cm into the sapwood. Measurements were taken every 15 s and averages were recorded every 5 min on CR1000 and CR800 data loggers (Campbell Scientific). Sap flux density was calculated using the calibration of Granier (1985).
Maximum stem area-specific hydraulic conductivity (K max ) and vulnerability curves (Equation 9) were derived from laboratory measurements on terminal branches collected from canopy trees of each of the FANG ET AL.
10.1029/2021JD035004 focal species. Following Wolfe et al. (2016), stem segments were bench-dried to stem water potentials that ranged −0.15 to −6.64 MPa, then measured for hydraulic conductivity. Stem water potential was measured with psychrometers in the laboratory following Wolfe (2017). Details can be found in the data set in Wolfe et al. (2021). Vulnerability curves were fit using nonlinear quantile regression with the package quantreg in R. Within the regression, an intercept term was included in Equation 5 to account for raw values of K max rather than assuming proportionality. K max was calculated as the intercept of the regression ( Figure S1).

Model Parameters and Sensitivity Experiments
Hydraulic trait parameters for HYDRO are listed in Table 1, which represent mean values observed for tropical forests based on a pantropical hydraulic traits data set synthesized by Christoffersen et al. (2016) from the literature and existing trait databases. We use a tree height of 24.9 m for trees observed at BCI for the simulations, which is reasonable as a global forest canopy height map in Lefsky (2010) shows that the 90th percentile height of tropical humid forests are greater than 20 m. The two parameters defining the root profile model adopted from Jackson et al. (1996) are r a = 7, r b = 1, respectively, resulting in 99% of roots in the top 4 m of soil. In ELM, the default ψ O and ψ C in Equation 1 are −66 and −255 m, respectively. The average observed LAI at BCI is ∼6.0 .
A total of 80 simulations were run using plant traits of a single species for the whole plot. Apart from the saturated xylem conductivity (K max ), water potential at which 50% of xylem conductivity is lost (P 50 ), and xylem vulnerability curve shape parameter (a) measured at BCI for the four canopy trees (Inga pezizifera, Gustavia superba., Miconia argentea, and Simarouba amara) Table 2, we also combined 29 sets of K max , P 50 , and a i in Barros et al. (2019), and 47 sets of K max , P 50 , TLP, wood density (WD, g cm −3 ) and specific leaf area (SLA, m 2 g −1 ) compiled in Xu et al. (2016) for individual species. For the Barros' data, parameter a was calculated using the sigmoidal function and the water potential at which 88% of xylem conductivity is lost (P 88 ) as reported in Barros et al. (2019). For the Xu data, a i and a gs were set to 4 and 6, respectively, following the code in the supplementary information in Xu et al. (2016). Data from Xu et al. were for seasonally dry tropical forest sites and data from Barros et al. included two tropical Amazon forests with contrasting seasonal precipitation. These data (Table S1)    the diversity of hydraulic traits on ET and gross primary productivity (GPP) response to VPD and soil water stress at BCI. As the model is the big leaf version, in which sub-grid tiling operates on the basis of plant functional types, we only consider the dominant broadleaf trees in our study. No biodiversity of trait sets is accounted for.

Calculation of G s
Canopy conductance (G s , m s −1 ) can be calculated from the eddy covariance data and sapflux measurements as shown in the following.

Inversion of the Penman-Monteith (PM) Equation
Eddy-covariance data were used to calculate G s by the following equation inverted from the PM equation ( where G s (m s-1) is canopy conductance, T is the transpiration (kg m −2 s −1 ),  (J kg −1 ) is the latent heat of water vaporization, Δ is temperature-dependent slope of the saturation-vapor pressure curve,  (Pa K −1 ) is the psychometric constant, A (W m −2 ) is the available energy of the canopy,  (kg m −3 ) is the density of dry air, C p (J K −1 kg −1 ) is the specific heat of air, e a (kPa) is the actual vapor pressure, e s is the saturation vapor pressure (kPa), VPD = e s -e a is the vapor pressure deficit, and G a (m s-1) is the aerodynamic conductance. Following the approach in Barros et al. (2019), the available energy A is assumed to be equivalent to the sum of sensible and latent heat flux. G a is calculated with the expression in Garratt (1994): where u* (m s −1 ) is the friction velocity measured at the eddy covariance height z (41 m), d (m) is the zero plane displacement height, z om (m) is the roughness length for momentum, k (-) is the von Karman's constant (0.41). Given canopy height of h = 24.9 m, d, and z om were estimated as d = 0.67 h, 0.123 h, respectively as a common practice in the absence of other information (Garratt, 1994). Ψ H is a correction for stability conditions.
Partitioning ET into transpiration and soil evaporation is not an easy task and is still an area of active research (Stoy et al., 2019). In ecosystems with high LAI, transpiration dominates ET because leaf evaporative area greatly exceeds ground area and because little energy reaches the forest floor. Therefore, T in Equation 15 was taken from latent heat flux for this analysis. G s was calculated using only the fluxes for local day times (6:00-18:00) when the canopy was dry (all data up to 12 h after precipitation were removed) to exclude evaporation from intercepted rainfall and minimize contribution from soil evaporation (Barros et al., 2019).

Estimated Gs Based on Sap Flux
Eddy covariance water flux measurements combine transpiration and evaporation, therefore the VPD effect on G s based on the inversion method above represents a whole ecosystem level response which includes contribution from soil (Scanlon & Kustas, 2012;Wilson et al., 2001). At BCI, LAI varies between 5 and 6.5 with lower values during the dry season Wirth et al., 2001). Thus, considering the large foliage evaporative area relative to the ground surface and that little energy reaches the forest floor, except for periods immediately after rainfall, the contribution of soil evaporation should be small.
In order to verify the effect of VPD on G s without having to worry about soil evaporation, G s was also estimated using sap flux measurements from individual trees (Ewers et al., 2001;Granier, 1987;Oren et al., 1999 FANG ET AL.

Empirical Relationship Between G s and VPD Under Varying Soil Moisture Conditions
The stomatal sensitivity is proportional to the magnitude of G s when VPD is less than or equal to 1 kPa (Oren et al., 1999 and references therein). Restricting the analysis to high solar radiation (R s > 500 W m −2 ) to minimize light effect on G s , the relationship between G S and VPD changes with soil moisture and is quantified using the following linear regression equation Oren et al., 1999): where G s, ref (m s −1 ) is a reference canopy conductance at VPD = 1 kPa. Parameter m describes the sensitivity of G s to VPD.
To estimate the parameters in Equation 20, the hourly soil moisture data during the dry season in 2016 were first binned into quantiles defined by the 0-15th, 15th-30th, 30th-50th, 50th-70th, 70th-90th, and 90th-100th percentiles. The parameters were then determined by linear regression of G s with the observed ln(VPD) following the procedure in Novick et al. (2016). To estimate the magnitude of limitations to G s due to soil water and VPD, we assumed the well-watered reference canopy conductance (G s,ww ) to be the average of conductance that meet two conditions: (1) soil moisture content exceeds the 90th percentile, and (2) 0.9 < VPD < 1.1 kPa during 2016 for the eddy covariance data analysis as the estimated G s may be influenced by soil evaporation. For the other approaches, G s,ww was assumed to be the maximum of those in 0.9 < VPD < 1.1 kPa. The total limitation on G s is calculated by subtracting each G s predicted by Equation 20 from G s,ww . The fraction (α VPD ) that limits G s due to VPD is m×ln(VPD) divided by the total G s limitation. The mean parameters were then reported for each bin. These processes were done for G s estimated from the measurements and that calculated internally from the simulations.

Seasonal Model Performance
At BCI, the monthly averaged VPD during 2013-2017 ranges from 0.25 to 0.56 kPa. The ET and GPP correlation between the simulation and the observation using HYDRO was low compared to the simulation using the heuristic approach. The root-mean-squared error (  This could be due to the lack of representation of soil macroporosity, which is important in tropical forests (Bodner et al., 2013).
Plotted against monthly VPD using the results from 2013 to 2016, there is more than 10% variation in ET, transpiration (T), soil evaporation (E) and GPP (  Xu et al., 2016). It is referred to as "Sim37" hereafter. In Sim37, K max is around the median, and P 50 and TLP are the least negative among the ensemble of parameters (Figure 2). From Sim37 we can see that the simulated transpiration and GPP from ELM and HYDRO are almost identical ( Figure 5). However, soil evaporation increases with the increase of VPD, and it doubles that from ELM at VPD = 0.6 kPa (Figure 5c). In general transpiration from model simulation plateaus while GPP decreases when VPD > 0.4 kPa. Relative to ELM, the larger simulated ET from HYDRO ( Figure 5a) is mainly due to the larger soil evaporation (Figure 5c) from the wetter soil (Figure 4). During the dry season, the ET simulated by ELM and HYDRO are both too strong compared to the observations while the simulated transpiration alone matches the observed ET very well (note the circles in Figure 5b are observed ET, lines are T from simulations), suggesting that the simulated soil evaporation may be excessive.

Diurnal Leaf Water Potential and Sap Flux
As described in Section 2.3, leaf water potential at BCI was measured only on March 12, 2016. There were also continuous sapflow measurements during 2016. The averaged diurnal cycle of observations of leaf water potential and sapflow (averaged from four species at the site) were compared to the simulated variables FANG ET AL.

10.1029/2021JD035004
13 of 20 as shown in Figure 6. In general, the diurnal pattern of the simulations follows that of the observations. Compared to the observations: (1) the diurnal variation of leaf water potential from a majority of the HY-DRO simulations is less negative (Figure 6a); and (2) the simulated sapflow is too high (Figure 6b). None of the members of the HYDRO ensemble of simulations can fully capture both the measured leaf water potential and sap flux. We note that the sapflow measurements can depend on many factors including the point where the probe is inserted if the xylem structure is irregular, the depth inside the xylem, and the reaction of the plant to the wound. Although the absolute magnitude of the flux can diverge substantially among different probes, the temporal pattern, as also predicted by most of the model simulations, remains robust.
The RMSE of leaf water potential ranges from 0.25 to 1.04 MPa (Figure 6c) and the RMSE of sap flux ranges from 1.13 to 18.73 cm h −1 . The R 2 between the observed LWP and the HYDRO simulations (0.36-0.60) is lower than the R 2 between the measured sap flux and the Hydro simulations which ranges from 0.88 to 0.99. Of the ensemble of HYDRO simulations, a simulation using one set of parameters (K max = 2.0 kg MPa −1 m −1 s −1 , P 50 = −2.12 MPa, a = 4, TLP = −1.76 MPa, SA = 0.005 m 2 g −1 , and WD = 0.52 g cm −3 ) gives the best model performance for sap flux (RMSE = 1.13 cm h −1 ) (the red line in Figure 6b). This simulation is referred to as "Sim58" hereafter. Even though one of the simulation (Sim60) with the lowest hydraulic conductance ( potential (RMSE = 0.25 MPa for Sim60, the green line Figure 6a), the simulated sap flux in the afternoon is not well correlated with the measurement. Neither of these two simulations can explain the seasonal variations of ET and GPP and the K max from Sim58 and Sim60 is smaller than that (3.96 kg MPa −1 m −1 s −1 ) in Sim37 that best reproduced the observed ET and GPP. For ET, T, E, and GPP shown in Figure 5, Sim58 and Sim60 are at the high end of the shaded gray area (not plotted).

VPD Importance on G s Limitation
Here the G s limitation is estimated using the two methods described in Section 2.5.1 and 2.5.2 to evaluate the role of VPD in plant response during the 2015-2016 El Niño drought. The simulated G s limitation is also analyzed and compared to the observation-based estimates. Using Equation 15, hourly G s was first calculated by inverting the Penman-Monteith equation using the eddy covariance (EC) data in the dry season. The soil water content derived from the TDR measurements was binned to derive the parameters for the G s versus VPD relationship (the purple line in Figure 7). The average parameters for the G s ∼ VPD relationship is G s,ref = 8.3 mm/s and m = 0.79. The average ratio of VPD limitation to the total limitation for G s (α VPD ) was calculated as 0.1. There is no clear trend of these parameters varying with soil moisture conditions. On the other hand, the parameters derived from G s calculated using the sap flux from two species (Gustavia superba, and Simarouba amara) exhibits smaller G s,ref (average of 0.55 mm s −1 ) and m (average of 0.54) higher α VPD (average of 0.23). The G s estimated from the sap flux of Inga Pezizifera and Miconia argenea were not accounted for as the former had a lot of gaps in the data set and the latter showed no VPD limitation.
If not stated otherwise, the following comparisons are between the HYDRO simulations with better performances (Sim37, Sim58, and Sim60) and ELM. The hourly G s data were computed internally in the model. Even though ELM simulates drier soil (black line in Figure 7) compared to the HYDRO simulations (gray lines and Sim37, Sim58, and Sim 60 highlighted in blue, red, and green in Figure 7), the HDRO simulation and the ELM simulation show decrease of G s,ref with soil moisture and the average G s,ref are comparable. Within the range of observed soil moisture, calculated average α VPD for Sim37, Sim58, Sim60, and ELM are 0.3, 0.36, 0.32, and 0.24, respectively. These fractions from HYDRO are higher than those (α VPD < 0.23) calculated based on measurements, i.e., VPD has higher importance constraining G s from model simulations than observations. Estimated parameter m and α VPD from Sim37 and Sim60 are within the range of those from the observations when soil moisture is less than 0.28. α VPD from ELM is comparable to those from the observations, but m is lower. VPD tends to be a dominant limiting driver for simulations with more negative P 50 and TLP. For example, the simulation with parameters K max = 3.03 (0.25-0.30), α VPD is greater than 0.5, VPD as a limiting driver to G s becomes more dominant. Among all simulations, parameter G s,ref is almost the lowest while m is almost the highest for Sim37 (highlighted in blue in Figure 7) that best explains the seasonal variation. The regression parameters estimated from Sim37 are comparable to those from HYDRO Sim60 (highlighted in green in Figure 7) that best explains the leaf water potential. At low soil moisture (VWC < 0.28), almost all HYDRO simulations have higher α VPD compared to that estimated from the observations (Figure 7c). HYDRO Sim58 that matches the sap flux measurement gives a high average α VPD (0.40).

Discussion
As the forest plot at BCI has been intensively studied, many measurements can be used for model parameterization development and evaluation. Using analysis from observations and numerical model simulations, we evaluated the relative importance of VPD as a driver for G s limitation at BCI during the ENSO event in 2015. From the regression model derived based on observations, we found that VPD is not a dominant driver for G s limitation at BCI. However, variations of the regression parameters derived from observations using different methods are large. Similarly, variation in the regression parameters from model simulations with different model parameter values is large. VPD is a more dominant driver for G s limitation from the simulations than observations.

Estimated G s and VPD Relationship From Different Observations
EC data and sap flux measurements have been commonly used to estimate G s . When we used EC data to determine the parameter m that indicates G s sensitivity to VPD, it was higher than the value determined from the sap flux measurements (average = 0.79 and 0.54, respectively). Despite the uncertainty in estimating G s and m, the VPD contribution to G s limitation is not dominant (α VPD < 0.4).
We chose the EC data when the canopy was dry (all data up to 12 h after precipitation were removed) to assume that most of the flux from EC was due to transpiration. However, sap flux (representing tree level transpiration) can only explain 60%-70% of the ET subset, i.e., not all ET from the EC data is due to transpiration only. We found the correlation between the hourly sap flux and soil water content of the top 15 cm was low (R 2 = 0.01) and thus they were almost decoupled emphasizing the role of water demand (VPD) on sap flow. Nevertheless, the measured soil water content can explain more than 10% of ET when combined with sap flux. These results suggest that the influence of soil evaporation was not completely removed when EC data were used to calculate G s .
The average G s,ref from the EC data and sap flux are 8.3, and 0.55 mm s −1 , respectively. They are lower than the reported 11.8 mm/s in Granier et al. (1996) for tropical rainforest and the monthly mean G s,ref (∼10 mm/s) in Barros et al. (2019). When EC data were used to calculate G s , we assumed the available energy absorbed by the surface to be equivalent to the sum of sensible and latent heat flux because of the lack of measured soil heat flux. This assumption is reasonable in that it accounts for problems with flux underestimation and energy closure that are known to affect turbulent fluxes obtained from eddy covariance. Our analysis using net radiation minus soil heat flux (approximated from the GradCal approach in Liebethal et al. (2005)) as available energy in the PM equation shows underestimated G s and lower α VPD ( Figure S2) compared to those using the EC fluxes. The leaf area to sapwood area ratio needs to be better quantified in order to better estimate G s using sap flux measurements from individual trees.

Simulated G s and VPD Relationship
As inferred from the G s,ref simulated by ELM and HYDRO (Figure 7a), there is no significant difference in the magnitude of G s,ref calculated by the model whether the stomatal sensitivity to leaf water potential is ignored or not. But the range of G s,ref is much lower than that from the EC data (by a factor of ∼5.1 on average). When soil moisture is less than 0.28 (corresponding soil water potential of −0.09 MPa) and within the range of the observations, the sensitivity of G s to VPD (parameter m) is relatively low using the heuristic approach compared to those estimated from the observations and most of the simulations using HYDRO. Depending on the metrics for model performance, estimated parameter m can differ significantly. m from Sim58 that better explains the measured sap flux is almost the lowest among all HYDRO simulations.
There is not a common set of hydraulic traits that can be used to explain the seasonal and diurnal measured variables. Regardless of which variables can be better explained by model simulations, VPD contribution to G s limitation from HYDRO is relatively high compared to the estimates from the measurements. However, most simulations also indicate VPD limitation is not as important as soil water stress (α VPD < 0.5).

Model Sensitivity to the Diversity of Hydraulic Traits
Monthly ET and GPP are most sensitive to plant hydraulic traits in March and April, with ET ranged from 2.6 to 3.4 mm d −1 and GPP from 6.2 to 8.8 gC m −2 d −1 . Midday leaf water potential and sap flux on March 12, 2016 can vary from −1.5 to −0.75 MPa and 15 to 52 cm h −1 , respectively. They impact the stomatal sensitivity to VPD (average m is between 0.23 and 0.61). However, their impact on soil water content is insignificant compared to that from different model representations.

Model Performance and Future Work
Measured ET, soil moisture, leaf water potential, sap flux at BCI were used to evaluate the model performance. Compared to the observations, monthly ET were underestimated and GPP were overestimated in the wet season. In the dry season, all model simulations overestimated ET and a majority of the model simulations overestimated GPP. Simulated transpiration was higher or equal to measured ET in the dry season, which suggests there is a need for parameterization of stomatal conductance in the model. Soil evaporation simulated by HYDRO is too strong, which exacerbates total ET compared to the observation. This is no surprise as HYDRO simulated wetter soil and the soil evaporation in ELM and CLM4.5 is dependent on the top layer volumetric soil moisture (Swenson & Lawrence, 2014). A better parameterization of soil evaporation (e.g., an alternative soil resistance parameterization in Swenson & Lawrence, 2014) can be tested in the future.
Our results showed that soil moisture can be better simulated when plant hydraulic processes and hydraulic redistribution are considered in the model. The fact that none of the model simulations with different parameter values was able to capture the observed high-water content in the wet season indicates missing processes in the model, such as the change of soil characteristics due to wetting-drying cycles.
The ensemble simulations did not show a single set of hydraulic traits that can produce results matching well with the observed diurnal leaf water potential, sap flow, and seasonal ET and GPP. This could be due to the footprint spatial-temporal variability of ET from the EC tower that cannot be accounted for in the model. Another cause is that the spatial-temporal heterogeneity of soil moisture experienced by each individual tree, which is on the scale of meters, differs from that of the plot scale. As demonstrated in Baker et al. (2017), model ET could be improved if soil moisture spatial variability was accounted for in land models. This can be addressed in the future version of ELM which accounts for subgrid topographic landunits and biodiversity of trait sets interacting with other abiotic drivers.

Conclusion
A mechanistic plant hydraulics process model has been included in the "big-leaf" version of ELM. Observations and model simulations are used to understand the role of VPD as a driver for Gs limitation in a tropical forest site at BCI. Including plant hydraulics in land models has a noticeable effect on simulating the ET and GPP, and the average depth of plant water uptake from soil, particularly under drought conditions. Although increasing model complexity by representing plant hydraulics does not noticeably improve model performance in ET compared to the observations and ELM, the model can better explain where water is being extracted from by the roots to meet transpiration demand and allow more informed interrogation of the hydraulic process throughout the soil-plant-atmosphere continuum.
In the long dry season of 2016 caused by the strong El Niño event, G s limitation of BCI plants was driven more by water stress than VPD, as indicated by the G s estimates from the site eddy covariance and sap flow measurements. Model simulations also indicate a larger role of water stress than VPD in the G s limitation, although generally the model indicates a more important role of VPD than suggested by the observations.
Despite the available data set at BCI and improvement in modeled soil moisture, model-data comparison indicates a need for more experimental and observational studies. During the dry season at BCI, seasonal ET, especially soil evaporation at VPD > 0.42 kPa, simulated using HYDRO and ELM are too strong. Both ET and GPP are more sensitive to the diversity of plant hydraulic traits than soil moisture. The fact that none of the model simulations with a wide range of model parameter values are able to match different measured variables simultaneously suggest model structural uncertainty in addition to parameter uncertainty. Hence, model-data comparison at BCI points to the need for improving parameterizations, particularly for stomatal conductance and soil evaporation, in ELM. More observation data are also needed to better constrain the plant hydraulic traits, which have significant impact on simulating ET, GPP, and the role of VPD in the G s limitation. Furthermore, the model runs were for a single plant functional type. Biodiversity model such as the FATES demographic model should be used for the evaluations in the future.

Data Availability Statement
Observational data and data for the model are available at http://doi.org/10.5281/zenodo.3752127.