Assessing Role of Drought Indices in Anticipating Pine Decline in the Sierra Nevada, CA

: Tree mortality in Sierra Nevada’s 2012–2015 drought was unexpectedly excessive: ~152 million trees died. The relative performance of ﬁve drought indices (DIs: SPEI, AI, PDSI, scPDSI, and PHDI) was evaluated in the complex, upland terrain which supports the forest and supplies 60% of Californian water use. We tested the relative performance of DIs parameterized with on-site and modeled (PRISM) meteorology using streamﬂow (linear correlation), and modeled forest stand NDVI and tree basal area increment (BAI) with current and lagged year DI. For BAI, additional co-variates that could modify tree response to the environment were included (crown vigor, point-in-time rate of bole growth, and tree to tree competition). On-site and modeled parameterizations of DIs were strongly correlated (0.9), but modeled parameterizations overestimated water availability. Current year DIs were well correlated (0.7–0.9) with streamﬂow, with physics-based DIs performing better than pedologically-based DIs. DIs were poorly correlated (0.2–0.3) to forest stand NDVI in these variable-density, pine-dominated forests. Current and prior year DIs were signiﬁcant covariates in the model for BAI but accounted for little of the variation in the model. In this ecosystem where trees shift seasonally between near-surface to regolithic water, DIs were poorly suited for anticipating the observed tree decline.


Introduction
Over the last two decades, droughts have been more extreme, extended, and hotter in the Pacific and Southwestern U.S. [1][2][3][4][5], challenging our current understanding of drought and its effects on forests and humans alike. Massive tree die-offs in the Sierra Nevada of California [6,7] and surprise at its magnitude [8,9] have drawn world-wide attention. Linking the level, duration, and scale of drought to impacts on hydrological and biological resources is a critical step in anticipating future forest resource response and loss [10].
Water availability is a demonstrably vulnerable source of water for these forested ecosystems, which rely on, buffer, and supply water for human use. The Sierra Nevada provides~60% of Californian water resources utilized by humans [11]. Water resources are evaluated via drought indices (DIs), and these are used to assign the level of drought severity and help communicate and allocate near-term water resources for direct human consumption, agriculture, and industrial uses. Accurate estimates of water resources and their effects on biological responses are also needed (ecological drought) [12]. Modeled meteorological data (PRISM, or Parameter-elevation Relationships on Independent Slopes Model) [13,14] are broadly used for a range of applications, and in particular to calculate DIs. The performance of PRISM in the foothills of the Appalachian Mountains was tested and the model performed well [15]. On the eastern slope of the Sierra Nevada, PRISM sufficient) water sources trapped in bedrock interstices [23,24] or pr streams, or rivers. Xeric microsite trees presumably experienced ground seeps, pockets of water trapped in bedrock interstices, or sto Overall, the Jeffrey pine investigated here averaged 32 + 7 cm in bo (1.37 m), 104 + 36 yr old at DBH, and 17 + 5 m in height. At each site the initial number of trees selected was 32 in 1998; about 10% of the time to successful bark beetle attack and/or drought. The level of ph stress for each tree at MFm and MFx was well-characterized in 1999 unpubl. data), and in 2013 from thermal imagery [45]. All trees used i alive in September 2017 at MF and HM, and in September 2019 at SC.

Work Flow
The procedure of this work was as follows. Monthly on-site meteo compiled from different sources and missing data were imputed by r Monthly on-site and modeled (PRISM) meteorological data were com data were examined for discrepancies at each site. Five commonly us eterized with on-site and PRISM meteorological data. For PDSI-base known pedological values, both suggested default values were teste mm). For each site and DI, sources of meteorological data were com correlation coefficients. We tested the relative performance of DIs para site and modeled meteorology using streamflow (linear correlation), stand NDVI and tree basal area increment (BAI) with current-and predicting BAI, current and lagged DIs were placed within the conte competition metrics using a random effects model. The effectivenes was evaluated using p-values of the t-test and the likelihood ratio test

Work Flow
The procedure of this work was as follows. Monthly on-site meteorological data were compiled from different sources and missing data were imputed by rule-based methods. Monthly on-site and modeled (PRISM) meteorological data were compared and modeled data were examined for discrepancies at each site. Five commonly used DIs were parameterized with on-site and PRISM meteorological data. For PDSI-based DIs requiring unknown pedological values, both suggested default values were tested (152 mm and 100 mm). For each site and DI, sources of meteorological data were compared using linear correlation coefficients. We tested the relative performance of DIs parameterized with on-site and modeled meteorology using streamflow (linear correlation), and modeled forest stand NDVI and tree basal area increment (BAI) with current-and lagged-year DI. For predicting BAI, current and lagged DIs were placed within the context of tree vigor and competition metrics using a random effects model. The role of DI's in predicting BAI was evaluated using the fitted models.

Collating Meteological Data
Meteorological data were extracted from existing records of daily precipitation, and minimum and maximum air temperature extracted from October 1933 to September 2020 for HM and MF and from October 1939 to September 2020 for SC. The data were obtained from the National Oceanic Atmospheric Administration National Center for Environmental Information Database (NOAA NCEI, [46]), Sequoia National Park archives (obtained from E. Meyer, Research Division, Sequoia-Kings Canyon National Parks), the Western Regional Climate Center (WRCC, [47]), and the National Atmospheric Deposition Program (NADP, [48]). The meteorological data (monthly precipitation sum, average mean air temperature), sources, rule-based methods for accommodating missing data, and metadata have been published [49]. Calculations of potential evapotranspiration (PET, [50]) and parameterizations of the five drought indices (see next section) at the three sites are also contained in posted data. Data [49]. The 800 m pixel data was not available for these sites.

