Hydraulically‐vulnerable trees survive on deep‐water access during droughts in a tropical forest

Summary Deep‐water access is arguably the most effective, but under‐studied, mechanism that plants employ to survive during drought. Vulnerability to embolism and hydraulic safety margins can predict mortality risk at given levels of dehydration, but deep‐water access may delay plant dehydration. Here, we tested the role of deep‐water access in enabling survival within a diverse tropical forest community in Panama using a novel data‐model approach. We inversely estimated the effective rooting depth (ERD, as the average depth of water extraction), for 29 canopy species by linking diameter growth dynamics (1990–2015) to vapor pressure deficit, water potentials in the whole‐soil column, and leaf hydraulic vulnerability curves. We validated ERD estimates against existing isotopic data of potential water‐access depths. Across species, deeper ERD was associated with higher maximum stem hydraulic conductivity, greater vulnerability to xylem embolism, narrower safety margins, and lower mortality rates during extreme droughts over 35 years (1981–2015) among evergreen species. Species exposure to water stress declined with deeper ERD indicating that trees compensate for water stress‐related mortality risk through deep‐water access. The role of deep‐water access in mitigating mortality of hydraulically‐vulnerable trees has important implications for our predictive understanding of forest dynamics under current and future climates.


Article acceptance date: 29 April 2021
The following Supporting Information is available for this article: Figure S1 Observed relationship between GPP and VPD at Barro Colorado Island. Figure S2 Seasonal variation in estimated standardised LAI for 29 species at Barro Colorado Island. Figure S3 Ensembles of soil hydraulic conductivity variation by depth for Conrad Trail Stream Catchment, Barro Colorado Island. Data from Godsey et al. (2004). Figure S4 Volumetric water content vs. TDR data for a vertical TDR probe used for TDR probe calibration. Figure S5 δ 2 H of soil water by depth measured at Barro Colorado Island. Data from Meinzer et al. (1999). Figure S6 Parameter sensitivity of ELM-FATES soil water content and evapo-transpiration. Figure S7 Parameter sensitivity of ELM-FATES stream discharge. Figure S8 ELM-FATES predicted soil water potential dynamics by depth over 1990-2018. Figure S9 Leaf hydraulic vulnerability curves fitted to observed " leaf vs. Ψleaf data for 21 tree species from Barro Colorado Island. Figure S10 Percent loss of leaf hydraulic conductivity curves for ERD species, including curves fitted to " leaf vs. Ψleaf as well as those based on scaling relationships with leaf mass area and wood specific gravity. Figure S11 Effective rooting depth (ERD) versus δ 2 Hxylem for each of the ERD model structures tested. Figure S12 Correlation matrix between ERD and hydraulic traits. Figure S13 Mortality rates for deciduous species across census intervals from 1985-2015 versus ERD. Figure S14 Modeled effective rooting depths of deciduous species versus time spent beyond critical hydraulic threshold. Table S1 ELM-FATES parameters used to generate ensembles, with description, prescribed global ranges, rationale for the choice of ranges and references, as well as ranges for best-fit ensembles. Table S2 QA/QC procedure applied to eddy covariance fluxes. Table S3 Above-ground hydraulic traits data from Wolfe et al. (2019) and Wolfe et al. (2021) used for comparison with ERD. Table S4 Leaf vulnerability curve parameters A & B, " max,leaf and Ψ 20,leaf by source (data or model) for the 29 species with ERD estimates.
Methods S1 Alternative structures for effective rooting depth model. Methods S2 Statistics for identifying best-fit ERD. Methods S3 Processing of forest data for growth estimates. Methods S4 Leaf Area Index calculations. Methods S5 Details of the ELM-FATES model. Methods S6 Details for ELM-FATES model parameterization.

Methods S7 ELM-FATES calibration.
Dataset S1 Microclimatic and flux tower data. Dataset S2 Leaf hydraulic conductivity and vulnerability to cavitation. Dataset S3 Wood specific gravity (WSG) and leaf mass area (LMA). Dataset S4 Leaf deciduousness categories. Dataset S5 Stream discharge from the Conrad catchment. Dataset S6 Volumetric water content. Dataset S7 Stem maximum hydraulic conductivity and vulnerability to cavitation. Dataset S8 Leaf turgor loss point. Dataset S9 Above-ground hydraulic safety margins. Notes S1 ERD model structure selection. Notes S2 Additional results for exposure to water-stress.

