Drought Sensitivity and Trends of Riparian Vegetation Vigor in Nevada, USA (1985–2018)

Dryland riparian areas are under increasing stress due to expanding human water demands and a warming climate. Quantifying responses of dryland riparian vegetation to these pressures is complicated by high climatic variability, which can create strong, transient changes in vegetation vigor that could mask other disturbance events. In this study, we utilize a 34-year archive of Landsat satellite data to (1) quantify the strength and timescales of vegetation responses to interannual variability in drought status and (2) isolate and remove this influence to assess resultant trends in vegetation vigor for riparian areas across the state of Nevada, the driest state in the USA. Correlations between annual late-summer Normalized Difference Vegetation Index (NDVI) and the Standardized Precipitation–Evapotranspiration Index (SPEI) were calculated across a range of time periods (varying timing and durations) for all riparian pixels within each of the 45 ecoregions, and the variability of these values across the study area is shown. We then applied a novel drought adjustment method that used the strongest SPEI–NDVI timescale relationships for each ecoregion to remove the influence of interannual drought status. Our key result is a 30 m resolution map of drought-adjusted riparian NDVI trends (1985–2018). We highlight and describe locations where impacts of invasive species biocontrol, mine water management, agriculture, changing water levels, and fire are readily visualized with our results. We found more negatively trending riparian areas in association with wide valley bottoms, low-intensity agricultural land uses, and private land ownerships and more positive trends in association with narrow drainages, public lands, and surrounding perennial water bodies (an indication of declining water levels allowing increased vegetative cover). The drought-adjusted NDVI improved the statistical significance of trend estimates, thereby improving the ability to detect such changes. Results from this study provide insight into the strength and timescales of riparian vegetation responses to drought and can provide important information for managing riparian areas within the study area. The novel approach to drought adjustment is readily transferrable to other regions.


Introduction
Across the globe, riparian ecosystems provide many ecosystem services and are valued for their role as critical hubs for the processing and transport of energy and matter, the biodiversity they support, and their scenic, economic, and recreational values [1][2][3]. In arid and semi-arid regions, riparian areas are especially valuable, but are sometimes points of conflict as competition for water between natural systems and humans increases due to limited supplies [4]. Hydrologic changes due to surface and groundwater development, alteration of stream channels and wetlands, land development, and non-native species introductions, in combination with changing climate, have increased pressure on wetland and riparian systems in arid and semi-arid regions in particular, often resulting in declines in their condition and extent [1,4,5].
As the driest state in the U.S., demand for water in Nevada is particularly high and supports local economies and communities, irrigation of crops, mining, and municipal uses, in particular [6]. Long-term declines in groundwater levels [7] and streamflow [8] have been observed in some places, reflecting the effects of groundwater and surface water development, reductions in snowpack, and increased evaporative demand due to drought and warming climate. The demand for water in the state will increase with population and as the climate warms, highlighting the need to quantify the status and trends of aquatic and riparian ecosystems so that future water allocation and management decisions can be based on the best available science.
Nevada's geographic position in the rain shadow of the Sierra Nevada results in a climate that is highly variable and precipitation amounts are strongly dependent on whether large frontal winter storms, termed atmospheric rivers, are able to penetrate over or around this mountain range [9]. The presence or absence of these large storms in a given year result in high precipitation variability relative to the average conditions [10], and is a dominant influence on the duration and magnitude of droughts [11]. Due to the water-limited hydroclimate in Nevada, there is high interannual variability in vegetation productivity [12], which can make it very difficult to isolate the impacts of land management changes or disturbances from variability due to interannual climatic fluctuations or drought. While the productivity and vigor of upland vegetation in water-limited regions is often strongly related to variability in the annual water balance [13,14], the role of interannual climate variability as a determinant of vegetation productivity in riparian areas is less clear given the potential for these systems to access and use water contributions from precipitation, surface water, and shallow groundwater, all of which vary across different spatial and temporal scales [15]. Separating the influences of interannual climatic variability from management, disturbance, or other human activities requires long-term monitoring data [16], which are generally not available due to the high costs and extensive areas that need to be monitored. Recent advances in cloud computing are helping to address these challenges by enabling efficient processing and analysis of large satellite remote sensing and climate datasets over multiple decades [17,18]. With a record of over 35 years, a spatial resolution of 30 m, and a temporal resolution of approximately 2 weeks, Landsat satellite data are unparalleled in terms of their value for characterizing long-term changes in the narrow riparian areas of arid and semi-arid regions [19].
In this study, we access the Landsat and climate data archives in the Google Earth Engine [17] to assess drought sensitivities and how they may obscure trends in riparian vegetation caused by other factors across the state of Nevada. Specifically, our objectives are to quantify the strength and timescale of responses of riparian vegetation to interannual variability in drought conditions and to use this information to remove interannual variability due to drought from longer term trends in vegetation that may be associated with disturbances or human impacts over a 34-year time period.