Parameterization of Drought Indices
All DIs were parameterized with on-site and PRISM meteorological data and reported on an annual ("12") and hydrologic year ("h," 10/1 prior yr to 9/30 reported yr) basis, hence "SPEI12h." All five of the DIs evaluated required calculation of PET. We used the Thornthwaite [50] equation here because it requires only temperature and latitude. The Penman-Monteith calculation is preferred in more moist environments [51,52]. For each of the three sites, the two parameterizations of DIs were tested through the strength of the correlation coefficient. The five DIs were similarly tested among themselves using on-site and then PRISM parameterization.
SPEI was parameterized with summed monthly precipitation, averaged monthly mean air temperature, and location on Earth (latitude as a constraint on radiant energy driving evaporation) to calculate potential evapotranspiration (PET). PET was subtracted from precipitation to quantify deficits and expressed as a distribution-based computation for the convenience of all positive numbers when modeling (log)BAI. Aridity Index (AI) used the same data to calculate the ratio between summed precipitation and mean air temperature [20].
PDSI was developed to quantify drought intensities for crop production and thus considered not only gain and loss of water, but also moisture stored in the soil [21], which SPEI and AI do not. PHDI is a modified version of PDSI and accounts for both development of and recovery from drought [53], with the primary modification of assessing drought termination. PHDI captures the impact of long-term droughts on hydrological balance. While the end of drought is marked by the ratio of PET to precipitation exceeding 0% for PDSI, this ratio is set to 100% in PHDI. scPDSI is another variant of PDSI introduced to improve consistency across regions [22]. For PDSI and derivatives (scPDSI, PHDI; [21,22]), we calculated a monthly drought measurement, the Palmer Z index, which was then summed annually. Because soil available water content (AWC) was not available for our sites, we used the two recommended default values for each PDSI and its derivations, 100 mm and 152 mm. The R package 'scPDSI' was used for these calculations [54]. The relative performance of current-or prior-year DIs was assessed with the strength of the correlation coefficient, parameterized with on-site or PRISM meteorological data, with a 7 yr hydrologic record (2014 through 2020) of the second-order portion of the Marble Fork Kaweah River. This was the only gauged river immediately associated with our sites, located 0.5 km upstream from the midpoint of the MF site. Streamflow data was recorded every 5 min, extracted from the United States Geological Survey (USGS), National Water Information System [55], and summed monthly.