Figure S1
Scatterplot of mean midday (11:00 am to 15:00 pm) gross primary productivity (GPP) and mean midday vapor pressure deficit (VPD) measured at the eddy flux tower in Barro Colorado Island over 2012-2017. Red line shows the fit from a polynomial model described in inset.

Figure S2
Seasonal variation in standardised LAI (range 0-1) for 29 species at Barro Colorado Island (panels). Species are color-coded by leaf habit as defined by expert botanists.

Figure S3
Generated 5000 curves of soil hydraulic conductivity variation by depth based on Latin Hypercube Sampling within 95% CI of observed data at depths of 12.5 cm and 60 cm for Conrad Trail Stream Catchment (Godsey et al., 2004), assuming a linear decline between those depths, and values before 12.5 and after 60 cm as that of 12.5 cm and 60 cm, respectively. Only 100 curves are shown for clarity. Note the log scale on horizontal axis.

Figure S4
Volumetric water content vs. TDR data for a vertical TDR probe (0-15 cm depth) at the base of the eddy flux tower near the 50-ha plot at BCI.

Figure S6
Parameter sensitivity of ELM-FATES soil water content (that is, VWC, a-b) and evapotranspiration (ET, c-d) is shown for two representative sample months, as sensitivity for a particular parameter varies depending on the chosen month. Horizontal axes show global ranges for a given parameter. Parameters are color coded and are consistent across panels. Each point represents a simulation among the 5000 ensemble-member runs. See Table S1 for description, global ranges and chosen ranges for parameters.

Figure S7
Parameter sensitivity of ELM-FATES stream discharge (QRUNOFF) is shown for four sample months (a-d), as sensitivity for a particular parameter varies depending on the chosen month. Horizontal axes show global ranges for a given parameter. Parameters are color coded and are consistent across panels. Each point represents a simulation among the 5000 ensemble-member runs. See Table S1 for description, global ranges and chosen ranges for parameters.

Figure S8
Daily mean soil water potential (& soil,# ) by depth ' (in m; panels) over 1990-2018 predicted from each of the 100 best-fit ensembles (colored lines). Vertical gray lines demarcate median dates of the tree demographic censuses at the BCI 50 ha plot.

Figure S9
Leaf hydraulic vulnerability curves fitted to observed " leaf vs. Ψleaf data for 21 tree species from BCI-in absolute terms (a, b), as well as, in terms of percent loss of conductivity (c, d). Each line represents a species, color coded by wood specific gravity (WSG) (a, c), or leaf mass area (LMA) (b, d), with blue tones indicating smaller values of WSG or LMA and red tones indicating larger values.

Figure S10
Percent loss of leaf hydraulic conductivity curves for species with ERD estimates. Each line represents a species, color coded by source of the curve, with those fitted to observed data for " leaf vs. Ψleaf in yellow (data) and those based on the scaling relationships with leaf mass area and wood specific gravity in gray (model).

Figure S11
Effective rooting depth versus δ ! $xylem for each of the ERD model structures tested (panels). For each ERD model structure only those species have an estimate of ERD that have passed the test of Pearson's ( > 0.71 during ERD parameterisation. The number of species with an ERD estimate as well as δ ! $ for each model structure are shown in inset, along with the goodness-of-fit () 2 ) and *-values for a linear relationship between ERD and δ ! $xylem. ERD models in panels a-e correspond to Eq. (S1)-(S6), respectively.

Figure S13
Mortality rates across census intervals from 1985-2015 for species (circles) of various deciduous leaf-habit categories (colors) plotted against ERD. ) 2 and significance levels for linear model fits are given in panel insets. Census-intervals with a superscript * and * * are shown for comparison with a corresponding figure for evergreen species (Fig. 8 in the main text) and indicate intervals in which significant mortality was explained by ERD among evergreen species at -= 0.1 and -= 0.05, respectively.

