Satellite-derived foresummer drought sensitivity of plant productivity in Rocky Mountain headwater catchments: spatial heterogeneity and geological-geomorphological control

Long-term plot-scale studies have found water limitation to be a key factor driving ecosystem productivity in the Rocky Mountains. Specifically, the intensity of early summer (the ‘foresummer’ period from May to June) drought conditions appears to impose critical controls on peak ecosystem productivity. This study aims to (1) assess the importance of early snowmelt and foresummer drought in controlling peak plant productivity, based on the historical Landsat normalized-difference vegetation index (NDVI) and climate data; (2) map the spatial heterogeneity of foresummer drought sensitivity; and (3) identify the environmental controls (e.g. geomorphology, elevation, geology, plant types) on drought sensitivity. Our domain (15 × 15 km) includes four drainages within the East Water watershed near Gothic, Colorado, USA. We define foresummer drought sensitivity based on the regression slopes of the annual peak NDVI against the June Palmer Drought Severity Index between 1992 and 2010. Results show that foresummer drought sensitivity is spatially heterogeneous, and primarily dependent on the plant type and elevation. In support of the plot-based studies, we find that years with earlier snowmelt and drier foresummer conditions lead to lower peak NDVI; particularly in the low-elevation regions. Using random forest analysis, we identify additional key controls related to surface energy exchanges (i.e. potential net radiation), hydrological processes (i.e. microtopography and slope), and underlying geology. This remote-sensing-based approach for quantifying foresummer drought sensitivity can be used to identify the regions that are vulnerable or resilient to climate perturbations, as well as to inform future sampling, characterization, and modeling studies.