NDVI
Normalized Difference Vegetation Index (NDVI) is a measure of vegetation response to environmental and meteorological conditions, (NIR−RED)/(NIR + RED), where NIR is near-infrared wavelengths and RED is visible red wavelengths. NDVI was extracted from 18 February 2000 to 29 September 2020 from the NASA Moderate Resolution Imaging Spectroradiometer (MODIS aboard Terra and Aqua satellites). The MOD13Q1 product was corrected for atmospheric and illumination effects and was produced using 16 d composite periods. NDVI values during snowy and cloudy conditions were removed, and the response was smoothed using the cloud computing platform Google Earth Engine [56]. For each site, the maximum number of 250 m × 250 m MODIS pixels (HM, n = 3; MFm, 5; MFx, 3; SC, 5) contained within the polygon of the site were extracted and averaged for the two sample dates in August reflecting the maximum effect of drought on vegetation or summed over the growing season at each site. The growing season length was defined by start-of-season mean air temperature >7 • C to end-of-season mean air temperature <0 • C, resulting in site-specific and sub-site (MFm and MFx) annual values (temperature choices influenced by [57][58][59]. Similar to the approach taken for streamflow, the relative performance of current or prior year DIs was assessed with the fitted linear model R 2 , parameterized with on-site and PRISM meteorological data, with mean August, or summed NDVI over the growing season.

BAI
Based on a previous analysis of a closely related species (Pinus ponderosa Dougl. ex Laws) also evaluated in dry pine-dominated woodlands to forests, an approach to modeling annual BAI was developed [60] testing current yr and up to 3 yrs prior, and 2 through 4 yr running averages of the DI chosen, as well as a number of crown vigor attributes (see below), insect and disease incidence, and tree-to-tree competition metrics. In the current analysis, we retested current, lagged, and running averages of DIs and biological covariates for significant explanatories that explained the greatest proportion of the variance in predicting BAI of Jeffrey pine. HM trees were last cored in 2006, and a number of trees had heart rot or resinous centers, so a shorter evaluation period was available . Trees at HM were analyzed separately from the other sites . We modeled BAI of 15 trees at HM using a 25 yr BAI record. MF (MFm, MFx) and SC trees were last cored in 2017 and the BAI chronology was evaluated from 1955 to 2016. We modeled BAI of MFm (n = 17), MFx (n = 16), and SC (n = 20) using a 61 yr BAI record. Because SPEI and AI are highly correlated, and the PDSI variants are highly correlated, we chose only SPEI12h and PDSI12h 152 for this analysis. Significance of SPEI or PDSI metrics and the proportion of variance explained in the full model were used to evaluate the two DIs parameterized with on-site and PRISM parameterizations.
Ring width in a cross section of a tree bole is a precise measurement, from which BAI was calculated [61]. Cores were taken within 10 cm of diameter at breast height (DBH, 1.37 m) above the average soil surface to tree center in early to late October, avoiding irregular bark or bole imperfections. Cores were mounted on blocks and sanded to 400 grit. Annual ring widths (early + late wood) were measured to 0.001 mm (Velmex, Inc., Bloomfield, NY, USA) using program J2x, Voor Tech Consulting (2008) and checked for missing rings [62] based on cross-correlation with a regional chronology (Merschel, unpubl. data). Dating accuracy was evaluated using the software program COFECHA (version 6.06P; [63,64]. Potential dating errors identified were visually checked, re-dated, re-measured, or re-collected as necessary. For cores that did not intersect the pith, we estimated the number of rings to the pith geometrically [65] for not more than 4 yrs. BAI for tree i was calculated from radial growth expressed on bole area in the year t minus bole area in year t − 1: where r t is bole radius at time t. Tree vigor can affect the capacity of the tree to respond to the environment. 'RANK' is a qualitative metric, comprised of several qualitative and quantitative measures (visual assessment of needle color, length, and retention (number of needle ages; retention within a needle age class), branchlet length and diameter), and tree crown structure described and developed for ponderosa pine [38] and applied to Jeffrey pine previously in MFm and MFx [66]. DBH/tree age, here referred to as 'bole vigor:' [barkless tree diameter in 1980]/[tree age in 1980 at DBH], a point-in-time measure, roughly a mid-point of this analysis. Competitive zone density (CZD) is a measure of the tree-to-tree competition of each target tree. In each of the four aspects (NE, SE, SW, and NW) to 20 m around the target tree, the DBH of the nearest conspecific tree (>10 cm DBH) divided by the distance from target to neighbor, center to center was measured, and the four measures (DBH/distance) were averaged [40]. This metric was superior to all others tested (variable radius plots; resilience capacity; trees per ha; basal area per ha; and number of trees per cluster [60]) in assessing individual tree response to conspecific competition in its immediate neighborhood.

On-Site vs. PRISM Meteorology
The correlation coefficient between on-site and PRISM precipitation and temperature was 0.98-0.99, respectively, for all three sites. The comparisons are given for MF in Figure 2, and those for HM and SC given in Appendices A and B. Although monthly PRISM precipitation and temperature data appeared to be well matched to on-site data, seasonal discrepancies were apparent. PRISM temperature data was higher than on-site measurements early in the hydrologic year (November to January) and lower from April to June. An extreme cold event was not predicted by PRISM. There was greater than a 6 • C discrepancy between modeled and on-site temperature during the third week of December 1992 and the first and second weeks of January 1993, which decreased these monthly averages of mean daily air temperature [67].
The discrepancy between on-site and PRISM precipitation was greatest during the period of greatest precipitation, December through April (Figure 2c), when 80% of the annual precipitation is received. Extreme precipitation events (1/1969, 1/1980, 2/1986, and 1/2017) were well outside of the 5% error, and PRISM averaged total annual precipitation for HM, MF, and SC was greater (113 cm, 119 cm, and 115 cm, respectively) than that of on-site data (110 cm, 114 cm, and 106 cm, respectively).

On-Site VS. PRISM Parameterization of Drought Indices
All DIs generally tracked known high and low precipitation years (MF: Figure 3, HM, SC: Appendices C and D). Linear correlation coefficients between on-site vs. PRISM meteorological parameterizations of DIs were very high for those evaluated, ranging from 0.87 to 0.99 (Table 1). Although PET and precipitation are components of the DIs investigated, SPEI12h and AI12h were dependent on precipitation to a greater extent than PDSI and its derivatives. SPEI12h had the highest correlation coefficients, and scPDSI had the lowest with 100 mm AWC was slightly lower than that at 152 mm. PDSI and derivatives require additional calculations above those needed for SPEI and AI. We suggest that scPDSI may have the lowest correlation coefficient because it is the most modified from the underlying meteorology (precipitation and temperature) and may be less suited to the hydrology of these sites. Despite strong correlations between the two parameterizations of DIs, PRISM parameterizations reflected the November through June temporal discordance in precipitation ( Figure 3). The findings for precipitation at MF also held for HM and SC (Appendices C and D).

On-Site vs. PRISM Parameterization of Drought Indices
All DIs generally tracked known high and low precipitation years (MF: Figure 3, HM, SC: Appendices C and D). Linear correlation coefficients between on-site vs. PRISM meteorological parameterizations of DIs were very high for those evaluated, ranging from 0.87 to 0.99 (Table 1). Although PET and precipitation are components of the DIs investigated, SPEI12h and AI12h were dependent on precipitation to a greater extent than PDSI and its derivatives. SPEI12h had the highest correlation coefficients, and scPDSI had the lowest with 100 mm AWC was slightly lower than that at 152 mm. PDSI and derivatives require additional calculations above those needed for SPEI and AI. We suggest that scPDSI may have the lowest correlation coefficient because it is the most modified from the underlying meteorology (precipitation and temperature) and may be less suited to the hydrology of these sites. Despite strong correlations between the two parameterizations of DIs, PRISM parameterizations reflected the November through June temporal discordance in precipitation ( Figure 3). The findings for precipitation at MF also held for HM and SC (Appendices C and D).

Correlation among Drought Indices
Linear correlation coefficients between paired DIs within a parameterization (onsite or PRISM) are summarized in Table 2. The correlation coefficients were high and ranged from 0.7 to nearly 1.0. SPEI12h, and AI12h had high correlation coefficients in both on-site and PRISM parameterizations (0.98). Each of the Palmer indices had small differences in correlation coefficients when comparing the two AWC default values. The correlation coefficients were higher between PDSI derivatives at the same AWC level (PDSI152, scPDSI152) and between DIs parameterized with PRISM meteorology, except for PHDI12h where the highest coefficient was obtained with the 100 mm AWC default value (0.989).

Drought Index Performance Relative to Short-Term Streamflow
The relative performance of the DIs was tested using the strength of the correlation coefficients between current year SPEI at MF and current year streamflow of the Marble Fork Kaweah River over a 7 yr period. SPEI12h and AI2h had the highest correlation coefficients (R > 0.98), using on-site parameterization (Table 3). PDSI and derivatives performed better with PRISM than on-site meteorology perhaps because of a similar spatial scale. Among the PDSI derivatives, PDSI12h performed the best, with slightly improved correlation coefficients using the default AWC of 152 mm. Temporally lagged streamflow (1 yr) was not significantly correlated with precipitation (R = 0.30).

Drought Index Performance Predicting Site NDVI
Using a model to predict NDVI, the relative performance of the DIs was tested, parameterized with either on-site or PRISM meteorology. A given DI is denoted by D t with the specified parameterization in year t. We compared the mean August NDVI at HM, SC, and the sub-site MFx against the reference site, MFm. Letting Y s,t be the average August NDVI at site s in year t, we fit the linear model: Here, { s,t } centered at zero as independent and identically distributed (i.i.d), Gaussian error terms. I s=S is the 0-1 indicator that the observation was at site S. Each DI was standardized (each variable was centered with unit standard deviation) which permitted comparison of coefficients across indices. We also standardized the year variable. Figure 4 shows the data together with the predicted average August NDVI as a function of SPEI, holding all other variables in Equation (2) at their mean values 0. The estimates and t-tests of the coefficients in Equation (2), fit with DIs parameterized with both on-site and PRISM meteorology, are shown in Table 4. The significance of the DI coefficients varied among the fitted models, and the 1 yr lagged DI was significant for SPEI and AI but not for PDSI 152 or PDSI 100 . Across the DIs and AWC levels tested, and the source of the meteorological data, the model could not differentiate MFm-or MFx-dominated subsites, the latter with quantified differences in tree physiological drought stress between them [44,45]. Only SPEI and AI with on-site parameterizations differentiated SC from MFm. HM was not statistically distinguished from MFm in any of the models. Nonetheless, the overall fit of the model among the indices is similar. The R 2 of the models ranged from 0.22 to 0.32. All models designed to predict NDVI with drought indices were similar to one another and of low predictive effectiveness. There was no predictive capacity of the DIs and site-specific, average August, or summed growing season NDVI.  To accommodate one tree with no measurable ring width in one year, 0.01 was added to enable a non-zero log transformation of BAI (log(BAI)) shown in Figure 5. Most HM trees were linearly increasing in BAI early in the chronology shown. Four trees had a different growth pattern from the others, with a peak between 1994 and 1997, suggesting reliance on a water source different from the others at the site. Of the four trees, two were growing out of a rock outcrop and two were in a swale, in deep soil. Adjacent trees in these two settings did not have the same growth pattern. Current year and temporally lagged DI were included in the model and were denoted by D t , Through a preliminary and exploratory leave-one-out multiple regression, we identified additional factors that modified tree response to the DIs: RANK (with variable R, and the coefficient β R ), DBH/AGE (H, β H ), and CZD (C, β C ). Temporal change in BAI that was not attributable to the above variables was incorporated in a linear term, YEAR. This variable was scaled to avoid computational issues. A tree random effect was also included to model tree to tree variability. Our model for tree i in year t was: where B 0,i and B 1,i are i.i.d. centered Gaussian variables. Thus, there was a per-tree random effect accounting for between-tree variability of size and its growth trend. We fit a linear mixed model to the data at HM using SPEI then PDSI, in onsite and PRISM parameterizations. There was no difference in overall fit between the two indices with the same parameterization, and within a DI, the two parameterizations; the marginal R 2 's were the same: 0.24.
Climate 2022, 10, x FOR PEER REVIEW 3.6.2. Predicting BAI at MF and SC: On average, MFm trees had slightly higher BAI (log(BAI)) than MFx; MFx BAI were similar ( Figure 6). For many trees at SC, BAI appeared to be indepen known years of limited precipitation. The following years had < 80% of long-term a precipitation: 1959-1961, 1964, 1966, 1968, 1970-1972, 1976-1977, 1981, 1987-199 1999, 2004, 2007, and 2012-2015. At SC, two trees had markedly different BAI grow terns, one within the annual flood plain of Stoney Creek, and the other on a conv slope. When these trees were removed from the model, there was little change in cance or in our interpretation of the model. The response of log(BAI + 0.01) g against SPEI12h or PDSI12h was similar ( Figure 6).
In the model below, we concurrently predicted BAI of MF and SC, including for site and MF subsites. The terms IA are dummy variables with the value of 1 or 0 corresponding to w or not condition A was met. The smoothing terms sMFx(t), sMFm(t), and sSC(t) are fu of (t), and B0,i and B1,i are tree-dependent, centered, random Gaussian terms. Each standardized to have mean 0 and unit standard deviation, enabling comparison o cients between indices.
Fitted models with site-specific parameterization for each DI are presented i 6. Current and 1 yr lagged DI were significant predictors for the on-site parameteri although the effect size was very low. Plots of the predicted log(BAI) as a function and PDSI, fixing all other variables to their mean values are presented in Figure the dynamic range of SPEI and PDSI, predicted log(BAI) did not increase substa For example, considering the PDSI at SC (Figure 7, lower right), the expected lo To assess the role of DI in the model, we compared the full model (3) to one omitting the 4 terms consisting of the DI and its time lags: The likelihood ratio test comparing Equation (4) to Equation (3) for the two DIs and two parameterizations suggested that there was insufficient evidence to reject the model omitting the DI variables (Table 5). Moreover, the reduced model omitting all DI variables had smaller AIC and BIC in all cases, suggesting the DI variables did not enhance predictions in Equation (4) using only tree attributes. Removing the four trees with different patterns of growth did not change the statistics appreciably, nor our interpretation. Perhaps because different trees may have been relying on different water sources, including the pixel-based meteorological parameterization (PRISM) created a more effective model. Finally, we evaluated model performance using just DI and temporal lags to predict BAI. We fit the following simple mixed linear model: The marginal R 2 for this model (Equation (5)) suggested little difference between indices (SPEI and PDSI on-site, marginal R 2 = 0.03, 0.05, respectively; and SPEI and PDSI PRISM, marginal R 2 = 0.04 for both) and demonstrated that the indices alone were inadequate for predicting BAI at HM.

Predicting BAI at MF and SC:
On average, MFm trees had slightly higher BAI (log(BAI)) than MFx; MFx and SC BAI were similar ( Figure 6). For many trees at SC, BAI appeared to be independent of known years of limited precipitation. The following years had <80% of long-term average precipitation: 1959-1961, 1964, 1966, 1968, 1970-1972, 1976-1977, 1981, 1987-1992, 1994, 1999, 2004, 2007   In the model below, we concurrently predicted BAI of MF and SC, including a factor for site and MF subsites. Let j(i) denote the site of tree i. The model for tree i at time t was: The terms I A are dummy variables with the value of 1 or 0 corresponding to whether or not condition A was met. The smoothing terms s MFx (t), s MFm (t), and s SC (t) are functions of (t), and B 0,i and B 1,i are tree-dependent, centered, random Gaussian terms. Each DI was standardized to have mean 0 and unit standard deviation, enabling comparison of coefficients between indices.
Fitted models with site-specific parameterization for each DI are presented in Table 6. Current and 1 yr lagged DI were significant predictors for the on-site parameterizations, although the effect size was very low. Plots of the predicted log(BAI) as a function of SPEI and PDSI, fixing all other variables to their mean values are presented in Figure 7. Over the dynamic range of SPEI and PDSI, predicted log(BAI) did not increase substantially. For example, considering the PDSI at SC (Figure 7, lower right), the expected log (BAI) increased from a minimum of~2.75 to a maximum of~2.9 while holding all other variables at their mean values. Tree-to-tree competition was significant for both DIs and parameterizations. Due to the variation among trees, differences among sites and subsites were not significant. All models had similar fits; the marginal R 2 for both DIs and parameterizations had the same value: 0.62. Thus, both SPEI and PDSI, parameterized with either on-site or PRISM meteorology performed similarly, and moderately well in predicting BAI at MF and SC. Leaving out 'bole vigor' from the otherwise full model reduced the predictive capacity of the model (marginal R 2 = 0.23).
To assess the effectiveness of SPEI and PDSI 152 alone in predicting BAI with onsite or PRISM parameterization without additional biological covariates, we fit an additional model: This model exhibited poor performance, with the marginal R 2 low (0.14) for b and parameterizations. For the on-site parameterization of SPEI, the estimate of the ard deviation of the random effect (B0,i) was 0.73. This was large relative to the r error, which had an estimated standard deviation of 0.50. Other indices and para zations had similar estimated standard deviations.  (7)) are fixed at their mean val crown vigor, bole growth vigor, and tree competition). This model exhibited poor performance, with the marginal R 2 low (0.14) for both DIs and parameterizations. For the on-site parameterization of SPEI, the estimate of the standard deviation of the random effect (B 0,i ) was 0.73. This was large relative to the residual error, which had an estimated standard deviation of 0.50. Other indices and parameterizations had similar estimated standard deviations.