Figure S14
Modeled effective rooting depth (ERD; horizontal-axis) versus time spent beyond critical hydraulic threshold (vertical axis) by census interval (colored bars) for deciduous species in the mortality analyses (Fig. S13). Each bar represents the average time species of the same ERD spent beyond species-specific critical hydraulic thresholds in a given interval, that is, the proportion of days for which Ψ soil,#$ERD was more negative than species Ψcrit, defined as Ψ20,leaf, and where ' is the soil depth matching species ERD. Standard error of the mean are shown over each bar. Note the squared y-axis scale.  (Zeng, 2001) that regulates the shape of the rooting profile. Range corresponds to this parameter specified for BATS (or IGBP) land cover classified Deciduous Broadleaf Trees (5.9 m -1 and Evergreen Broadleaf Trees (7.4 m -1 ) as given in Table 1 of (Zeng, 2001 (Zeng, 2001) that regulates the depth of the rooting profile. Chosen parameter range of b is derived using this equation so as to fit the observed range of rooting depth (dr) of 2 -18 m for Tropical Deciduous Forest (mean ± Standard Error (SE); 3.7 ± 0.5, n = 5 trees; min = 2, max = 4.7) and Tropical Evergreen Forest (mean ± S.E.; 7.3 ± 0.5, n = 3 trees and 3 communities; min = 2, max = 18 m) combined (Canadell et al., 1996). Besides the direct observation of roots at 18 m included by (Nepstad et al., 1994) in Paragominas, eastern Amazonia that is included in the above study; in Tapajos, eastern Amazonia water extraction by roots was also inferred up to 18 m. (Davidson et al., 2011) fates_smpsc For reference, a storm with 12.5 mm hr -1 rainfall intensity has a 0.2 probability of occurring in any given rainfall event.
fpi_max Maximum interception fraction of precipitation 0.05 0.44 0.06 0.43 unitless Based on throughfall data from (Zimmermann et al., 2010) and precipitation data for BCI from STRI Physical Monitoring program, defined for precipitation events greater than 10 mm.  Below we present all the alternative models for effective rooting depth (ERD). Daily average growth, ' ( !,#|% for species ), in the census interval * is described as follows, in which & leaf , and thus +,leaf , is driven by soil water dynamics at depth . for a given hydrological realisation ℎ: where nt are the total number of days in census interval *, the superscript * indicates that the variable has been standardized between 0 and 1, and 7 is the model error term. VPD ; is the predicted GPP from a locally derived polynomial relationship between gross primary productivity (GPP) and VPD, so as to account for the non-linear impact of VPD on growth. See the ERD model description in the main text for other details.

Methods S2 Statistics for identifying best-fit ERD.
For each ERD model structure (Eqs. (S1)-(S6)), species-specific best-fit ERD was identified as follows: for example, for the chosen model (Eq. (S3)), for a species ) and hydrological realization h for soil water potential Ψ soil,% dynamics, we obtained Akaike's Information Criterion (AIC) and Goodness-of-fit (A 2 ) for each modeled soil layer depth .. We used AIC as the first selection criterion to ensure that the fitted growth values from Eq. (S3) are within the confidence interval of the error distribution of observed growth (as the same model Eq. (S3) was being compared for each h, number of model parameters were the same). For each h, we retained those depths with the maximum likelihood of growth (similar AIC based support, see below) filtering out the rest. However, AIC does not ensure that fitted values will have a positive trend with growth. We therefore used A 2 between observed and fitted growth values as a goodness-of-fit measure and retained depths associated with a relatively high threshold A 2 value of 0.5; after filtering out depths associated with negative correlations B, if any. Finally, among these we chose the depth with the highest A 2 per h. Results were not sensitive to inclusion of this final step of retaining only a single depth as the best-fit (not shown). For a species ), we defined effective rooting depth (ERD) as the median of the best-fit depth over all h, which can be described as: where r is the correlation, and A 4 the Goodness-of-fit, between observed growth ' and estimated growth ' ( . ℒ is the likelihood function, 0 N is the maximum-likelihood estimate, ' the observed growth with P observations (or, census intervals *). Q = (1, .., 100), the 100 best-fit hydrological realizations and R = (0.4, 1, 1.7, 2.9, 4.7, 7.8, 13) m. We removed finely resolved soil layers <0.21 m and averaged Ψ soil,% for depths 0.21, 0.37 and 0.62 (=0.4) m as sufficient resolution for model validation.