Study Area
Our study area consists of riparian areas throughout the state of Nevada, USA. Nevada is predominantly within the Northern, Central, and Mojave Basin and Range Ecoregions of North America, with a small portion covering eastern Sierra Nevada. The terrain is dominated by parallel north-south mountain ranges separated by arid valleys associated with horst and graben faulting of the Basin and Range physiographic region. The Sierra Nevada mountains along the western edge of the state create a rain shadow that makes Nevada the driest state in the country. Annual average precipitation and the proportion falling as snow varies with latitude and altitude, ranging from about 100 mm in the lowlands of the south to 250 mm in the north. Precipitation at higher altitudes may exceed 1 m per year, with a large proportion falling as snow.
Riparian areas are defined based on a 30m resolution Map of Riparian Vegetation for Nevada (MRVN; Figure 1; https://www.dri.edu/project/mrvn/; [20]). Riparian areas were mapped by an algorithm that identified where anomalously high NDVI occurred for a given elevation and aspect. Most riparian vegetation occurs at high elevations and in the northern portion of the study region ( Figure 1) where precipitation and surface water availability are highest. Dominant trees and shrubs include mesic-adapted species such as Fremont, black, and narrowleaf cottonwood (Populus fremontii, P. trichocarpa, P. angustifolia.), aspen (Populus tremuloides), and multiple willow species (Salix spp.). Common herbaceous species include grasses, rushes, and sedges as well as low-intensity agricultural pasture and hayland species [21] that are commonly intermingled with natural riparian vegetation and respond to water availability similarly to their natural surroundings [20]. Riparian areas are less commonly found in the Mojave region, which includes the southern and more arid parts of the state. In this region, riparian vegetation is often associated with localized sources of groundwater and typically includes xeric-adapted phreatophyte species such as honey and screwbean mesquite (Prosopis glandulosa, P. pubescens.), tamarisk (Tamarix spp.), and velvet ash (Fraxinus velutina). Riparian areas are defined based on a 30m resolution Map of Riparian Vegetation for Nevada (MRVN; Figure 1; https://www.dri.edu/project/mrvn/; [20]). Riparian areas were mapped by an algorithm that identified where anomalously high NDVI occurred for a given elevation and aspect. Most riparian vegetation occurs at high elevations and in the northern portion of the study region ( Figure 1) where precipitation and surface water availability are highest. Dominant trees and shrubs include mesic-adapted species such as Fremont, black, and narrowleaf cottonwood (Populus fremontii, P. trichocarpa, P. angustifolia.), aspen (Populus tremuloides), and multiple willow species (Salix spp.). Common herbaceous species include grasses, rushes, and sedges as well as low-intensity agricultural pasture and hayland species [21] that are commonly intermingled with natural riparian vegetation and respond to water availability similarly to their natural surroundings [20]. Riparian areas are less commonly found in the Mojave region, which includes the southern and more arid parts of the state. In this region, riparian vegetation is often associated with localized sources of groundwater and typically includes xeric-adapted phreatophyte species such as honey and screwbean mesquite (Prosopis glandulosa, P. pubescens.), tamarisk (Tamarix spp.), and velvet ash (Fraxinus velutina).

Landsat Data Processing
Interannual changes in the vegetation status of riparian areas were measured using normalized difference vegetation index (NDVI; [22]) data from Landsats 5 and 8, covering the period 1985 to

Landsat Data Processing
Interannual changes in the vegetation status of riparian areas were measured using normalized difference vegetation index (NDVI; [22]) data from Landsats 5 and 8, covering the period 1985 to 2018. The years 2011 and 2012 were excluded from this analysis due to the failure of the Thematic Mapper (TM) sensor on Landsat 5 and data dropouts from the scan line corrector mirror failure of the Enhanced Thematic Mapper on Landsat 7. A change in spectral bandwidths between the Landsat 5 TM sensor and the Landsat 8 Operational Land Imager (OLI) affects the continuity of NDVI measurements, and this was addressed using the adjustment to OLI NDVI values described in Huntington et al. [7]. The NDVI value for each year was computed as the median of cloud-free, good quality land surface reflectance values acquired between July 1 and September 30, with cloud and quality screening based on per-pixel image metadata. This represented a period when streamflow persistence [23] and mesic-associated vegetation [24] in this region tend to be most sensitive to interannual variability in precipitation, and when riparian vegetation undergoes slow phenological change due to the perennial availability of shallow groundwater. This seasonal window also was advantageous for minimizing the potential influence of cloud cover and use of the median value suppressed the influence of infrequent summer precipitation events or Landsat data artifacts such as sun glint anomaly and data loss [20].

Deriving Drought-NDVI Relations at the Ecoregional Scale
Annual drought effects on vegetation can vary significantly over space or over time depending on weather patterns, plant species adaptations and ecological interactions, in addition to human or natural disturbances that affect water availability [25]. To avoid convolving the effects of local disturbances on the relationship between interannual drought status and NDVI, we established a drought-NDVI relationship based on the NDVI signal for all riparian areas within each of the 45 Level IV ecoregions within the state of Nevada, as defined by the U.S. Environmental Protection Agency (Figure 1; [26]). Annual median late-summer (July-Sept) NDVI for each year in the archive were averaged across all riparian pixels within each ecoregion and then transformed into a percentage change from the mean of the period-of-record (1985-2018). As seen in Table 1, the area of ecoregions varied greatly, possibly raising concerns with the modifiable areal unit problem [27]. In this analysis, sample size could affect the stability of the estimates for ecoregion means. The smallest region contained 11 Landsat pixels (13ac in Table 1) which we judged to provide a sufficiently stable estimate of mean for the trend analysis. The EPA ecoregions were designed to provide a high level of ecological homogeneity, and this should provide a level of consistency in ecological response regardless of geographical extent of the ecoregion.
We quantified relationships between each ecoregion's annual percent change in NDVI and Standardized Precipitation-Evapotranspiration Index (SPEI; [28]); a multi-scalar drought index representing the difference between precipitation (PPT) and potential evapotranspiration (PET). We selected this index because preliminary analyses indicated that NDVI was more strongly related to the difference between PPT and PET than to PPT or PET alone (see Supplementary Figure S1) and because strong relationships have already been found in Nevada for SPEI compared to hydrologic droughts (i.e., streamflow and reservoir levels [23]) and compared to wildland fire potential and dead fuel moisture [29]. PPT and PET (based on Penman Monteith reference evapotranspiration [30]) data were obtained from the 4 km spatial resolution GRIDMET dataset [31], which is based on a sparse observation network in Nevada and has the potential for error at the weather scale such as light (drizzle) or extreme precipitation events [32] and temperature inversions that are common in complex terrain [33]. We acknowledge this limitation on the accuracy of GRIDMET. However, using climatic time scales (one month and longer) and standardized values reduces uncertainty in GRIDMET as regional climatologies are well captured. Moreover, NDVI, PET, and PPT exhibit a complementary relationship according to well-known complementary theory where atmospheric water demand, water use, and vegetation vigor are inherently linked [7,18]. Prior to analysis, we resampled the GRIDMET to 30 m resolution with bilinear interpolation to match the resolution of Landsat. Similar to the treatment of NDVI data, we averaged the PPT and PET across all riparian pixels within each ecoregion. The use of multi-scalar drought indices such as the SPEI provide opportunities to simultaneously characterize multiple drought attributes that are most relevant to vegetation responses, including consideration of both precipitation and atmospheric water demand, as well as their timing, magnitude, and duration [27]. Iterative analyses that assess relations between NDVI metrics and SPEI calculations at varying timescales have been used to identify the duration and timing of drought that vegetation is most sensitive to, which vary among vegetation types and across gradients of climatic conditions [13,[33][34][35]. Similar to these studies, we performed exploratory data analysis to determine the duration and timing (i.e., temporal scaling) of SPEI that is most predictive of the interannual NDVI response over the 1985-2018 time period. This involved iterating regressions between late summer NDVI and one-to twelve-month SPEI aggregation periods with SPEI end months (lag) ranging from Oct-Sept of the same water year as the NDVI. For each ecoregion, we calculated the difference between PPT and PET (PPT-PET) for each end month and aggregation period combination ( Figure 2). We then applied a non-parametric standardization to these values [34] using the SDAT MATLAB toolbox [35] and regressed them against the observed percent difference in NDVI from its long term (1985-2018) mean across that ecoregion. Since the difference from mean and SPEI are both zero-centered distributions, this was the same as fitting a slope between the two variables as a multiplicative effect. We used a least squares bisector [36] Type II regression [37], which is an "error in variables" method and is appropriate when there is measurement uncertainty in the independent (climate) variable. Sources of uncertainty include the fact that the climate data is modeled, bias corrected, and downscaled to interpolate point data and SPEI incorporates an estimate of evaporative demand. In this circumstance, the assumption of all errors being attributed to the dependent variable in ordinary least squares regression would result in a deflated estimate of the regression slope [38,39].
Following the iteration of these regressions, we identified the timing and duration (temporal scaling, e.g., February-May) of SPEI with the strongest relationship to NDVI for a given ecoregion and used this value in subsequent analyses ( Figure 2). In addition, we assessed variations in the strength and optimal scaling of the SPEI-NDVI regressions among ecoregions and as a function of ecoregional aridity (measured as the average annual water balance) to compare our results to those obtained in a global study of vegetation responses to drought time-scales [13].

Quantifying Trends in Drought-Adjusted NDVI
We applied the ecoregional SPEI-NDVI relationship at the pixel scale to derive a drought-adjusted NDVI ( Figure 2). This was accomplished by dividing the original NDVI values for each pixel by one plus the predicted SPEI effect for that year and ecoregion. Because the SPEI-NDVI relationship is modeled as a multiplicative effect on percent change, the application of this relationship at the pixel scale translates to a percentage change in NDVI which accounts for the fact that more densely vegetated pixels with higher NDVI have the capacity to exhibit a larger response to drought than more sparsely vegetated pixels.
We applied the non-parametric Mann-Kendall trend test [40,41] to timeseries of drought-adjusted NDVI values to assess the presence and significance of monotonic trends from 1985-2018 on a per-pixel basis ( Figure 2). Sen's slope estimator [42] was used to quantify the magnitude of trends and a threshold of p < 0.05 was applied to determine statistical significance. We also quantified trends in ecoregional SPEI using these metrics to provide a context for interpreting the drought-adjusted results. The direction and magnitude of trends in SPEI over the course of the study period, in part, determine the degree to which drought-adjusted trends in NDVI will differ from trends in the original NDVI time series. For example, in cases where SPEI is trending positively (i.e., toward wetter conditions) over time, the drought adjustment will reduce annual NDVI values during those wetter years to neutralize the effects of wetter conditions and provide an estimate of the expected NDVI in an average year based on the ecoregional SPEI-NDVI regression relationship. Overall, this results in a drought-adjusted trend that is less positive in magnitude relative to the trend slope of the original, unadjusted NDVI time series. Conversely, in cases where SPEI is trending negatively (toward drier conditions), the drought adjustment will result in a positive adjustment to the trend in annual NDVI values to neutralize the effect of drier conditions. To demonstrate and evaluate the effect of the drought adjustment, we compared differences in the slope magnitudes and statistical significance of original and drought-adjusted NDVI trends for the full time period (1985-2018), during which there was generally no trend in SPEI, and for a demonstration time period during which there was a significant decreasing trend (i.e., wetter to drier) in SPEI (1995-2013) across much of the study area. This latter comparison was made to highlight the utility of the drought-adjustment approach for accurately quantifying vegetation trends amid trending climatic conditions.
We assessed spatial patterns of trend direction and magnitude for the full 1985-2018 time series based on the distributions and central tendencies of slope values in relation to several landscape-scale datasets, including (1) riparian land cover types and land use that were defined using the Landfire Existing Vegetation Type classification (Group level) [21]; (2) land ownership using the Federal surface management dataset [43]; (3) annual precipitation [31]; (4) major landforms [44], and (5) hydrographic types, including perennial and intermittent streams and water body types from the NHDPlus dataset [45]. In the latter case, riparian vegetation was associated with NHD hydrographic types as a function of distance, with vegetation within the perimeter of NHD water bodies being associated with the water body type as the first priority, vegetation within 50 m of perennial streams being associated with perennial streams as the second priority, vegetation within 10 m of intermittent streams associated with intermittent stream types as the third priority, and vegetation at larger distances from NHD features designated as 'other.'

Drought Timescale-NDVI Relations by Ecoregion
Riparian NDVI responses to SPEI differed substantially across the region in terms of the amount of variability explained, sensitivity, and time scale (Table 1; Figure 3a-d). Coefficient of determination (R 2 ) values, representing the fraction of interannual variability in NDVI explained by SPEI, ranged between 0.03 and 0.82. R 2 s were higher in the Central and Northern Basin and Range ecoregions relative to the Mojave ( Figure 3e) and varied with the long-term average annual water balance (R 2 =0.28; Figure 4). The slopes of the SPEI-NDVI relations, signifying the percent change in NDVI, given a one-unit change in SPEI (i.e., sensitivity), ranged between −0.04 and −0.37. With a few exceptions, ecoregions in the southern part of the state tended to have lesser slopes, while the northern part of the state exhibited a mix of greater slopes corresponding with mid-and low-elevation basin, valley, and desert ecoregions and lesser slopes corresponding with higher-elevation forested ecoregions (Figure 3b).
Optimal SPEI aggregation timescales were greater than 7 months for most ecoregions, though the timing and aggregation periods varied among Level 3 ecoregions (Table 1; Figure 3c,d). For example, nearly all Northern Basin and Range ecoregions had optimal aggregation timescales greater than 8 months, with the start months corresponding closely with the start of the water year (i.e., in October or November). Most Mojave ecoregions had optimal aggregation timescales in the range of 6-9 months, with wide variation in start and end months. July was consistently the optimal end month for higher elevation Mojave ecoregions, while lower elevation Mojave ecoregion end months ranged between December and May. Optimal aggregation timescales were consistently 7 months ending in May for all Sierra-influenced ecoregions. There was wide variation among the Central Basin and Range ecoregions with aggregation durations ranging between 3 and 12 months and end months ranging between April and September, but most started between October and December and ended in May or June. In contrast to results from Vicente-Serrano et al. [13], there was essentially no relation between aggregation timescale and long-term average annual water balance (R 2 =0.02; Figure 3e, Figure 4). There was little trend in optimal time scale SPEI over the duration of the study period, with statistically significant trends occurring in only two ecoregions (13d and 80e; Table 1). The optimal time scale SPEI exhibited slightly negative trends for most ecoregions, indicating a tendency toward drier conditions, but trends were slightly positive in the Northern Basin and Range and other ecoregions in the northern part of the state (Table 1; Figure 3F).    Standardized Precipitation-Evapotranspiration Index (SPEI) timescale and late summer NDVI (as % difference from long term average) and (b) the optimal aggregation timescale at which the correlation between SPEI and NDVI is maximized.  Standardized Precipitation-Evapotranspiration Index (SPEI) timescale and late summer NDVI (as % difference from long term average) and (b) the optimal aggregation timescale at which the correlation between SPEI and NDVI is maximized.

Drought-Adjusted Trend Results
Overall, 22.8% of all riparian pixels had significantly positive drought-adjusted trends, 13.5% were significantly negative, and 63.7% had no statistically significant trend. We used results to identify 'hotspots' of change in the study region ( Figure 5). For example, the trend map revealed impacts in some areas that could readily be attributed to disturbances such as mine dewatering (Figure 5b), vegetation defoliation by the tamarisk leaf beetle (Figure 5e, 5h), fire (Figure 5f), and low-intensity agricultural land uses (Figure 5a). Of all statistically significant pixels, 3% were within 10 km of known occurrences of the tamarisk leaf beetle [46], 11% were within mapped fire perimeters [47], and 17% were riparian pixels classified as 'Agricultural Pasture and Hayland' [21]. Other locations along major river corridors and adjacent to lakes and reservoirs also showed large negative (Figure 5a, 5d) and positive changes (Figure 5c) that did not appear to be associated with the above-mentioned disturbances. In these latter cases, possible reasons for these observed trends are discussed in Section 4.2.
Broader scale patterns of change in riparian vegetation in association with long-term average climate, geomorphology, land ownership, and land cover were often subtle, but were observed between drainage channels (narrow valley landforms that occur in rugged terrain) and plains (flatter areas comprised of broad valley bottoms and plateaus), with the former exhibiting positive central tendencies (mean, median) in slope values, and the latter showing negative central tendencies in slope values and wider ranges of values, indicated by higher standard deviations (Figure 6a, Table 2).
Because precipitation covaries with landform, we made comparisons between low and high precipitation (using 12 cm annual precipitation, a natural break in the data, as the threshold), for each landform. This comparison revealed little difference in the distribution and central tendencies between these two precipitation categories (Figure 6a, Table 2), indicating that trends vary more with landform than with annual precipitation amounts.
Riparian areas classified by Landfire as 'Introduced Vegetation' and 'Agricultural Pasture and Hayland' had negative central tendencies, with the former having a relatively high standard deviation and the latter having the most negative central tendency (Figure 6b, Table 2). Among the major land ownership types within the state, Bureau of Land Management, U.S. Forest Service, and Other non-private land ownerships had positive central tendencies, while Private land ownerships had negative central tendencies (Figure 6c). Private and Other land ownerships had larger standard deviations relative to Bureau of Land Management and U.S. Forest Service ownerships (Figure 6c).     Differences in the distributions of slope values for riparian vegetation associated with different types of streams, rivers, and water bodies were evident (Figure 6c). Slope distributions for perennial water bodies such as lakes/ponds and swamp/marshes had the largest standard deviations and positive central tendencies, while central tendencies for intermittent water bodies and playas were negative. Intermittent and perennial streams/rivers had smaller standard deviations relative to water bodies, with perennial rivers having a larger standard deviation than intermittent. Both types of streams/rivers were skewed towards more positive slope values.

Comparison of Raw and Drought-Adjusted NDVI Trends
The drought adjustment resulted in an overall increase in slope magnitudes, making some negative slopes more negative, in particular ( Figure 7). This occurred primarily in the northeast corner of the state (Figure 8), where SPEI exhibited a slight wetting trend over the course of the study period (Figure 3f). In the south and central-eastern part of the state, where SPEI exhibited a slight drying trend over the course of the study period (Figure 3f), the drought adjustment resulted in more positive slope magnitudes (Figure 8). The drought adjustment also resulted in a higher frequency of pixels with statistically significant results due to an overall lowering of p-values (Figure 7). While 36% of all riparian pixels had statistically significant drought-adjusted trends (22.8% positive and 13.5% negative), only 30% of pixels (19% positive and 11% negative) had significant trends based on the raw, unadjusted time series. The largest differences between the raw and drought-adjusted trends occurred in the northeast corner of the state (Figure 8), where the relationship between SPEI and NDVI was strongest (Figure 2).  Because there was little trend in SPEI during 1985-2018, the differences between the raw and drought-adjusted trends over the entire time series were relatively small. In contrast, the differences between the raw and drought-adjusted trends for the 1995-2013 demonstration time period during which there was a decreasing trend in SPEI over much of the study area that was quite large (Figure 8). Moreover, in the absence of the drought adjustment, trends are consistently negative across the state due to the drought effect, and localized areas of known anthropogenic disturbances ( Figure 5) are less distinguishable from strong and variable trends in surrounding areas (Figure 8).  18year (1995-2013) subset are presented here to demonstrate the value of the drought adjustment concept, with the acknowledgment that this time period is insufficient to reliably quantify trends.
The optimal aggregation timescale for SPEI also varied considerably across the study region. Based on the range of annual water balance values for the study area, aggregation timescales for upland vegetation may be expected to range from four to nine months [13], with vegetation in more arid regions responding to shorter aggregation time scales. We did not observe this pattern in this riparian study; rather, the optimal SPEI aggregation period was relatively consistent across the aridity gradient and was greater than seven months for ¾ of the ecoregions, which covered most of the state. Thus, in more arid ecoregions in particular, the aggregation timescale for riparian areas is longer relative to that reported for upland vegetation with similar aridity [13], likely reflecting the greater water availability in riparian systems, although we acknowledge that the differences between  18-year (1995-2013) subset are presented here to demonstrate the value of the drought adjustment concept, with the acknowledgment that this time period is insufficient to reliably quantify trends.

Drought Timescale-NDVI Relations by Ecoregion
While several studies have examined the strength and optimal timescales of relations between drought metrics and NDVI in semi-arid regions at the global scale [13], and in the western US [48][49][50] specifically, these studies have been focused on upland vegetation communities. To our knowledge, this study is the first to identify the drought time scales most relevant to riparian systems in our study region, which, despite being controlled by the same climatic regime, exhibit substantially different responses to drought given differences in hydrology and water availability relative to surrounding upland regions.
The strength of the relationship between SPEI and NDVI varied substantially across the study area, with stronger relationships, indicated by R 2 values, observed in more northerly ecoregions. This contrasts with climate correlations for upland areas of Nevada from McGwire et al. [49] where correlations were very strong across central and southern areas and weaker in the northeast. At a global scale, the strength of the relationship between SPEI and NDVI varies across gradients of aridity, measured as the annual water balance (water year PPT-PET), with arid and semi-arid systems generally exhibiting stronger relationships [13,51]. Our study offered an opportunity to test the applicability of this relationship to riparian vegetation across a large gradient of annual water balances ranging from 500 (PPT > PET) to nearly −2000 (PET > PPT). Interestingly, we observed the opposite relationship from [13,51], with more arid regions exhibiting smaller R 2 values and less arid regions exhibiting larger R 2 values. While we acknowledge several differences between our study and these, including differences in the NDVI metric (maximum NDVI in [13], average growing season NDVI in [51] vs. July-Sept, representing post-peak NDVI used here), scale of the study, and differences in climate patterns and species compositions, all of which merit further investigation, we posit some additional hypotheses for these differences here. First, we note that riparian areas in the more arid southern (Mojave) region of the state tend to be rarer and have smaller areal extents. The weak relationship could be the result of an inability to detect the vegetation signal due to the overall rarity and small size of these areas. Conversely, differences in vegetation composition and potentially stronger groundwater dependency of riparian vegetation in the Mojave region, which is indicated by a lack of perennial streams on the landscape, may be driving the overall lower sensitivity of these areas to interannual variability in drought status. Some of the desert wash riparian features mapped in the south were mesquite and tamarisk, both of which are xeric phreatophytes that tap groundwater and tend to have low sensitivities to interannual climate variability [52,53]. Additional support for this latter reasoning includes the vastly different timings of peak NDVI between upland and riparian vegetation in Mojave ecoregions as compared to others in the state (Supplementary Figure S2). In the Mojave ecoregions, peak NDVI for upland vegetation typically occurs in April, closely following peak precipitation in winter, while riparian NDVI peaks in September, 5 months later. In contrast, peak riparian NDVI exhibits much shorter lags, on the order of 1-2 months, behind upland vegetation in other ecoregions in the state. The much later timing of peak NDVI in riparian vegetation as compared to the winter-dominated precipitation climatology and green-up of upland systems in the Mojave is consistent with the hypothesis that these ecosystems are not as reliant on precipitation and generally have access to consistent groundwater resources that can sustain them through to the end of the growing season.
The optimal aggregation timescale for SPEI also varied considerably across the study region. Based on the range of annual water balance values for the study area, aggregation timescales for upland vegetation may be expected to range from four to nine months [13], with vegetation in more arid regions responding to shorter aggregation time scales. We did not observe this pattern in this riparian study; rather, the optimal SPEI aggregation period was relatively consistent across the aridity gradient and was greater than seven months for 3 4 of the ecoregions, which covered most of the state. Thus, in more arid ecoregions in particular, the aggregation timescale for riparian areas is longer relative to that reported for upland vegetation with similar aridity [13], likely reflecting the greater water availability in riparian systems, although we acknowledge that the differences between our study and [13] described above may also influence this result. In our study, variations in aggregation time scales and end months likely reflect slight differences in climatic regimes and snowfall influences among different ecoregions, with longer durations in the northeastern part of the state reflecting the influence of early summer precipitation common to this region [49], for example. In contrast, the northwestern and central parts of the state, where the winter-dominated precipitation regime is most consistent, had start months near the beginning of the water year and end months in May or June (Table 1; Figure 3c,d).
Overall, understanding the strength, sensitivities, and timescales of drought-vegetation relationships in riparian areas, and how they are distinct from those in surrounding upland regions, provides important information for management, which is often strategically focused on riparian areas due to the ecosystem services they provide. These relations can further help to distinguish those riparian areas that are likely to be more immediately vulnerable to changes in the annual water balance from those that are more strongly groundwater dependent and may serve as important hydrologic refugia [54] during extended drought periods or under longer-term changes in climate, provided that these groundwater resources are not depleted due to human extraction of water resources. The optimal aggregation timescale and end months represent the integrated influences of precipitation, evaporative demand, surface water, and groundwater availability, all of which vary on different timescales, on vegetation. Thus, the SPEI metric serves as a useful tool that encapsulates these influences, which are otherwise difficult to evaluate in isolation of each other.

Drought-Adjusted Trend Results
The results presented here demonstrate the utility of moderate resolution Landsat data for increasing understanding of the sensitivity and changes of riparian vegetation in response to climate and management. While several studies have assessed climate sensitivities and trends in vegetation using remotely sensed information [48,55,56] at regional scales, many of these are based on coarser resolution datasets and, as such, do not permit the detection of localized disturbance nor do they capture riparian features, which are often narrow and of limited areal extent. The value of Landsat's 30 m spatial resolution, multi-decadal dataset for characterizing change in riparian systems has long been recognized [19,57], but the potential to apply these data across large spatial extents has only recently been demonstrated [24,54,[58][59][60] and is largely attributable to the recent advances in cloud computing facilitated by the Google Earth Engine [17]. The results from this study can help to distinguish areas of disturbance and longer-term changes in hydrologic regimes from those primarily affected by interannual variations in water availability. Here, we highlight and discuss specific examples of common types of changes that were detected in our study. We follow this with discussion of larger scale patterns of change that were observed.

Defoliation by the Biological Control Agent, Diorhabda Carinulata
The tamarisk leaf beetle, Diorhabda carinulata. was introduced in the western U.S. as a biological control agent for tamarisk, a phreatophytic non-native shrub that has successfully invaded many river corridors and areas of shallow groundwater across the western U.S. The leaf beetle was introduced to the Walker (Figure 5h) and Humboldt Rivers in 2001, with the objectives of reducing evapotranspiration and increasing streamflows [61]. It was successfully introduced to areas surrounding the lower Virgin River (Figure 5e) in 2010 [62]. The defoliation impacts of the leaf beetle are evident at all of these locations, with broad swaths of dense riparian vegetation in these areas showing consistently negative NDVI trends in the Virgin River, reflecting the more recent activity here, and a combination of positive and negative trends in the Walker, reflecting their longer term periodic presence and associated defoliation and recovery over time. The beetle has the potential for rapid dispersal [63] and its presence is not well-documented in Nevada. Given the strong signals observed in our dataset at locations of known beetle occurrences, maps and data provided by this study have the potential to indicate where new invasions may be occurring.

Declines in Surface Water Extents of Perennial Lakes and Wetlands
Increased anthropogenic and atmospheric water demands are resulting in surface water declines in endorheic basins worldwide, with many global studies pointing to locations in Nevada as exemplars of this pattern [64][65][66][67]. This is evident in locations such as Stillwater National Wildlife Refuge (Figure 5g) and Lake Mead (Figure 5e), among others observed in our dataset (Figure 6d). Because water is strongly absorbing in the near infrared range, NDVI values for surface water tend to be very low, thus declining water levels resulting in exposed soils and subsequent vegetation encroachment are expressed as positive trends in NDVI. These cases of positive NDVI trends in response to declining water availability highlight the fact that the interpretation of positive vs. negative NDVI responses is not always straightforward and highlights the need to incorporate additional sources of information to interpret causal factors relating to NDVI trends.

Changes to Grazing Management
Historically, livestock grazing was considered to be one of the greatest threats to riparian areas in the western U.S. and was the most researched threat to these systems in Nevada [68]. The best documented changes to grazing management and the associated impacts to riparian vegetation are for Elko County in the northeastern portion of the state. Up until the early 1990s, most riparian lands were grazed by livestock year-round, with little opportunity for rest and recovery of vegetation [69]. With widespread shifts from year-round grazing to rest-rotation regimes, remarkable recovery of riparian areas [69][70][71] is readily observed as widespread positive NDVI trends in the Upper Humboldt River watershed (Figure 5c), which have been further bolstered by beaver reintroduction and riparian restoration activities across this region [7,69,72,73].

Mine Water Management
Water withdrawals for mining account for 207 million gallons of water per day in Nevada, almost all of which are from groundwater [6]. The majority of these withdrawals are for the purpose of dewatering active mine pits and are redistributed, rather than directly consumed, as exhibited in Figure 5b, where multiple effects of water redistribution have resulted in rising water tables and subsequent vegetation growth in the Boulder Flat area of the Humboldt River basin. Here, water transferred from a dewatered mine into a nearby reservoir has infiltrated through the volcanic geology below. This infiltrated water, coupled with the replacement of localized groundwater pumping for irrigation with dewatered mine water have served to increase spring flow and raise water tables [74], resulting in the greening of areas that were not previously occupied by riparian vegetation. Declining trends are also observed in riparian corridors adjacent to the dewatered mine (Figure 5b), which may have resulted from lowered water tables in this area. Extensive changes in riparian vegetation are observed as positive and negative trends along the length of the upper Humboldt River, where mine dewatering and redistribution have become integrated with other aspects of water management, including irrigated agriculture [74]. Given multiple competing water demands and the requirement to conjunctively manage surface and groundwater in this basin, this region has been a focus for recent hydrologic modeling studies. Indications of vegetation changes provided here may provide a valuable complement to this work and may be used to monitor effects of changing water management in response to mine closures in the future.

Agricultural Land Use
Seventy percent of water use in Nevada is for agriculture [6] that is distributed along many of the alluvial valley floors across the state. Declining NDVI trends in natural riparian vegetation adjacent to irrigated agricultural lands and within lands designated as pasture/haylands are common throughout the state, with the Quinn (Figure 5a) and White (Figure 5d) Rivers providing two examples. Along the Quinn River, groundwater levels have exhibited substantial declines [75] and water rights in the basin are overallocated by over 50% of the estimated supply [76]. With 1/3 of hydrographic basins in the state over appropriated and limited historical groundwater level records, data provided by this study may help to identify places where declining streamflow and groundwater levels are affecting riparian and other groundwater dependent vegetation.

Broader Patterns of Change
Our results indicate that trends in late summer NDVI were statistically significant (p < 0.05) in one-third of riparian pixels in the state. Of these, positive trends are more widespread relative to negative trends, especially those located in mountainous (drainage/channel) landform types (Figure 6a). The paradox of increasing trends in semi-arid vegetation across the globe amidst warming temperatures and associated increases in evaporative demands is gaining increasing attention as researchers have sought to interpret these trends and understand the underlying causes. Possible explanations include the effects of increased carbon dioxide in the atmosphere, which may function to increase plant water use efficiency and soil moisture [77,78], the increases in growing season length due to climate warming [55], or the changes in community composition due to expansion of woody vegetation [79,80]. The significant changes to grazing management that occurred early in the study period (see references in Section 4.2.3) may also contribute to positive trends in NDVI that are observed in this study.
In contrast to the overall more positive trends observed in riparian areas of mountainous regions, riparian areas located in the flatter terrain of wide valley bottoms or plateaus exhibit a wider range of, and overall more negative trend slopes ( Figure 6a, Table 2). Greater dynamism in NDVI in flatter terrain could simply be associated with a more concentrated human footprint in these areas. Urban, agricultural, and industrial land and water uses may be more likely to have direct or indirect effects in lowland regions as the effects of upstream diversion of surface waters, declines in groundwater levels, land use change, and other types of anthropogenic disturbances aggregate with declining watershed position [81]. Indeed, our results indicate tendencies toward more negative trend slopes on private lands and in association with low-impact agricultural uses such as hay and livestock production. A large proportion of riparian areas in Nevada are classified as this land cover type, and although these do not function as natural riparian vegetation would, they do serve important roles for sustaining wildlife in this region [82], suggesting a need to focus efforts on improving vegetation conditions in these areas. In a similar vein, large ranges of positive and negative trends adjacent to perennial water bodies suggest the need for increased management attention toward addressing the loss of surface water bodies and associated wetlands, which may already be affecting changes in migratory water bird habitat use [83].
While there is strong evidence to suggest that anthropogenic impacts of agriculture may explain the more negative NDVI trends observed in valley landforms, these differences could also be influenced by changes in water availability that are driven by greening trends at higher positions within watersheds. In the absence of increasing water use efficiencies that might be caused by carbon dioxide fertilization [78], greening trends and increased water use [84] due to woody plant encroachment are likely to reduce water availability downstream [85,86].

Comparison of Raw and Drought-Adjusted Trends
In this study, we introduce a novel approach to disentangling the influence of interannual variations in drought status from other drivers of variability and change in riparian areas. Unlike prior analyses that assessed trends in the residuals of the interannual climate-NDVI regression relationship at the scale of individual pixels [87][88][89], our approach leverages information from all riparian pixels within an ecologically homogenous region to establish an expected relationship between drought status and NDVI. In doing so, our approach serves to improve disturbance detection in cases where the drought-NDVI relationship for individual pixels masks other changes at the pixel scale [90].
The importance of deconvolving the interannual drought signal on a per ecoregion basis is illustrated by the spatially adjacent classes of lowland 'Upper Lahontan Basin' (13z) versus upland 'Lahontan Sagebrush Slopes' (13k) in Table 1. The same temporal aggregation of November to June was identified for both classes, and the regression R 2 s were very similar and strong. However, the green-leaf response to drought in lowland riparian vegetation, as indicated by the regression slope, was twice as strong as the upland riparian areas. The range of vegetation responses to drought status, seen in Table 1, suggest that attempts to quantify trends or remove transient drought effects by using a per-pixel approach that is applied uniformly across a large region may provide sub-optimal, or even misleading, results.
The drought adjustment increased the statistical significance of trend tests by eliminating noise caused by interannual climate variability over the course of the entire study period, during which there was little trend in climate. Although the drought adjustment was relatively small in this case, the improved ability to identify statistically significant changes will improve the capacity to automatically detect meaningful changes that algorithms or analysts may then interpret. The utility of applying the drought adjustment is more fully realized during the 1995-2013 demonstration period of trending drought status, where virtually all riparian areas exhibited trends in response to drought conditions. In this case, the drought adjustment very clearly highlighted and distinguished areas of disturbance or change from those that were merely responding to the drought. During such times, when vegetation trends have the potential to be masked by trends in climate, this approach is particularly important for teasing apart human-driven causes of change from those controlled by climate variability.
It is important to note that the Mann-Kendall method and results described here are aimed toward identifying monotonic trends in vegetation, and thus have the potential to miss more complex changes or impacts that occurred during the study period. The timing of disturbances, impact onset, and speed of recovery are highly important to change detection using this method and we acknowledge this as a limitation to the study. Furthermore, the riparian areas used for this study represent the estimated contemporary footprint of these ecosystems and thus, our trends do not incorporate information on riparian systems that have been lost due to expanding urbanization or intensive agricultural land use.

Conclusions
The objectives of this study were to quantify relations between riparian vegetation and drought status and to use this information to isolate the influences of interannual variability of drought conditions from longer term trends in vegetation over a 34-year time period. Comparison of our results to similar studies of vegetation response to drought in arid and semi-arid upland systems [13,49] highlights that riparian areas tend to respond at longer timescales and exhibit a range of explanatory powers (R 2 s) that appear to be more strongly driven by gradients of groundwater and surface water dependence than by the gradients of aridity that control the same relationship in upland systems [13]. Because riparian areas are often a specific management focus in arid and semi-arid regions, understanding these distinctions may be useful for better anticipating responses to drought variability and climate change.
We developed a novel method to assess drought-adjusted trends that helps to deconvolve effects of localized disturbances on natural relations between drought and vegetation and applied this method to moderate resolution Landsat data to identify locations, patterns, and potential causes of riparian change across the state of Nevada. Although riparian changes across the state often represent multiple causal factors that are difficult to tease apart in the absence of focused studies, the observations contained in our results provide the basis for new hypotheses to address the roles of longer term climate as well as land and water uses on riparian change in Nevada. In addition, the mapped data product resulting from this study provides a novel long-term perspective for managers seeking to understand the status, trends, and sensitivities of riparian areas and can provide a valuable resource for prioritizing management actions. While our inferences here are limited to our study area, the methods described here are readily transferrable to other regions.