Discussion
Two types of DIs were evaluated, one based on meteorological and physics-based inputs (SPEI and AI), and the other on meteorological-, physics-, and pedogenic-based inputs (PDSI and derivatives scPDSI and PHDI) in mid-elevation, complex terrain. Relatively few comparisons have been made between on-site and modeled meteorology in these landscapes. While on-site and PRISM monthly mean temperature and summed precipitation were highly correlated, periodic discrepancies in temperature were identified at our sites on the western slope of the Sierra Nevada. PRISM mean temperature was higher than that observed from November to January, and lower for April, May, and June. Higher temperatures in the winter increase ET but would have a low impact on vegetation. Lower temperatures in late spring would decrease calculated DIs and underestimate the level of drought severity when it occurs. In another study across a mid-to high elevation basin on the eastern slope of the Sierra Nevada [16], PRISM maximum temperatures were lower relative to on-site values in early spring (March, April), and mid-summer (July, August). Of interest is that there was an overlap but not an identical match of underestimated monthly temperature between the western (this study) and eastern slopes [16] of the Sierra Nevada, and that the impact of the temperature discrepency would be greater in the latter study as water deficits increase through the summer. During the period of highest precipitation, from November to March, the estimated error variance of PRISM precipitation was high in our study. Extremes, both higher and lower than received, in on-site precipitation exceeded 5% of that estimated by PRISM. Overall, the annual precipitation estimated by PRISM was 3% to 8% greater than received across the three sites. The temperature and precipitation discrepancies of PRISM relative to on-site meteorology both contribute to an underestimation of DI severity during the growing season.
An extreme cold event was not predicted by PRISM (>5 • C discrepancy between modeled and on-site temperature) in December 1992 and January 1993 [67]. Extremely low temperatures on the western slope of the Sierra Nevada are rare but can occur with strong equatorial warming and broadening and deepening of the now permanent year-round lobe of the "winter" arctic jet stream. In December 1992, 'frigid' arctic air was reported as far south as southern California, as were exceptionally low minimum temperatures for the first and second weeks of January 1993 [67]. One of the effects of variability in the winter jet stream is that aberrant storm tracks are driven further south along the Pacific coastline [68]. The arctic air mass could have come in over the top of the lower elevation marine air, moving the cold air further up the mountain, as the adjacent, lower elevation Central Valley just east of Visalia had typical daily minimum air temperatures during those events around 0 • C (PRISM daily data for 36.3815, −119.0839, 132 m; explanation suggested by R. Neilson, pers. comm.). An approach to better predicting extreme atmospheric conditions has been proposed [69].
The Marble Fork Kaweah River record was short but included years of low precipitation (<60% of the long-term average) and two of the highest precipitation years of the 90 yr record (>150% of average). Current year DIs were significantly correlated to streamflow. The physics-based DIs, SPEI, and AI, had the highest correlation coefficients (>0.98). PDSI and derived indices had moderate correlation coefficients (0.72 to 0.87), with default AWC values of 152 mm performing better than 100 mm. Prior year DI, regardless of parameterization or type, was not correlated to current year streamflow perhaps because of the short distance between the headwaters and the measurement location, reducing residence and response time [70]. In this environment with highly pulsed precipitation (90% of the annual precipitation is received between November and April) and unknown residence time, forecasts of annual streamflow can be expected to be elusive, although advances in seasonal forecasting have been made [71].
We developed a simple model using current and 1 yr lagged DIs to test their performance in predicting NDVI in the three sites. NDVI is commonly used as a proxy for vegetation productivity, a sensitive indicator of hydrologic drought [72], and an indicator of vegetative degradation from other sources of disturbance [73,74]. NDVI may be particularly applicable to sites dominated by a single species, as investigated here. MODISderived NDVI was available for 20 yrs and, except for the 1977 drought year, contained the lowest and all but the highest precipitation year (1969) in the 86 yr precipitation record constructed [49]. In our study, current and lagged year SPEI and PDSI 152 were significant explanatories (p ≤ 0.01) of site NDVI, but the linear models had low R 2 , whether NDVI was represented as a monthly average during period of greatest drawn-down of site water availability (August), or as an integration of site-specific, summed growing season NDVI. All indices and parameterizations performed similarly in tests of predicting NDVI (e.g., low R 2 ). Precipitation five months prior to the growing season (January, coinciding with the period of highest Sierra Nevada precipitation) was significantly correlated with NDVI (also MODIS-derived) [75]. The authors suggested that a long growing season and greater carbon gain in the prior year could have also influenced gross net primary production. In another study, NDVI of individual tree crowns from a fixed-wing platform was one of a number of significant metrics in differentiating low vigor from 'not low' vigor crowns in a closely related pine [76]. However, the standard deviation of crown NDVI explained more of the variation than mean crown NDVI (1st vs. 7th most important variable; Random Forest; [76]).
Jeffrey pine in mesic (MFm) and xeric microsites (MFx) could not be differentiated with MODIS-derived NDVI. Morphological crown response to drought (~30% loss of needles mid-crown in the 2007 drought; Grulke, unpubl. data at MF), physiological drought stress [44], elevated foliar oxidative stress [66], and as detected with thermal aerial imagery in 2013 ( [45]) clearly differentiated trees in the two subsites at MF. NDVI has been successfully used to validate the use of DIs in vastly different ecosystems (agricultural [77]; chaparral, [78]; forest [79]; and the African continent [80], but in the Californian mixed conifer zone with discontinuous woodlands and forest cover, neither mean August NDVI nor site-specific NDVI summed over the growing season could be used as a validation metric for DI parameterization (however, see [77]).
We then evaluated the role of DIs and the source of their parameterizations in predicting Jeffrey pine BAI. The HM site had a shorter tree ring chronology than MF and SC and was analyzed separately. At HM, none of the five DIs, with either parameterization, at either of the AWC default values, were significant explanatories of BAI. The HM site was downslope (<0.5 km) from Crescent Meadow, which may have acted as a reservoir to buffer annual precipitation inputs, or lack thereof. It is possible that the 1987-1992 drought may not have been long or extreme enough, nor was the single drought year 2002 extreme enough, to reduce the putative supply from the meadow that might have been reflected in BAI. Proximity (distance) to a meadow accounted for 8% of the total variability in canopy water content (using Random Forest, [81]) of giant sequoia (Sequoiadendron giganteum Lind. & Buch.), a species that occurs adjacent to the HM site. In the same area, the underlying geology, and proximity to a stream were significant predictors of crown water content of both giant sequoia and Jeffrey pine in 2015 and 2016 [82]. However, based on isotopic analysis [83], several deciduous desert tree species did not appear to utilize water from an adjacent stream. The trees sampled and presented here were all alive in 2019, 14 yrs after end of the tree ring record analyzed at this site. Patterns of tree BAI at HM were variable, but not more variable than at the other two sites. BAI of four HM trees had a different growth pattern suggesting reliance on a different water source than others in the stand. A more complex model was developed for MF and SC to predict BAI from current and temporally lagged DIs (SPEI and PDSI 152 ) and biological co-variates (tree crown vigor, bole growth vigor, and tree-to-tree competition). The co-variates tested for Jeffrey pine were selected from an exhaustive analysis of a closely related species, ponderosa pine [60], and were all re-tested in this analysis of Jeffrey pine. Four of the five explanatories were significant using either DI, and the model R 2 was moderate (0.62). The fifth explanatory, tree crown vigor, was not significant with either DI or parameterization source, possibly explaining why although DIs were significantly correlated to NDVI, they explained little of the variation. Temporal lags in precipitation were also found to be significant explanatories in a model forecasting the likelihood of tree mortality due to drought and insects (bark beetle, wood borers) in the Sierra Nevada province [84], as well as for temporally lagged SPEI for pine ring width in another summer-dry forest (for Scots pine, (Pinus sylvestris L.), and black pine (Pinus nigra Arn.), in Spain; [85]). Black pine was sensitive to short term drought, but at higher elevation was less so [86]. Jeffrey pine is known to be more drought tolerant than closely related ponderosa pine [87], and perhaps is less sensitive to short term drought as in stone pine (Pinus cembra L.) in the same study [86]. Despite the possibility of long vs. short term drought sensitivities, within-growing season drought significantly affected ring width (index) in black pine. Wood properties, rather than the ring width analyzed in this study, have also been shown to be more responsive than ring widths to climate signals [88], applicable to non-treeline locations.
The hydrology of mid-elevation sites on the western slope of the Sierra Nevada is complex. Jeffrey pine accesses moisture in shallow soil in the spring, but when near-surface horizons dry out, it likely relies on soil pockets in bedrock interstices, trapped moisture under cobbles and boulders, bedrock itself [42,43], surficial and underground seeps, streams, and rivers (however, see [83]). Uphill meadows may also provide a water supply that buffers the impact of droughts [81], as may be the case at our HM site. Individual trees in the same stand may access different water sources, which may vary with residence time, flow or percolation rates, and as well as a possibly unknowable amount of storage in soils and regolith (Cr) [25]. Adjacent trees did not necessarily have the same pattern of growth. The converse was also true: trees in very different microsites (swales with deep soils vs. trees emerging from exposed, cracked bedrock) exhibited similar patterns of BAI. In the central Sierra Nevada, no relationship between Jeffrey pine basal area and soil thickness (A + C, nor Cr) was found [88]. These observations suggest that trees access different sources of water and or there are temporal lags in tree response to similar sources of water.
Considering the complex hydrology and largely unknown or unknowable pedology of these sites suggests that SPEI or AI would be more appropriate indices. However, for predicting BAI in our sites, there was no detectable difference between SPEI and PDSI 152 as they exhibited the same correlation coefficient with BAI. In this study, DIs alone (e.g., no biological co-variates), whether parameterized with on-site or modeled meteorological data, were not effective in predicting BAI (R 2 = 0.14), and thus in anticipating tree mortality. Surprisingly, despite the different scales of source meteorological data (microsite, point to landscape, and pixel), there was little to no statistical difference in model performance [89].
The Sierra Nevada and Eastern Nevada Province [84] had very high tree mortality in the 2012-2015 drought, an 'unprecedented drought' [90]. Although giant sequoia grows in moist sites with deeper soils [82,91], and Jeffrey (and ponderosa) pine grows in thin soils overlying weathered bedrock [24], not all trees sampled of either species responded to the drought [81]. In giant sequoia, lower canopy water content occurred on the margins of groves, at lower elevations further from streams, or in areas of high stand density [81]. Scots pine was more susceptible to drought-induced decline and mortality on forest edges [92]. In our study, we observed high within-site variability BAI in the investigated sites. Adjacent trees in apparently similar microsites had differing BAI patterns, and trees with similar BAI growth patterns occurred in disparate microsites (streamside, dry convex slope). Similar to giant sequoia, some Jeffrey pine were not responsive to drought and others exhibited increased growth regardless of calculated drought severity (e.g., resilience to drought; [93]). In this study, SPEI and PDSI 152 with on-site or modeled meteorological parameterizations, performed similarly (R 2 = 0.62, MF/SC combined) in predicting BAI only when biological covariates were included in the model. Point-in-time bole growth rates and tree-to-tree competition were significant covariates. Point-in-time bole growth rate, here barkless DBH as a function of tree age at a mid-point in the analysis, was a significant co-variate. In a much larger analysis, larger trees had a disproportionate competitive advantage in mixed age conifer forests [39], and trees with higher bole growth rates were less likely to respond to drought defined by PDSI. The important role of stand-and tree-level competition in bole growth rates has long been established [94,95], and 'neighborhood' competition (CZD) was also a significant co-variate in our model. Trees may succumb quickly or at length to a physical or biological assault. Manion's [35] 'death spiral' suggests mortality in 5 to 8 yrs with multiple insults in long-lived trees. In south central Oregon, ponderosa pine was found to take 40 yrs to succumb after an inciting drought, based on steadily declining BAI to nil (Grulke and Mershel, unpubl. data). Jeffrey pines died in the same year on the eastern edge of the Transverse Range in 1999, the first year of a multi-year drought (1999-2002) (Grulke, pers. obs.). In that range, 25% to 75% tree mortality occurred in stands with greater than 50% tree cover [96]. In the Sierra Nevada, the 'background' tree mortality from 2010 to 2014 was 2.7 million trees annually, preceded by a decade of already lower precipitation and higher temperatures than long term records. The 2013, 2014, and 2015 mortality rates were 29, 62, and 27 million trees, respectively, dropping to 19 million trees in 2018 [6], and still elevated 3 to 6 yrs after prolonged 'hot' drought as ground water deficits lag precipitation deficits [36].
We can monitor declining canopy water content [97] and pre-emptively anticipate physiological tree drought stress by gaging within-crown thermal differences [45] or oxidative stress (PRI, xanthophyll cycle [98]), but the analysis presented here suggests that biological co-variates, bole growth vigor and tree-to-tree competition, aid in predicting declining BAI. Only the latter co-variate can be assessed remotely. Continuity of mature tree canopies (e.g., greater cover) may be assessed with LiDAR, which is increasingly monitored and available over forested regions. At the California state-wide level [79], a change metric was successfully identified using NDVI that was associated with increased tree mortality 6 to 19 months later [99]. Crown vigor was correlated to 'bole vigor' in an analysis of ponderosa pine [60] but was not a significant variant in our model predicting BAI, whether or not bole vigor was also included. The predictive capacity of our model dropped from moderate to poor when bole vigor was not concurrently considered.
Anticipating the meteorological conditions that will lead to extensive tree die-offs in mid-to upper elevations is needed, as loss of vegetation cover will lead to lower on-site water retention and 'flashier' downslope flows. DIs alone do not appear to be particularly effective in predicting tree bole growth, an attribute that declines as trees succumb. Conditions and thresholds have been defined that lead to hydraulic failure in trees [100], but the approach may not be scalable in mid to upper elevation forests that rely on regolithic water in the last third to half of the growing season. Estimating the rate of decline, year over year in NDVI [79,100], oxidative stress [99], and or crown or sub-crown thermal imagery [45] could be used as a harbinger of tree decline and mortality, and perhaps with sufficient lead time, for mitigation. We can anticipate loss of water resources with DIs, but at present we only monitor the degradation of mid elevation, open forest Sierran ecosystems with lower water availability and higher temperatures.

Conclusions
An 86 yr meteorological record was constructed for mid-elevation sites in the Sierra Nevada, useful for hydrologists, meteorologists, and biologists in future investigations. Correlations between on-site and modeled environmental data within a single DI were high, suggesting that modeled data, ubiquitously used today, can be reliably and interchangeably utilized in research applications in complex terrain. However, periodic discrepancies were discovered in the monthly mean temperature, and the estimated error variance of monthly summed PRISM precipitation exhibited heterogeneity. A similar study on the eastern slope of the Sierra Nevada also observed discrepancies in monthly temperatures, but for different time periods. Correlations among meteorological-and physics-based drought indices (SPEI and AI) and meteorological-, physics-, and pedologically-based drought indices (PDSI and derivatives) were high, suggesting that a similar outcome could be expected from the use of any one of these indices. Correlations of these indices with local streamflow were high, with the highest correlation provided by SPEI and AI. All current year DIs evaluated in this study were unable to predict NDVI in three sites with variable tree densities across this landscape. Drought indices (SPEI or PDSI) accompanied by biological covariates (bole growth vigor, tree-to-tree competition) were moderately predictive of BAI. Although DIs were significant, their effect sizes were low. DIs alone were poor predictors of BAI. This conclusion is robust without regard to which index or parameterization was used. A threshold of some level of meteorological-and physics-based DI such as SPEI or AI may aid predicting biological responses to multiple years of extreme drought but may not permit time for management to mitigate its effects at scale. Rates of change of remotely-detectable tree drought stress, crown vigor, and stand density (canopy continuity) may be a means of identifying high risk locations as well as providing lead time for management.