Methods S3
Processing of forest data for growth estimates.
We used tree diameter recensus data from the Barro Colorado 50-ha plot (Condit et al., 2019). We selected trees for which the height of the diameter measurements were unchanged across all censuses (1990,1995,2000,2005,2010,2015). For each individual tree (S), we calculated absolute growth rates ' T ),# (cm yr -1 ) for each of the five census intervals * as the difference between tree dbh measured in successive censuses prorated by the inter-census interval based on the exact date of measurements. We eliminated positive outliers as any growth rates amounting to >75 mm yr -1 , as this was the highest rate observed with confidence for the fastest growing species (Condit, 2017). Negative outliers were removed based on known measurement error, U 6 = 0.006214V + 0.903, where d is dbh. Any stem whose later dbh measurement fell <4U 6 the earlier measurement was removed (Condit, 2017). As tree growth rates vary with size (Muller-Landau et al., 2006), to keep growth rates comparable across census intervals, we obtained residuals R ),# from a dbh model of tree growth based on a tree's dbh at the beginning of a census interval (R ),# = ' T ),# − W(V ),# )). The diameter model of growth, W(V ),# ), a B-spline based polynomial equation with five degrees of freedom, was fitted to growth rates from all trees pooled together as there was insufficient species-specific data for large trees (A 2 = 0.27; X-value < 0.01). We retained only those trees that had R ),# estimates for all census intervals (total trees = 972) and standardized (centered and scaled) individual tree time series of R ),# . For each species ) we then obtained growth time series ' ! as the time-series of medians of R ),# for the five census intervals * to use in the ERD model. ' ! were obtained for only those species (4 = 29) with complete records on at least three trees (median 10, max 111 trees per species).

Methods S4 Leaf Area Index calculations.
We generated an estimate of the mean seasonality curve for normalized leaf area index (LAI * , range 0-1, unitless) for individual species as follows. We interpolated species-specific weekly leaf-fall data collected at two similar sites in BCI (1985-2020) to a daily estimate of leaf mass fall per unit ground area (Wright & Cornejo, 1990). We estimated living canopy leaf mass per unit ground area present at time t, L[t], from mean leaf longevity in days, LL, and leaf mass fall per unit ground area at time t
For 13 species, an estimate of leaf lifespan was available from direct observations at a dry and wet site in Panama (Osnas et al., 2018). For the remaining species leaf lifespan was gapfilled with the following preference order depending on data availability: genus-level mean (n = 1), family-level mean (n = 7), site-level mean (n = 6 deciduous species using data from the site Parque Natural Metropolitano situated on the drier side of the rainfall gradient along the Isthmus of Panama). For one species, leaf lifespan was predicted from leaf mass area (LMA) (YSWZ)X[4 = 3.18 + 2.72,\]; A 2 = 0.65; X = 0.0027, 4 = 11).
Although deciduous leaf-habits are generally expected to be fully deciduous, at BCI there is a large within-species variation in the timing of leaf-drop and leaf-out and not all individuals of a deciduous species reach full deciduousness each year (Condit et al., 2000). Our reconstruction of LAI * seasonality is based upon species leaf-fall data, which captures the within-species variation in the timing of leaf-drop. So, in our reconstruction the seasonal decline in LAI * appears to be more spread out and does not reach zero than that is one expected for an idealized individual representing each of the deciduous leaf habits. Greater uncertainty in leaflifespan for deciduous species may have compounded this issue further. We did not normalize the deciduous species' LAI * curves to reach zero as it may further distort the shape of the LAI * curve shifting a larger period under a lower LAI * .