Introduction
Ecosystems in headwater catchments are important for water resources, because they influence hydrology through evapotranspiration (ET) and nutrient cycling (e.g. Lukas et al 2015, Maxwell andCondon 2016). Recent global-climate-model ensembles predict increased temperature and earlier snowmelt in western North America (Higgins andShi 2001, Diffenbaugh et al 2013). Additionally, some studies predict reduced spring precipitation and increased latesummer monsoon precipitation in the future (Seth et al 2011). Together, these changes would increase the length of time between snowmelt and summer monsoon, or the 'foresummer' part of growing seasons (Rauscher et al 2008, Swain andHayhoe 2015). Low snowpack years with earlier snowmelt would expose plants to potentially longer and drier periods before the onset of monsoonal precipitation. Combined with predicted increasingly warmer temperatures, this foresummer period could become more drought-like.
Recently, (Sloat et al 2015) documented the importance of this foresummer drought period by combining a watering manipulation experiment with 11 years of long-term monitoring data at the Rocky Mountain Biological Laboratory (RMBL) in Gothic, Colorado, USA. They found that peak and cumulative net ecosystem productivity (NEP) is negatively correlated with the severity of drought conditions in the primary growing season (June). They concluded that NEP will not increase in the future, despite the increase in temperature and longer growing seasons. This is consistent with other studies that reported water limitation of ecosystem productivity in the Rocky Mountain regions (e.g. Lamanna 2012, Williams et al 2012, Harte et al 2015 and in the western USA (Berner et al 2017). These regions are thus in contrast with other regions where water is not a limiting factor, and therefore early snowmelt lengthens the growing season, and increases the rate of peak and cumulative ecosystem productivity (e.g. Euskirchen et al 2006;Ernakovich et al 2014).
Here, we address the key challenge of scaling up such plot-scale experiments in order to quantify overall ecosystem and/or plant productivity at the scale of watersheds. Ecosystems in mountainous regions are particularly heterogeneous, influenced by steep and complex terrains. In these systems, plant types can vary on a spatial scale of 50-100 m (Zimmermann and Kienast 1999). Slope and aspect affect solar radiation, which in turn influences energy balance and soil moisture (Körner 2007). Soil moisture is also affected by plant types, soil types, and other factors (e.g. Mohanty et al 2000). In addition, snow accumulation and melting-which leads to infiltration and provides a critical water storage mechanism for ecosystems in the growing season ( (Anderson et al 2014, Painter et al 2016. In this study, we propose a new sensitivity metric for the Rocky Mountain region, and to map this sensitivity using historical satellite and climate data. The historical records of satellite datasets are valuable in evaluating the long-term trends and responses of ecosystems to climate variations, and also in inferring their future responses to changing climate (e.g. Zhao and Running 2010, Seddon et al 2016, Knowles et al 2017, 2018, Stocker et al 2019, Dong et al 2019. Given the water-resource limitation in this region during early growing seasons, herein we define foresummer drought sensitivity as the sensitivity of peak plant productivity to foresummer drought conditions. We then quantify this sensitivity based upon the historical records of the Landsat-derived normalized difference vegetation index (NDVI) at 30 m resolution, which is known to be strongly correlated with plant productivity (e.g. Tucker et al 1985, De Jong et al 2011, Dong et al 2019. We represent the drought condition based on the June Palmer Severity Drought index (PSDI) in the same manner as (Sloat et al 2015). In contrast to other satellite-based studies, our drought sensitivity is based on plot-scale experiments and the associated system understanding shown in (Sloat et al 2015). In addition, we use a machine learning approach to investigate environmental controls on spatially heterogeneous sensitivity, including elevation, geomorphology, and geology, using publicly available spatial datasets. Such data-driven analysis provides useful insights into underlying processes, enhances our ability to predict future trajectories, informs mechanistic ecohydrological models, and also facilitates site characterization and sampling plans.

Study area
We consider an approximately 15-km-by-15-km domain near Gothic, Colorado, USA (figure 1). (Hubbard et al 2018) provides a detailed site description. The domain is part of the Elk Mountain Range in the Rocky Mountains, with elevation ranging from 2800 m to~4000 m (figure 1(a)). The major land cover types are rock outcrop (12%), evergreen forest (29%), deciduous forest (18%), and grassland (30%; figure 1(b), NLCD 2011). The domain includes four drainages (East River, Washington Gulch, Slate River and Coal Creek) and four of the five experimental plots used in (Sloat et al 2015).
Historically, snow precipitation starts in October to November, and the first bare-ground date ranges from May to June. (Carroll et al 2018) analyzed the peak snow distribution in April, 2016 based on the NASA Airborne Snow Observatory, and found that the snow depth varies, ranging from 0 to 2.36 meters depending on the elevation, aspect and plant cover type. At the Butte Snow Telemetry (SNOTEL) station  (Sloat et al 2015), and the red square is the Butte SNOTEL location. The red region in (b) is the developed area.
(figure 1), the historical average of peak snow-waterequivalent and first bare-ground date are 400.5 mm and May 21st, respectively.

Palmer drought index and climate data
We used the June PDSI of Colorado Division 2 from NOAA (www.cpc.ncep.noaa.gov/products/analysis_ monitoring/cdus/palmer_drought/) to represent foresummer drought conditions. PDSI is computed based on precipitation, temperature, and division constants (such as soil water capacity). Although the limitations of PDSI have been recognized (Alley 1984, Dai et al 2004, Trenberth et al 2014, it is still the most widely used index for drought conditions (e.g. Dong et al 2019). Note that although the data source is different from (Sloat et al 2015), we assume that the general climate variability is captured similarly by both PDSIs. In addition, we evaluated other drought indices (text S1): Standardized Precipitation Index (SPI; Mckee et al 1993), and Standardized Precipitation Evapotranspiration Index (SPEI; Beguería et al 2010, Vicente-Serrano et al 2012. We also used snowmelt timing (i.e. first bareground date) and June mean air temperature from the Butte SNOTEL station (elevation 3097 m; www.wcc.nrcs.usda.gov/snow/). We used the homogenized SNOTEL temperature data provided by Oyler (et al 2015). All the data values are included in table S1 and figure S1 (available online at stacks.iop.org/ERL/15/084018/mmedia). We also confirmed that the average June precipitation is significantly lower than the other months (table S2).
These images were processed, including the atmospheric correction by the LEDAPS method (http://ledaps.nascom.nasa.gov/). We computed NDVI at each pixel, and annual peak NDVI (i.e. the maximum value at each pixel) in each year. Finally, we downloaded these annual peak NDVI images for further analysis. Since two Landsat paths overlapped over this domain, the repeat cycle was 8 d, which contributed to minimizing the effect of cloud coverage.
We first extracted peak NDVI at the pixels corresponding to the observation plots in (Sloat et al 2015) to investigate the relationship between peak NDVI and June PDSI, snowmelt timing, and June mean air temperature. We then defined the foresummer drought sensitivity as the slope of peak NDVI as a linear function of June PDSI. The slope represented the change in peak NDVI given the change in June PDSI. We also computed the average and standard deviation (SD) of annual peak NDVI at each pixel. In addition, we analyzed the relationship between NDVI and leaf area index (LAI) based on the groundbased measurements collected in 2019 (text S2 and figure S2).

Random forest analysis for environmental controls on drought sensitivity
We investigated key controls on foresummer drought sensitivity, based on other spatial data layers used in the hydrological modeling study within this domain (Pribulick et al 2016, Foster andMaxwell 2019). The Random Forest (RF) method is a machine-learning method developed by (Breiman 2001) to predict responses based on mixed numerical and categorical predictors, and to identify important predictors for given responses (Hastie et al 2001). RF generates a large number of regression trees from bootstrapped subsampled data, and averages over all the trees. RF is known to work well with correlated predictors similar to ridge regressions (Hastie et al 2001). In environmental applications, (Bachmair and Weiler 2012) used RF for identifying key controls on hillslope hydrological dynamics.
We defined a regression of foresummer drought sensitivity as a function of environmental variables. Using Topotoolbox (Schwanghart and Kuhn 2010), we computed topographic metrics based on the digital elevation model (DEM) from the National Elevation Dataset (30 m resolution, U.S (Geological Survey 2002)), including slope, topographic wetness index (TWI), bedrock-weighted upslope accumulated area (UAAB), and topographic position index (TPI). TWI is the log of flow accumulation area divided by slope, and TPI represents the local-scale variation of topography after topographic trend (i.e. the moving average of 100-m scale) is removed (Gillin et al 2015). Since solar radiation is known to be important for high-elevation mountain regions (Körner 2007), the annual sum of hourly potential solar radiation (including direct, diffuse, and reflected) was calculated from DEM, based on (Hebeler 2016) and (Kumar et al 1997).
For geology, we used the digitized geological map from the USGS National Geologic Map Data Base (Pribulick et al 2016). We grouped geological classes into six main classes: shale, igneous rock, alluvial, glacial, landslide, and unconsolidated deposits. Although a soil map was available, it was uniform except for outcrop regions and was thus not informative. We assumed that the geological map and topographic metrics could capture the variability in soil properties, since (Bailey et al 2014) and (Gillin et al 2015) documented the strong correlations between topographic metrics and soil properties.
Within the RF algorithm, the importance ranking of predictors was created by (1) setting aside a subset of data as a testing set (i.e. out-of-bag data), (2) predicting the drought sensitivity and computing the accuracy (i.e. out-of-bag error), and (3) computing the increase in the mean-squared-errors (MSE) of prediction after permuting each predictor (i.e. randomly assigning the predictor values from the data values). In addition, we created partial dependence plots to visualize the dependence of sensitivity on each predictor. We used R's randomForest package (cran.r-project.org/web/packages/rpart/index.html). The number of trees was equal to 800, which was enough to achieve convergence. The number of candidate variables at each split was the number of variables divided by three, and the minimum size of terminal nodes was five.

Results
At the plot locations in (Sloat et al 2015), Landsatderived peak NDVI is positively correlated with June PDSI and snowmelt timing (figures 2(a) and b, table S3), which is consistent with their findings for peak NEP (note that the data range of PDSI is different because of the differences in the PDSI sources). Increased drought conditions and earlier snowmelt are associated with decreased peak NDVI. Similar to peak NEP, the two subalpine-zone plots (elevation 3115 m and 3380 m) have higher peak NDVI than the montane-zone plots (elevation 2710 m and 2815 m). In addition, peak NDVI is negatively correlated with June mean temperature (June T) at all the locations (figure 2(c) and table S3). These findings are consistent when we use the other drought indices (SPI and SPEI; figure S3) and the linear regression without the extreme years ( figure S4). There are some differences between peak NDVI in figure 2 and peak NEP shown in (Sloat et al 2015). In figure 2, the slope values of peak NDVI are distinctly different between the subalpine and montane plots. At the subalpine plots, peak NDVI has lower dependency on June PDSI, snowmelt timing, and June T. By contrast, in figure 3 of (Sloat et al 2015), the slope values are similar at the four plots.
The average peak NDVI (figure 3(a)) is spatially heterogeneous over the domain, ranging from 0.2 to 0.9 in the vegetated area. The heterogeneity is related to both elevation and vegetation type (figure S5). Within the grassland area ( figure S5(a)), the overall trend of peak NDVI increases with elevation up tõ 3100 m, and then decreases. The deciduous forest region (i.e. Populus tremuloides or aspen) has higher average peak NDVI than the other vegetation types ( figure S5(b)). The evergreen forest has lower peak NDVI on average across the watershed (figure S5(c)), and also smaller spatial heterogeneity compared to the other vegetation types. Year-to-year variability of peak NDVI is represented by SD at each pixel (figure 3(b)). The region with the higher average peak NDVI (figure 3(a)) does not necessarily correspond to the one with high SD ( figure 3(b)). Higher elevation portions of the East River drainage, for example, have high peak NDVI on average, but low SD.
The foresummer drought sensitivity of peak NDVI (figure 3(c)) is positive in 94.1% of the vegetated area, although it is highly heterogeneous across the domain. Although the sensitivity map is fairly similar to the SD map ( figure 3(b)), the spatial heterogeneity of sensitivity is more pronounced than the SD. The southern (or lower elevation) part of the East River watershed (lower elevation and grassland areas) is particularly sensitive to the June drought condition. The northern (or higher elevation) part of the East River watershed has lower sensitivity, although the average peak NDVI is high ( figure 3(a)). The spatial heterogeneity of sensitivity depends heavily on plant types and elevation (figure 4). Summary statistics (table S4) show a clear dependency of drought sensitivity on plant types, confirmed by Tukey's pairwise comparison test (p-values < 1 × 10 −4 ). Grasslands ( figure 4(a)) have higher sensitivity than the other Figure 2. Landsat-derived annual peak NDVI as a function of (a) June Palmer Drought Severity Index (PDSI), (b) first bare-ground date, and (c) average June temperature at the long-term observation plots in (Sloat et al 2015). In (a)-(c), the line is based on linear regression. In (a), PDSI less than −4.0 is extreme drought, −3.0 to −2.0 is severe to moderate drought, −1.9 to +1.9 is normal, and +2.0 above is unusual to extreme moist conditions. The correlation coefficients are shown in table S3. plant types, particularly at lower elevation, and also has larger spatial heterogeneity across the domain. Elevation dependency (figure 4(d)) in the grassland region is much more apparent than in the SD map ( figure S6(d)). Evergreen forests exhibit the lowest sensitivity to the foresummer drought, although the sensitivity is still positive in 92.6% of the area. In addition, sensitivity is spatially less heterogeneous without significant elevation dependency (figure 4(f)).
We applied the RF analysis to foresummer drought sensitivity within the grassland region. Although we have the results in other plant types (table S5 and figures S7 and S8), we focus our discussion on the grassland region, because the grassland region has (1) higher spatial heterogeneity than other plant types, (2) the locations corresponding to the long-term plots in (Sloat et al 2015), and (3) the ground-based LAI-NDVI relationship (figure S2). The coefficient of determination (R-square) is 0.58, with the p-value less than 10 -15 . In the importance ranking (table 1), elevation is the strongest predictor for foresummer drought Table 1. Parameter importance ranking from the random forest analysis; the parameters influencing the spatial heterogeneity of foresummer drought sensitivity within the grassland region. The importance measure (i.e. %MSE) is normalized by one for elevation, so that it represents relative importance compared to elevation. The shaded cells indicate the top four in the importance ranking.  4(e)). Net potential radiation, topography position index (TPI), geology, and slope follow in the ranking. The three topographic metrics (TWI, UAAB and curvature) are relatively weak predictors. In addition, we have analyzed the datasets in different resolutions up to 600 m, which showed the same predictors as the 30 m resolution results (texts 3; table S6). Partial dependence plots are shown for the top four predictors in the importance ranking (figure 5). The dependency on elevation (figure 5(a)) is approximately linear, which is consistent with the elevation trend in figure 4(d). The dependency on net potential radiation (figure 5(b)) is nonlinear, with the effect more pronounced for the higher radiation regions. The dependency on TPI (figure 5(c)) is close to a step function, such that the regions having higher than the overall elevation gradients (i.e. microtopographically elevated) have higher drought sensitivity. With respect to geology (figure 5(d)), the igneous rock region is associated with decreased drought sensitivity, while glacial, landslide, and unconsolidated deposits are associated with increased sensitivity. We also investigated the correlations among the predictors such as elevation with radiation and aspect ( figure S9).

Discussion
At the long-term study plots in (Sloat et al 2015), the satellite observations of peak NDVI are consistent with peak NEP such that (1) peak NDVI is positively correlated with June PDSI and snowmelt timing, and (2) peak NDVI is greater at the subalpine plots than at the montane plots. These consistent responses confirm that plant dynamics are water limited in this region and that early snowmelt decreases plant productivity, as suggested by previous studies (Harte et al 2015, Sloat et al 2015. The subalpine plots are considered less water limited, given deeper snowpack and later snowmelt. In addition, we find that peak NDVI is negatively correlated with average June temperature. Higher June temperature is known to be associated with earlier snowmelt and higher ET , which exacerbates the foresummer drought condition and has a negative impact on plant growth. There are differences between peak NDVI and NEP. While peak NEP responds similarly to June PDSI across the elevation gradient in (Sloat et al 2015), satellite-derived peak NDVI is less sensitive at the subalpine plots. This could result from the fact that NDVI represents only aboveground plant dynamics, while NEP includes soil respiration. (Sloat et al 2015) found that soil respiration was less affected by watering experiments, suggesting that soil respiration was less sensitive to droughts. Although NDVI has been used for upscaling NEP (e.g. Sturtevant and Oechel 2013), the applications in mountainous regions may not be straightforward, owing to nonlinear responses along the elevational gradient.
We considered the potential effect of the NDVI saturation at the high LAI region. Although this grassland region is not a high biomass region (Gao et al 2000, Huete et al 2002, Gu et al 2013, NDVI at the two high elevation locations (figure 2) are as high as 0.85. In the NDVI-LAI relationship (figure S2), we observe possible saturation in NDVI above~0.85. We fitted these datasets with a linear and second-order polynomial function for NDVI < 0.85. Since the R 2 and BIC are comparable, we may use the linear function up to NDVI = 0.85. We would note that peak NDVI is less than 0.85 in 92.4% of the domain even in 1995, when peak NDVI is the highest. In parallel, we considered the potential effect of such saturation on our results. If the effect of the decreased sensitivity were due to NDVI saturation, the sensitivitydefined as the slope of peak NDVI as a function of PDSI-would decrease in high NDVI regions. In figure 4(d), the sensitivity decreases as the elevation increases, while the average peak NDVI (figure S5(d)) also decreases above~3100 m. Therefore, we may conclude that the decreased sensitivity at high elevation does not result from the saturation effect. In addition, we examined the correlation between sensitivity and average peak NDVI, which was not significant (correlation coefficient of −0.026).
Satellite-derived NDVI allows us to extend our plot-scale understanding to the watershed scale. Landsat 5 has provided long-term historical images at high-enough resolution to distinct plant types and gauge the effect of topography and geology on sensitivity. In the subalpine zone (around 3000-3300 m), peak plant productivity is high on average, and has small interannual variability without significant dependency on the regional-scale June PSDI, possibly because this zone has more snow and water compared to lower elevations. In addition, our results show that the magnitude and spatial variability of drought sensitivity is clearly dependent on plant type, with the grassland regions having higher sensitivity and higher spatial heterogeneity. This dependency on plant type is considered to result from rooting depth as well as geographic location-evergreen trees tend to occupy north-facing slopes that have more snow accumulation and higher soil moisture. At the same time, foresummer drought sensitivity is predominantly positive within the evergreen forest regions, suggesting that the increased drought severity-due to early snowmelt and/or increased spring temperature-would decrease forest productivity, which is consistent with a basin-scale study by (Knowles et al 2018). In addition, we observe that 6% of the region, primarily at high elevation, has increased peak NDVI in drought years. It suggests that the high elevation regions are temperature-limited rather than water-limited. This is consistent with (Dong et al 2019), which found that MODIS-based NDVI increased at higher elevation in drought years.
In this study, we defined the foresummer drought sensitivity by the slope of the linear regression between peak NDVI and June PDSI. With respect to sensitivity measures, there are slope-based and variance-based measures in general to represent sensitivity (e.g. Morris 1991, Saltelli et al 2008, Wainwright et al 2014. Several studies have used variance or variance-based metrics to represent ecosystem sensitivities (e.g. Seddon et al 2016). In our results, we find that the slope-based sensitivity measure is more informative in this type of analyses, since we can identify positive or negative changes associated with foresummer drought conditions.
We consider that PDSI and other drought indices (such as PSI and SPEI) represent the regional-scale climatic variability being driven by precipitation and temperature. Vincente-Serrano et al (Vicente-Serrano et al 2012) found that ecosystem responses (i.e. tree-ring growth and wheat yield) to droughts were captured by SPEI, as well as other drought indices. In our analysis, the response of peak NDVI was consistent to each other among these three drought indices (PDSI, PSI, SPEI), confirming the impact of water limitation on plant productivity over this region. At the same time, we highlight that this study focuses on the spatial variability of foresummer drought sensitivity at the local-scale (30 m), so that we can resolve the effect of topography, plant type, and geology. While these local-scale environmental characteristics could be viewed as secondary factors, recent studies have found that such characteristics (e.g. geology) are important for subsurface water storage (Markovich et al 2016) and the resilience of ecosystems (Rempe and Dietrich 2018).
The RF analysis enables us to identify key environmental controls on foresummer drought sensitivity within each plant cover type. Elevation and net radiation are the two most dominant factors, possibly because they control surface energy balance and snow accumulation and melting. We would note that the higher-elevation hillslopes tend to be north-facing in our domain, which could amplify the effect of reduced drought sensitivity in high-elevation regions. The topography position index (i.e. indicator of microtopography) and slope are known to control soil moisture (Mohanty et al 2000, Gillin et al 2015. Falco et al (2019) found a significant correlation between slope and soil moisture in the grassland regions within this domain. In addition, the results show the importance of underlying geologic composition on drought sensitivity. Well-drained soil developed upon glacial and landslides deposits are likely to shed near-surface soil moisture (associated with plant rooting) rapidly after snowmelt. In contrast, the region underlain by shale and igneous rocks has lower drought sensitivity. Rempe et al (2018) found that fractured bedrock can retain water during droughts and is less affected by year-to-year variability. Having shallow bedrock may provide resilience to droughts.
This study used publicly available PDSI and Landsat data to estimate foresummer drought sensitivity of peak plant productivity in headwater catchments. We did not explicitly include other datasets, for example, the spatial distribution of snow accumulation/snowmelt and precipitation, since these factors are difficult to map in space and time (Lettenmaier et al 2015). Instead, we assumed that the topographic metrics are reflective of snow and precipitation patterns, given that the effects of topography on these patterns have been well documented in many studies (e.g. Anderson et al 2014). Both these assumptions, and our approach, open the door for upscaling plot-scale analyses and understandings to a large area, using publicly available datasets. We acknowledge that the overlapping coverage of Landsat paths was advantageous for our domain to minimize the impact of cloud coverage. At the same time, our analysis based on different spatial resolutions found that the effects of key drivers (i.e. elevation and radiation) were consistent up to the resolution of several hundred meters, which would suggest that we may use lower-resolution high-frequency satellites such as MODIS. Remote-sensing-derived drought sensitivity can be a useful metric for identifying the regions that are resilient or vulnerable to climate perturbations and long-term climatic shifts, as well as for identifying key underlying processes.

Acknowledgments
This material is based upon work supported by the U.S. Department of Energy

Data Availability Statement
The data that support the findings of this study are openly available. All the datasets in this study are publicly available through the data source specified in the manuscript.