Methods S5 Details of the ELM-FATES model.
In ELM, soil water is predicted from a multi-layer model, in which the vertical soil moisture transport is governed by infiltration, surface and sub-surface runoff, gradient diffusion, gravity, canopy transpiration through root extraction, and interactions with groundwater (see Oleson et al. (2013) for process details and equations). FATES simulates size-structured tree community complete with individual tree rooting profiles. In our use of reduced-complexity configuration of ELM-FATES with a single plant functional type, ELM removes total community canopy transpiration as water-uptake from the soil layers in proportion to the total root fraction in each layer in the 1-D column. Note that we use the soil water potential dynamics obtained from ELM-FATES as an ecosystem-level outcome to drive our species-level, empirical ERD model. We thus believe that there is no conflict/circularity between tree-level rooting profiles or water uptake processes in ELM-FATES and the ERD model.
We ran ELM with FATES vegetation, in which the ELM model simulates interception, throughfall, canopy drip, infiltration, evaporation, surface runoff, subsurface drainage, redistribution within the soil column, and groundwater discharge and recharge so as to simulate changes in canopy water Δ^c an , surface water Δ^s fc , soil water Δ^l iq,% " at depth . ) and water in an unconfined aquifer Δ^7 (omitting processes relevant to snow, wetlands or lakes). Conservation of these terms is expressed as follows: Δ^c an + Δ^s fc + 5 Δ where, `r ain is rainfall, a 9 is transpiration of the forest, a : is ground evaporation, `o ver is surface runoff, `h 2osfc is runoff from surface water storage, `drai is subsurface drainage (all in mm H2O) and Δ* is the time-step (here, 1 hr).
Justification for the use of static stand structure: in static stand structure model, the model is started from an observed forest inventory and recruitment, growth and mortality are deactivated. Holding the stand structure static rather than allowing it to emerge from physiologically-mediated competition between trees removes a set of internal model feedbacks and thus the only differences between model simulations are driven by direct physiological effects of perturbed model parameters. This allows us to focus on the hydraulic and hydrological questions explored in this paper.

Methods S6 Details for ELM-FATES model parameterization.
In the Conrad Trail Stream catchment, in which the 50 ha plot is situated, soil hydraulic conductivity decreases with depth from 12.5 and 60 cm (Godsey et al. (2004); Table 1). We generated 5000 LHS samples within the 95% CI for the two depths and assumed a linear decline between each paired sample, while for depths <12.5 cm and >60 cm, we assumed the same value as that for depths 12.5 and 60 cm, respectively (Fig. S3).
For greater accuracy, instead of the ELM default pedo-transfer functions, we estimated parameters of soil retention curves using existing data for gravimetric water content (GWC) vs. Ψsoil for depths of 0.15, 0.4 and 1 m (Kupers et al., 2019). We converted GWC to volumetric (VWC, c) using a measured bulk density value of 0.8 g cm -3 (see below). The Campbell (1974) empirical equation was fitted for all depths pooled together, as exploratory analyses indicated that differences in water retention among depths were nominal, where Ψ 6 is the air-entry matric potential, θ ! is the saturated volumetric water content, and g, an index for soil pore-size distribution, respectively. Ψ 6 , θ ! and g that best-fit the data were 200 (mm H2O), 0.51 (cm 3 cm -3 ) and 10 (unitless), respectively.

Methods S7 ELM-FATES calibration.
We calibrated ELM-FATES over 2012-2018 against three key fluxes and states in Eq. (10), namely, (1) evapotranspiration ET from the flux tower by the 50 ha plot (⋍ a 9 + a : in Eq. For ELM-FATES calibration we calculated an objective function for each of the 5000member ensemble by equally weighting standardized RMSE between observations and simulations across all fluxes and states mentioned above. We describe the objective function i as: based on monthly sums for both ET and `run and daily averages for θs. ^l iq,% is the average of daily ^l iq,% , where . ∈ R; R = (0.01, 0.03, 0.06, 0.12) m. Superscript * indicates that RMSE is standardized between 0 and 1 across the 5000 simulations. Model fit to soil moisture dynamics by depth was one of our key priority. Preliminary analyses showed that fitting the model to all states and fluxes simultaneously captured soil moisture dynamics for shallower depths but not at 1 m depth. To ensure model fit at 1 m depth, we thus imposed a two-step procedure, wherein we first chose the 300 parameter ensemble members that minimized RMSE * (θ 1 m ,^l iq,1 m ) (< 0.2), and among those 300 ensemble members we chose 100 with the lowest values for the objective function i.
Dataset S1 Microclimatic and flux tower data.
The tower used for these measurements is 41 m above ground, on a plateau on BCI at the north west corner of the 50ha plot (Pau et al., 2018). The eddy covariance system includes a sonic anemometer (CSAT3, Campbell Scientific, Logan, UT) and an open-path infrared CO2/H2O gas analyzer (LI7500, LiCOR. Lincoln, NE). Hi-frequency (10Hz) measurements were acquired by a datalogger (CR1000, Campbell Scientific) and stored on a local PC. Data were processed with a custom program using a standard routine described in (Detto et al., 2010). QA/QC criteria for removing erroneous values are listed in Table S2. GPP was derived from daytime values of NEE by adding the corresponding mean daily ecosystem respiration obtained as the intercept of the light response curve (Lasslop et al., 2010). The light curve was fitted on a 15-days moving window using a rectangular hyperbolic function (runs with friction velocity less than 0.4 m s -1 were excluded). In order to compute daily time integrated budgets, gaps were filled using Artificial Neural Network (Papale & Valentini, 2003) with hydro-meteorological inputs as predictors (soil moisture, solar radiation, temperature, VPD and air pressure, all measured on the same tower). To train the network, the dataset was randomly divided in a training set (70%), a validation set (15%) and a test set (15%). A two-layer feed-forward network with 10 sigmoid hidden neurons and linear output neurons, was trained using the Levenberg-Marquardt algorithm until the mean square error (MSE) of the validation set stopped improving (Hagan & Menhaj, 1994). Performance, in terms of MSE, was evaluated using the test set at the end of the training. This procedure was repeated 100 times to produce 100 estimates of GPP. Training multiple times generates different results due to different initial conditions and random sampling of the three sets. Ensemble was obtained as weighted average from the 100 ANN predictions using the MSE of the test set as weights according to: @AA ? /\:a )? ∑ 1 ? /\:a )? The ANN was implemented using the Neural Network Toolbox in Matlab 2014a.

Dataset S2
Leaf hydraulic conductivity and vulnerability to cavitation.
& leaf was measured under steady state and high irradiance following the Evaporative Flux Method (FEM) (Sack & Scoffoni, 2012) using a flowmeter (J. Zailaa et al., unpublished). Briefly, the leaf was excised under filtered deionised water and attached to a water source (graduated cylinder) using clear plastic tubing. The system includes tubing of a known resistance, separated upstream and downstream by two sensors measuring the transpirationdriven water flow pressure in the EFM system (Sack et al., 2011). & leaf is calculated as the ratio of the transpiration rate (Eleaf; mmol s -1 ) to the difference in flow pressure (MPa) upstream and downstream from the resistance tubing, and further standardized to 25°C and by leaf area. At stability (CV < 5%), the leaf is left to equilibrate in plastic bags with moist paper towels and measured for Ψleaf using a Scholander pressure chamber. To construct vulnerability curves for each species, four different functions were fitted to species-specific & leaf vs. Ψleaf data using the anneal optimization function in the R package likelihood: exponential without an intercept (Sack & Scoffoni, 2012). The differences among AICc values for these functions were not large within species. We chose the exponential function without the intercept to fit observed data across species, and as a general function to scale-up to the community, because it had the lowest AICc for most species, captured the form of the observed data across all species well, and had the lowest number of parameters, that is, two.
We used WSG, the wood specific gravity after drying at 100∘C, and LMA, the leaf mass per unit area measured for the leaf lamina excluding the petiole and for compound leaves the petiolules for leaves receiving direct sunlight as described in Wright et al. (2010).
Among the 21 species with observed data for (& leaf vs. Ψleaf) used to build the community-level models of leaf vulnerability curves, two lacked data for WSG and five in LMA. We filled three gaps in LMA by substituting species-level means for LMAdisc, the latter defined as the mean leaf mass per unit area measured for a 1.483 cm 2 leaf disc taken to avoid veins (g m -2 ) for leaves receiving direct sunlight, in a linear relationship between LMA and LMAdisc (LMA = −2.25 + 1.05LMA -)!I ; Adj. A 2 = 0.87; X < 0.001, 4 = 202). We filled the two remaining gaps in LMA by substituting LMA measured during the & leaf campaign, LMAcampaign, in a linear relationship between LMA and (LMA = 25.87 + 0.97LMA I7>J7):+ ; Adj. A 2 = 0.59; X < 0.001, 4 = 16). We chose not to fill the two data gaps in WSG as available data were poor predictors of WSG.

Dataset S4 Leaf deciduousness categories.
At Barro Colorado Island SJW and Osvaldo Calderón have drawn on 62+ years of experience in Panama to score each species as one among four leaf habits (Meakem et al., 2018): evergreen species (4 = 17, for species with ERD) keep most but not all of their canopy in the dry season; obligate deciduous species (4 = 25) shed all their leaves in the beginning of the dry season; facultative deciduous (4 = 22) species drop leaves steadily as the dry season develops, and both facultative and obligate deciduous species leaf out in the beginning of the wet season; brevi-deciduous species (4 = 23) shed all their leaves in the dry season, but leaf out in a few days to six weeks under ongoing dry season.
Dataset S5 Stream discharge from the Conrad catchment.
The Conrad weir is located on the gently sloping western side of Barro Colorado Island. The weir consists of a 90 degree 'V' notch with a rectangular upper section. The weir was constructed in several stages between 1993 and 1996. Stream stage is recorded at 5-minute intervals. The weir drains a catchment of approximately 40 ha of which at least 90% consists of the central plateau area of the island. Peak flow usually follows peak rain by approximately 2-3 hours. Large storms are sometimes preceded by small, initial peaks generated by the rain falling on the steeper slopes close to the weir. Streamflow data are analyzed by software custom-built by author SRP. The program is used to correct for blockages of the weir 'V' by fallen vegetation (a common occurrence), sensor drift, small gaps and other problems. Discharge is then calculated using an empirically derived rating curve calculated for each 1mm stream stage.

Dataset S6 Volumetric water content.
We collected volumetric water content (VWC) data (2012-2018) from three locations near the eddy flux tower with time domain reflectometry (TDR, CS616, Campbell Scientific) probes inserted vertically ranging over the first 15 cm of soil. The relationship between the output from the TDR probe and VWC was calibrated based on measurements of gravimetric water content (GWC) on several occasions from representative soils near the vertical probes (0-0.15 m) (Fig. S4). GWC was converted to VWC based on the observed mean bulk dry density of 0.8 g cm -3 near the TDR probe (4 = 5, cylinder volume = 1250 cm 3 ). Beginning in 2016 we additionally installed three horizontal probes at depths of 0.15, 0.4 and 1 m and we applied the same calibration equation to these data (2016)(2017)(2018) as that for the vertical probe. These data are archived at https://ngt-data.lbl.gov/ (DOI pending).
We also leveraged existing GWC measurements by Kupers et al. (2019) at the depths of 0.15, 0.4 and 1 m in summer 2015 (February, March and April) and 0.15, 0.4 in summer 2016 (March), with 1299 samples in total that covered all soil types and habitats in the 50ha plot. The sampling periods were six, five, ten and eight days long, respectively. The 2016 dry season was associated with the 2015-2016 El Niño (see Kupers et al. (2019) for further details). We converted these GWC measurements to VWC assuming a bulk density of 0.8 g cm -3 . c 0.15 m , c 0.4 m and c 1 m refer to the plot-wide VWC averages by depth over each sampling period. Middle date of a sampling period is used for comparison with model simulations.

Dataset S7 Stem maximum hydraulic conductivity and vulnerability to cavitation.
Values of stem maximum hydraulic conductivity and vulnerability to cavitation were obtained from the dataset of Wolfe et al. (2021). The entire dataset includes measurements on 26 tree species collected at three sites in Panama. The raw data, R script for the analyses, and data summaries are archived in the NGEE-Tropics archive (Wolfe et al., 2021). The data were previously reported by Dickman et al. (2019) and Wu et al. (2020). Here, the seven species that overlapped with our assessment of ERD were included in analyses (Table S3).
Methods for stem maximum hydraulic conductivity and vulnerability to cavitation closely followed Christman et al. (2012) and Wolfe et al. (2016). Briefly, canopy branches were collected from 4-10 trees of each species. Segments measured for hydraulic conductivity were removed from the branches and placed in an apparatus similar to that described by Sperry et al. (1988) except that flow rate was measured with graduated pipettes instead of a balance. Each segment was perfused with a 10 mM KCl solution filtered to 0.2 mM at three pressure heads: two that ranged 0.5-2 kPa and at zero pressure. Following Torres-Ruiz et al. (2012), hydraulic conductivity was calculated as the slope of the flow rate versus pressure gradient (i.e., pressure head divided by segment length) among the three pressure heads. Stem area specific hydraulic conductivity (& stem ) was calculated by dividing hydraulic conductivity by the cross-sectional area of the stem. For each & stem segment, three measurements of stem water potential were made on stem sections that were located adjacent to the & stem segment using psychrometers as described by Wolfe (2017).
To assess vulnerability to cavitation, the branches were allowed to dry in the lab for 0 -180 hours before they were measured for stem water potential and & stem . For each species, & stem was plotted against stem water potential and a Weibull curve was fit through the 90th percentile of the points, following Wolfe et al. (2016). The intercept of the curve was interpreted as maximum & stem and the water potential where the Weibull curve predicted 88% loss of maximum & stem was interpreted as Ψ88,stem.

Dataset S8 Leaf turgor loss point.
For each tree species measured for & stem and Ψ88,stem, two to six leaves were measured for Ψtlp with pressure-volume analysis following Koide et al. (1989). Leaf water potential was measured on leaf discs with psychrometers (J.R.D. Merrill Specialty Equipment, Logan, Utah, USA). This method has been shown to give results similar to those of the pressure chamber method (Nardini et al., 2008). To construct the pressure-volume curves, leaf discs were repeatedly measured for water potential and mass. Between measurements the discs were bench dried for 5-30 minutes. Dry mass was determined after oven drying the discs at 70°C.

Dataset S9
Above-ground hydraulic safety margins.
Above-ground hydraulic safety margin was calculated as the difference between minimum leaf water potential Ψmin and Ψ88,stem. Values for Ψmin were obtained from the dataset of Wolfe et al. (2019), which was collected as part of the 2016 ENSO campaign. Pre-dawn and diurnal leaf water potentials were measured on a monthly basis from Feb to May 2016 at San Lorenzo and Parque Natural Metropolitano, Panama, for which leaves were sampled using a canopy crane. Measurements in BCI were available only for March 2016, for which leaves were sampled with a pole pruner, as there is no canopy crane on BCI. In total 25 species were sampled, with no species re-sampled on another site. The raw data and data summaries are archived in the NGEE-Tropics archive . Here, the six species that overlapped with our assessment of ERD and Ψ88,stem were included in analyses. We defined species Ψmin as the minimum of all the diurnal measurements for a species after taking an average value of measurements on two leaf samples for a given diurnal measurement (Table S3). Notes S1 ERD model structure selection. The best ERD model structure explained a large fraction of the variance in δ 2 Hxylem. This model included an effect of VPD and not LAI (Eq. (S3)). Two of the alternative models (Eqs. (S4) & (S6)) with LAI had a strong relationship with δ 2 Hxylem as well, but one that was based on a smaller set of species (4 = 4 or 5), because fewer species passed the criteria of sufficient species-level goodness-of-fit for these models (<17 species vs. 29 species for the best model, Fig. S11). For the rest of the analyses, we therefore used ERD estimates from the less parsimonious model (Eq. (S3)).

Notes S2
Additional results for exposure to water-stress.
Overall, species with ERD deeper than 2.9 m did not cross species Ψcrit for the entire 25yr period for which Ψ soil,% was modeled . Soil depths deeper than 2.9 m remained above -0.06 MPa (-0.02, -0.14; median, 95% CI) through all the dry seasons and drought periods, which is above Ψcrit of the most sensitive tree species (-0.17 MPa).