Climate Driven Vegetative Composition Changes in the Big Pine Creek Watershed Using Spectral Mixture Analysis and Time Series Analysis of Landsat Surface Reflectance Data over a 30 Year Period

Alpine ecosystems are crucial laboratories for the study of how changing climatic variables will impact local species assemblages. The steep elevation gradients in these regions provides for analysis of several ecotones within a small area. The biomes that inhabit these areas are particularly susceptible to changing environmental parameters since many exist at the limits of their ranges [1]. Since alpine ecotones represent bioclimatic transitions, species compositional change is high and susceptible to slight alteration in bioclimatic regimes [2]. While many studies have identified biotic response to climate change over large regions, the response at the local and individual ecosystem level are necessary to understand population dynamics that underlie range shifts [3,4].


Introduction
Alpine ecosystems are crucial laboratories for the study of how changing climatic variables will impact local species assemblages. The steep elevation gradients in these regions provides for analysis of several ecotones within a small area. The biomes that inhabit these areas are particularly susceptible to changing environmental parameters since many exist at the limits of their ranges [1]. Since alpine ecotones represent bioclimatic transitions, species compositional change is high and susceptible to slight alteration in bioclimatic regimes [2]. While many studies have identified biotic response to climate change over large regions, the response at the local and individual ecosystem level are necessary to understand population dynamics that underlie range shifts [3,4].
Existing research has focused on the response of individual species, often overlooking important biotic and abiotic interactions that drive community assembly. All the life forms within a local community interact with each other and their physical world forming a complex intricate fabric that identifies the characteristic traits of that assemblage. The predicted trend in climate induced range shifts is for increased extinctions at the warm boundaries and species expansions at the cold range limits [3]. However, in alpine regions, the loss of space with elevation will lead upslope migrating species into a summit trap which will drive extinction rates higher [5,6].
Since the response rate to altered environmental conditions varies among each member of the local assemblage, climate change will drive significant alteration of the interactions between the individual components and the overall functioning of the local community. Since alpine vegetation tends to be long lived, [2], changes in the timing and availability of resources can have significant negative impacts on individual species survival rates while at the same time providing opportunities for competition to allow replacement species to prosper [7]. An ongoing study of the Global Observation Research Initiative in Alpine Environments (GLORIA) site in the European Alps demonstrates this process as species richness has shown a 12% increase in only a 10 year period [8].
Spectral characteristics measured with remote sensing instruments such as the Landsat 5 Thematic Mapper (TM) and Landsat 7 Enhanced Thematic Mapper (ETM + ) enable us to analyze ecological properties of vegetation. Vegetation has characteristic spectral responses such as low red reflectance due to chlorophyll absorption and high near infrared (NIR) reflectance due to the reflectance of the internal structures of the canopy [9]. Changes in surface reflectance can thus be correlated with variation in vegetative cover and plant health. Since the constituents of the plants vary over their phenological cycle, it is also possible to identify the various stages of the cycle such as spring flowering and fall senescence. Changes in the timing of these cycles can serve as an indicator of climate change.
Soil also demonstrates unique spectral characteristics depending on properties such as its moisture, organic matter content and texture [10]. Lower soil moisture content, a possible indicator of water stress in vegetation, would cause higher surface reflectance in the midwave infrared (MWIR) region that can be detected using Landsat data [11]. Higher temperatures combined with lower humidity levels will increase evapotranspiration resulting in less soil and vegetation moisture which will place additional burden on ecosystem vegetation. These effects are heightened in regions experiencing historic droughts such as the southwestern United States [12]. Jackson et al. [10] found *Corresponding author: Sawyer PS, School of Environmental and Public Affairs, University of Nevada Las Vegas, USA, Tel: (702) 295-4221, (702) 630-0131; E-mail: sawyerp2@unlv.nevada.edu that plant water stress will decrease the NIR response while increasing the reflectance in the red region of the spectrum. In addition to altering the spectral reflectance properties of the vegetation, plant stress can alter the geometry of the plant through processes such as drooping or wilting, resulting in a higher soil fraction component of the response signal [10]. Todd and Hoffer [13] found that reduced vegetation moisture content tends to increase visible and MWIR reflectance [13].
Vogelmann et al. [14] examined trends in spectral response in a time series study of the San Pedro Parks Wilderness area in New Mexico for the years 1992 through 2006. Higher elevations were shown to be spectrally stable except for areas infested with western spruce budworm. Some of the lower elevation shrub regions had declines in their short-wave infrared (SWIR)/NIR ratios as did patches of conifer trees suffering from high mortality rates [14]. Loss of available moisture significantly impacts forest growth and overall ecosystem health (Williams et al., [15]). Higher temperatures may also promote pest infestation. Williams et al. [15] found that bark beetle populations increased during warmer periods, especially in forests already suffering moisture deficits induced by higher temperatures. This study determined that maximum temperature (T MAX ) is an ideal surrogate for determining vapor pressure deficit induced forest stress [15].
The ecological response to elevated temperatures and CO 2 levels is complex and will be affected by other factors such as water and other nutrient resource availability. In cold alpine regions where water availability is not limiting, higher temperatures are expected to increase the habitable zones for several species, allowing for upslope migration and increased vegetative cover. In alpine regions, higher temperatures combined with increased atmospheric CO 2 levels will increase photosynthesis resulting in increased biomass; provided other essential resources are not limited [16]. Conversely, where water is limited, higher temperatures will increase plant stress resulting in reduced vegetative cover [17].
Remote sensing using multispectral imagers such as the Landsat 5 Thematic Mapper (TM) and the Landsat 7 Enhanced Thematic Mapper (ETM + ) provide a wealth of data that can be used to monitor for changes in the environment. Large scale regional change are clearly evident from the 30 meter resolution imagery these instruments provide. However, at this resolution, important details within each pixel remain hidden. For remote sensing applications, unless the image is over human controlled agricultural plots, most TM or ETM + image pixels will include several components that cannot be discerned from the raw data. For each pixel, the radiance measured by the sensor in each wave band is composed of a mixture of reflectance energies given off by each of the individual components within that pixel.
In order to elicit the sub-pixel information needed to assess vegetative composition change, we need to employ spectral mixture analysis (SMA). The basic theory of SMA is that in any given pixel, a limited number of dominant components contribute the overwhelming majority of the radiance measured by the sensor. These components are called endmembers (EM). The simplest SMA technique is called Linear Spectral Mixture Analysis (LSMA). For LMSA, the fractional coverage of each EM is proportional to its contribution the overall radiance value of the pixel. This can be expressed mathematically using equation (1) Where Ri is the spectral reflectance for band i of a pixel, fk is the fraction of endmember k within the pixel, Rik is the known spectral reflectance of endmember k within that pixel in band i, ei is the error for band i, and n is the number of endmembers in the pixel [18]. LMSA can be used to unmix pixel spectra from both multi spectral and hyperspectral data [19]. Two methods for solving for fk have been used; constrained and unconstrained. In the constrained method, sum of the fractions must equal 1 as shown in equation (2) [18].
In the unconstrained method, f k is not required to sum to 1 which means the solution will not equal the actual percent cover of each EM [18]. The error for each band e i is defined as the root mean square error (RMSE) expressed as equation (3), Where m is the spectral band [20]. Models of each pixel are developed by varying the fractions of each EM. The model which produces the lowest RMSE is considered the best-fit and those EM fractions are recorded as the solution [20]. The RMSE represents the difference between the measured and modelled value. The RMSE should be in the range of the noise level of the absolute reflectance. For Landsat TM data, this is 0.56%. The analyst can select an arbitrary threshold for the RMSE value for determining when a given pixel has been successfully unmixed; typically 2% [21].
LSMA makes several important assumptions including; each EM is spectrally unique, each EM is constant over the entire spatial extent of the analysis, and the radiance measured by the sensor is a linear combination of all the EM contributions [22]. Typical EM selection includes an EM for soil, one for photosynthetic vegetation, and one for litter or dead vegetation. There are several methods for determining EM spectra including the use of existing spectral libraries, identification of a pixel within the image that is 100% fractional coverage for a particular EM, or obtaining actual field measurements of the EM at the sample site. One of the primary sources of inaccuracy in LMSA is variation of a particular EM spectra within the entire scene (Song, 2005). This limitation on LMSA led to the development of Multiple End Member Spectral Mixture Analysis (MESMA). MESMA allows for the EM spectra to vary pixel by pixel.
MESMA differs from LSMA in that while an individual pixel may contain a limited number of EM, (typically 2 or 3), the entire image may contain many different EM. This allows each pixel to be uniquely analyzed providing much greater definition for the entire image [23]. Since obtaining actual field spectra is often impossible, the use of within-image pure EM pixels is commonly used. These pure EM pixels can be selected using the Pixel Purity Index (PPI) technique. This method is available as a tool in the Exelis Visual Information Solutions (ENVI) software package called the n-Dimensional Visualizer. This tool calculates a score for each pixel based on repeated projections of the pixel data onto a randomly oriented vector which intersects the mean of the data cloud. The result is identification of pixels containing the highest purity of a specific EM [20].
Pre-processing of remote sensing imagery is often performed to reduce data set dimensionality. One of the most common is Principal Components Analysis (PCA). Since the bands in multi-spectral data are often correlated, PCA transforms correlated data into a reduced number of easier to interpret uncorrelated data that are called principal components [24]. When applying PCA to two images, the first few principle components will remain unchanged while the later components will contain the change information. A threshold value is assigned to identify which pixels have changed and which have not [25]. PCA can also be used as a noise reduction technique since the latter principle components consist mainly of noise, they can be zeroed out. PCA can also be used in change detection. When two successive images are processed using PCA, the discontinuities appear as latter components.
A modification of the PCA technique has been developed called the Minimum Noise Fraction (MNF) transformation. Whereas standard PCA applies successive linear transformations in order to maximize variance, the MNF transform was designed to maximize the signal to noise ratio (SNR) by applying linear transformations that successively minimize the noise fraction. The MNF transformed bands with the least SNR are then used for imagery analysis. The MNF transform is equivalent to PCA when the noise variance is the same in all bands [26].
One of the underlying assumptions of PCA is that the SNR declines as the principle component number increases. However, this is not always the case which means that potentially valuable information contained in the latter principal components is often discarded. The MNF transform is performed by first de-correlating and rescaling the noise using the data from the noise covariance matrix then applying a standard PCA analysis to the noise-whitened data. The noise is separated by only using coherent portions of the data [27]. The primary objective of the MNF transform is to identify a small number of noisefree components which can help in the selection of endmembers. Once an MNF transform has been applied, a PPI index is calculated for each pixel to identify the best EM pixels [23].
In previous studies, we examined the spectral response at numerous sample sites in the Big Pine Creek watershed to determine how those sites have changed over the last three decades. In our first study we examined the average spectral response across the watershed and found that both the visible and NIR responses were declining [28]. Vegetation indices are useful tools to analyze for the increase or decline in vegetative surface cover. However, these indices are based primarily on the ratio between the visible and NIR bands. A decline in the visible band surface reflectance from increased visible light absorption combined with an increase in NIR reflectance from higher surface complexity are indicative of increased vegetative surface cover. Likewise, increased surface reflectance in the visible range combined with a decline in the NIR reflectance from less surface complexity is an indicator of a decline in surface vegetative cover. While the simple ratios are useful in identifying increases or declines in surface characteristics, these indices do not provide clear information on what is taking place when the visible and NIR surface reflectance change in the same direction. In order to determine what is occurring at sites where both the visible and NIR spectral responses are trending in the same direction, we need to decompose the spectral responses of those sites using spectral mixture analysis.
In this study we explore ecosystem response to recent climate change by performing a spectral mixture analysis of 16 sample sites which demonstrated similar directional trends in their visible and near-IR spectral response and analyzing trends in sample site composition using time series analysis of Landsat surface reflectance data. We apply a statistical approach to determine trends in the data that are indicative of changing vegetative composition. We present this information by first describing the study area and the data used in the analysis, we then discuss the research approach and methods used to collect and process the data, followed by our results and conclusions. We hypothesize that sites where the visible and NIR spectral responses are changing in the same direction will demonstrate compositional changes that account for the spectral response trends.

Study Area and Data
This section describes the study area and the data used in the analysis. Figure 1 below shows the Big Pine Creek watershed located in California's Eastern Sierra Mountains. Big Pine Creek is a major tributary to the Owens River which is a significant source of fresh water for Los Angeles. The Owens River valley straddles the Great Basin and Mojave deserts with vegetation consisting primarily of pine forests at higher elevations and xeric species at lower elevations. Areas bordering streams and the Owens River are primarily grass dominated meadows [29]. Elevation within the watershed increases from East to West with the higher regions dominated by barren rock and woodlands with the lower regions dominated by mixed desert shrubs.

Study area description
The Big Pine Creek watershed ecosystem owes its existence to snow melt and melt-water from the Palisade Glacier. In addition to being the southern-most glacier in the United States, it is also the largest glacier in the Sierras with a surface area of 1.3 km 2 . It was formed about 3,200 years ago, reaching a maximum extent as recently as 170 years ago [30]. It has been generally in retreat ever since. The Big Pine Creek watershed drainage area covers approximately 82 km 2 and its average flow is 1.8 m 3 /s. Measurements taken in the 1980's indicate that the creek is a gaining stream at the lower elevations in contrast to most other Owens River tributaries which are losing streams [31]. Since all of the living species within this watershed depend on the glacier and snow melt for their survival, the impact of temperature and precipitation variations on the biodiversity of the Big Pine Creek watershed is the focus of this study.
In previous spectral studies of the Big Pine Creek watershed, we examined 105 sample sites. Three sites for each of the top ten predominant land cover classifications present in the watershed and three sites at 100 meter elevation gradients from 1200 meters above sea level to 3600 meters above sea level. At each elevation, a densely vegetated site, a moderately vegetated site and a sparsely vegetated site were selected. While many of the sites demonstrated clear trends in their spectral reflectance consistent with declining or increasing vegetative surface cover, numerous locations exhibited same direction trends in their visible and NIR reflectance bands making interpretation of what is taking place difficult. In this study 16 of those sample sites which were accessible for in situ sampling were chosen for detailed spectral mixture analysis to elicit sub-pixel information which could provide us with evidence of species compositional change not discernable from the multi-spectral 30 meter resolution Landsat imagery.

Data
The data in this study includes Landsat surface reflectance data obtained from the USGS Earth Explorer web site, ground truth spectra from 116 surface cover samples collected in situ throughout the study area during July 2014, and modeled endmember spectra and abundance values derived from the Landsat data using the ENVI 5.1 software package.
Surface reflectance data: The Landsat program has been providing earth observation remote sensing data to the scientific community for four decades. The first Landsat satellite was placed in orbit in 1972 with Landsat 7 remaining operational today. Landsat 5 was only recently taken off-line. The latest generation satellite, Landsat 8, was launched on February 11 th , 2013 and is now operational. Data for this study includes imagery from both the Landsat 5 TM and Landsat 7 ETM + sensors. Table 1 details the six reflective bands of the Landsat sensors, covering the visible (blue, green, and red), near infrared (NIR), shortwave infrared (SWIR), and mid-wave infrared (MWIR) regions of the spectrum as well as their ecological applications. Bands 5 and 7 are sometimes referred to as SWIR1 and SWIR2. However, in this study, we refer to band 5 as the SWIR and band 7 as the MWIR. Band 6 covers the thermal infrared (TIR) and the data from this region is not used in this study. These descriptions are retrieved from the Northern Arizona University Infrared Spectrometry Laboratory. (http://www.cefns.nau. edu/seses/llecb/Spectrometer/RemoteSensing.html).
The Landsat imagery used in this analysis was acquired for 30 dates in the month of July from 1984 through 2013. Most of the imagery used in this analysis is from Path 42, Row 34 with four of the images from Path 41, Row 34. Both image ID ground swaths cover the entire study area. The imagery acquisition date and time for the data used in this analysis are listed in Table 2 below. This surface reflectance data product for each of these imagers was obtained from the EarthExplorer web site operated by the United States Geological Survey (http://earthexplorer.usgs.gov/). Since the period of maximum leaf area index generally occurs in the mid-June to mid-August time frame [32], only imagery in the July time frame was considered for this analysis in order to minimize the impacts of the phenological cycle on the reflectance data.
Meteorological data: Meteorological data referenced in this study was obtained from the University of Oregon's Parameter-elevation Regressions on Independent Slopes Model (PRISM) web (http://www. prism.oregonstate.edu/). According to its website, PRISM data are modeled estimates based on point data and a digital elevation model and is available for the entire continental US at 4 km resolution. All sample sites in this study fall within 11 PRISM grid cells as shown in Figure 2.

Research Approach and Methods
This section contains a description of the research approach and the methodology used to collect and process the data.

Research approach
This study examines how the surface cover in the watershed has varied over the last 30 years at 16 sample sites in which the spectral reflectance trends in both the visible and NIR bands change in the same direction. This is accomplished by performing a spectral mixture analysis of each of the sample sites for each year in the study, then performing a time series trend analysis of the endmembers identified.

Research methods
The research methodology consists of surface reflectance data collection; endmember determination; data processing; and statistical analysis. Each step is described below.
Data collection: In order to perform a temporal study comparing the physiological changes over time at each of the sample sites, surface reflectance values for each year of the study period were obtained from the USGS Climate Data Record (CDR) archive. The climate data referenced in this study were obtained from the PRISM web site by entering the geographic coordinates of each sample site into the data base and downloading the PRISM data set for each site over the 30 years of the study period.
A total of 116 endmember samples were collected in the field. These samples fell into nine broad categories including six types of photosynthetic vegetation, non-photosynthetic vegetation such as litter, and two barren surface cover types, soil and rock. Due to the remoteness of the site and its rugged terrain, taking field spectral measurements was not practical. Therefore, the samples were sealed in plastic bags and stored on ice for transport back to the lab where their spectra was measured with an Analytical Spectral Devices, Inc. (ASD) 0.35 to 2.5 µm Flexscanspectroradiometer. The 1 nm bandwidth spectra generated by the ASD instrument was rescaled to match the Landsat spectral bands using the spectral resampling application in the ENVI 5.1 software.
Surface reflectance data: USGS surface reflectance data is generated from a software package known as the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS). The surface reflectance data is computed by applying an atmospheric correction to the raw Landsat imagery [33]. This atmospheric correction uses the Second Simulation of a Satellite Signal in the Solar Spectrum (6S) radiative transfer model to account for various atmospheric column constituents including water vapor, ozone, and aerosol optical thickness [34].
The LEDAPS process uses average daily lamp brightness history to obtain calibration coefficients based on acquisition date. These   [34]. The LEDAPS process converts at-sensor radiance to top-of-atmosphere (TOA) by an algorithm that incorporates solar irradiance derived from the MODTRAN model, bandpass, earth sun distance and solar zenith angle (Masek [34].The LEDAPS atmospheric correction assumes particle scattering and gaseous absorption can be decoupled [34]. LEADAPS applies Moderate Resolution Imaging Spectroradiometer (MODIS) atmospheric correction routines to Landsat data that correlates surface reflectance is with TOA reflectance using (4), where ρ s is the surface reflectance, T g is the gaseous transmission, T R+A is the Rayleigh and aerosol transmission, ρ R+A is the Rayleigh and aerosol atmospheric intrinsic reflectance, and S R+A is the Rayleigh and aerosol spherical albedo [34].
The 6S radiative transfer model is used to derive surface reflectance from (6) with the input of aerosol optical thickness (AOT), atmospheric pressure and water vapor. The Total Ozone Mapping Spectrometer (TOMS) carried by the Nimbus 7, Meteor 3 and Earth Probe satellites provides the ozone concentration data. For the 1994 through 1996 time period when TOMS data was unavailable, vertical sounder data from the National Oceanic and Atmospheric Administration (NOAA) was used. Rayleigh scattering is adjusted to local conditions using surface pressure data from NOAA's National Center for Environmental Protection (NCEP) [34].
This atmospheric correction methodology uses a dark dense vegetation procedure developed by Kaufman et al. [35] to determine AOT from the imagery. This technique is based on the assumption of a linear relationship between surface reflectance in the visible bands and the surface reflectance in the short wave band (2.2 µm where surface reflectance is not affected by the atmosphere) based on the physical correlation between bound water absorption and chlorophyll absorption. Using this procedure to calculate surface reflectance in the visible bands then allows for the determination of AOT by comparing the TOA reflectance to the surface reflectance. Since this technique only determines the AOT in the blue region, a continental aerosol model is used to determine AOT in the other spectral regions. The AOT, atmospheric pressure, ozone and water vapor data are then processed by the 6S model to convert TOA reflectance to surface reflectance [34]. The USGS CDR data set provides us with observed surface reflectance values for each of the six reflectance bands for all sample sites in each year of the study.
Spectral endmember determination: Spectral mixture analysis for this study was performed using the spectral hourglass wizard toolkit in the ENVI 5.1 software application. The spectral hourglass wizard provides a systematic process for determining spectral endmembers within a given region of interest. The first step in the process is application of a minimum noise (MNF) transform to the reflectance data. This transform determines the inherent dimensionality of the   data and generates an eigenvalue plot for each MNF band. In order to improve spectral processing results, MNF bands containing only noise are discarded to preserve linear independence among the components. The next step in the process is to derive endmembers from the image. The spectral hourglass wizard calculates a pixel purity index (PPI) by repeatedly projecting n-dimensional scatter plots on a random unit vector. A PPI image is created with the value of each pixel corresponding to the number of times the pixel was determined to be extreme. This process identifies the purest pixels in a scene. Spectral endmembers are then determined using the n-dimensional visualizer which locates and identifies the most extreme spectral responses in a data set. At this point, the spectral hourglass wizard produces a spectral unmixing abundance image which assigns abundance values to each spectral endmember in each pixel.
The ENVI Spectral Hourglass Wizard generated a set of endmember spectra for each sample site. In order to determine what each of the endmember spectra represented, the spectra was compared to the ground truth spectra. The abundance of each of the ENVI generated spectral endmembers was determined for each date in the study period and a trend analysis was performed to determine how the fractional coverage of each of the endmembers has changed over the last 30 years.

Statistical trend analysis:
The non-parametric Mann-Kendall (MK) trend test is used to establish the presence of trends in the spectral endmembers over the last 30 years. This analysis essentially determines if a set of values (y) are increasing or decreasing over time. Mann-Kendall analysis looks at the sums of the signs of the differences between successive data points and calculates a score or "S" statistic with the following properties: for S<0 (values are decreasing over time); for S>0 (values are increasing over time). The magnitude of the S-statistic is a measure of the strength of the trend. For a sample size of 30, S values of ±111 indicate a statistically significant trend with a p value of <0.05. This means the null hypothesis of no-trend in the data can be discarded with the risk of committing a Type II (rejection of a true null or H 0 ) error at less than 5%. The MK S-statistic is calculated using (5) wheren is the number of observations and y i (i=1…n) is the value at time T i and y j (i=1,…, n) is the value at time T j [36]. Variance in the S statistic is calculated as ( ) n(n 1)(2n 5) 18 Var S − + = (6) This variance assumes there are no tied pairs in the data. If tied pairs are identified, the software program applies a continuity equation which assumes a normal distribution for S with a zero mean. The variance is used to determine the probability (p) of obtaining a value of S greater than that calculated for the given number of data points when no trend is present. The probability statistic is determined from the Z score which is defined as: Where n is the number of observations. Kendall's tau is similar to the correlation coefficient in linear regression. The magnitude of the trend is determined using the Sen's slope estimation with confidence intervals defined as the upper and lower estimate for the mean value of the slope. Sen's slope is determined by calculating the slope at each data point and taking the median of those slopes as the magnitude of the trend as shown; These calculations are carried out in Excel using the XLSTAT addin statistical application. This program generates the S statistic as well as the probability (p) value which is used to quantify the statistical significance of the trend. The confidence factor (risk of rejecting a true null) is defined as (1-p) * 100%.

Results and Discussion
Ground truth endmember spectra A total of 116 surface cover samples located throughout the study area were collected in situ in July 2014 including three to four samples from each of the 16 sample sites used in this study. The surface cover types are classified as either photosynthetic vegetation (PV), nonphotosynthetic vegetation (NPV), or barren surface. There were six broad classes of PV present: broad leaf trees including aspen, birch, and cottonwood; narrow leaf trees including willow and woodrush; needle leaf trees including pine and fir; sage bush; stem shrubs including juniper, hardhack and sedges; and leafy shrubs including manzanita, monkey flower and desert peach (Figures 3 and 4).
The NPV consisted primarily of litter and dead vegetation. Barren surface included various soil and rock types. All of the ground truth samples were categorized into those nine surface types and a composite spectra of each surface cover was created. Figure 5 shows the composite ground truth spectra of all the PV types collected in the field. Figure  6 shows the composite ground truth spectra of non-photosynthetic (NPV) along with composite soil and rock spectra. Reflectance values have a scaling factor of 10,000 to match the Landsat surface reflectance data set.
For each sample site, the endmember spectra produced by the ENVI software application were compared against the ground truth endmember spectral library to determine their classification. The abundance values for each endmember type, (PV, NPV, and Soil/Rock), were compiled for each sample site for each year and a trend analysis was performed to determine how the surface cover has changed over the last 30 years. Figures 7 and 8 provide an example of this procedure. In Figure 7, sample site #1 ENVI generated spectral endmembers for photosynthetic vegetation from the Landsat 5 TM imagery for the year 2000 are compared against the spectra produced from vegetation samples collected at that site. At this particular location, vegetation is primarily California sage (Artemisia californica) and Big sagebrush (Artemisia tridentata). Although the raw reflectance values are significantly different, the shapes of the spectra allow us to clearly identify that the ENVI generated spectra nD Class #3 matches the ground truth spectra for Big sagebrush while the nD Class #4 matches the ground truth spectra for California sage.
ENVI Generated Spectra for Sample Site #1 PV Ground Truth Spectra for Sample Site #1 PV In Figure 8, sample site #13 ENVI generated spectral endmembers for non-photosynthetic vegetation (litter) and soil from the Landsat 5 TM imagery for the year 1994 are compared against composite litter and soil ground truth samples. As with the vegetation spectra, the soil and litter spectra differs in their data values, but their distinctive shapes are clearly distinguishable.
In general, the ENVI generated spectra data values were much lower than those produced by the ASD instrument. The spectra measured in the laboratory are contact samples, meaning there is no atmospheric column between the sample and the detector. This contrasts with the full atmospheric column between the study site and the Landsat sensor orbiting 700 km above the earth.

ENVI Generated Spectra for Sample Site #13 Soil and Litter Ground Truth Spectra for Composite Ground Truth Soil and Litter
In addition to the spectra, the ENVI software generated abundance images of each modeled endmember in every pixel. The constrained option was selected so that the abundance values reflected fractional coverage of each endmember within each pixel. The software also produced an RMS image indicating the root mean square error associated with the derived endmembers for each individual pixel.
Average error values for all 30 years of the study are shown in table four for each site. Table 3 shows the trends in surface cover for those 16 sites where the trends in spectral response would be consistent with a change in surface composition, (visible and NIR trending in the same direction). In most instances we see trends in the PV and NPV going in opposite directions. For example, 7 of the sites show declines in PV along with increases in NPV while at 4 sites, the PV is increasing while the NPV   At sites where the PV and NPV are moving in the same direction, we see barren surface going in the opposite direction as one would expect. Where the PV and NPV are receding, we see increased barren surface contribution to the observed surface reflectance signal. Likewise where PV and NPV are increasing, the barren surface contribution is in decline.
An interesting observation from this data is that the trends in barren surface do not necessarily correlate to the Tasseled Cap transformation for brightness (TC B ). Table 4 shows that for 6 of the 16 sites, the TC B value is trending in the opposite direction from the barren surface cover. Similarly, the Tasseled Cap greenness (TC G ) transformation would be expected to trend in the same direction as the PV fractional coverage. However in 9 of the 16 sites we find the opposite. The Tasseled Cap wetness (TC W ) index is often correlated to increases in surface vegetation cover as higher vegetation cover means more water content in the scene. However, in 7 of the 16 sites we see the PV trends going in the opposite direction from the TC W trends.
These results highlight the caution one must use when interpreting the meaning of the Tasseled Cap transformations. This is especially important when attempting to interpret what is taking place when the visible and NIR spectral responses are trending in the same direction. Although most of the instances cited do not involve statistically significant trends, at sample site #13, we see statistically significant positive trends in PV while at the same time we see statistically significant declines in TC G and statistically significant increases in TC B . For a detailed explanation of how we derived the vegetation indices and Tasseled Cap transformations, see Sawyer and Stephen [28].
The real advantage that spectral mixture analysis has over simple spectral reflectance derived vegetation indices is the ability to discern changes in the endmembers within an individual pixel. This information is essential in determining how climate change is impacting local species assemblages. Table 5 shows the average abundance of the six predominant PV types over the last 30 years at each of the 16 sites in this study. This table shows that our study site vegetation consists primarily of shrubs with some deciduous and conifer trees. Most of the trees are located within a hundred meters of the Big Pine Creek while the shrubs are ubiquitous throughout the study area.
Tree species identified in the study area include leaf species such as Willow, Woodrush, Cottonwood, Birch, Poplar and aspen, along with conifers including Jeffery Pine and Red Fir. Two sages were identified; Big sagebrush and California sagebrush. Numerous shrub species were identified. The shrubs were classified as stem type for those with needle like leaves or sedges and as leaf type for those with broad leaves. Stem type shrub species include Juniper, Golden Hardhack, Green Ephedra, Black Greasewood and sedges. Leaf type shrubs are primarily Manzanita, Monkey flower and Desert Peach. Table 5 shows the average fractional surface cover of each of the PV types at each site for the 30 years of the study period as derived from the spectral mixture analysis of the Landsat imagery. Table 6 shows how ENVI Generated Spectra for Sample Site #1 PV Ground Truth Spectra for Sample Site #1 PV Figure 7: Comparison of ENVI generated spectra for sample site #1 with ground truth spectra of the vegetation collected at that site.

ENVI Generated Spectra for Sample Site #13 Soil and Litter
Ground Truth Spectra for Composite Ground Truth Soil and Litter the individual endmembers have trended at each of the sample sites where the spectral analysis is consistent with compositional change. Table 7 provides a description of what the trend values are showing. In this summary, trends are defined as large if their |S| value exceeds 60, moderate for 20<|S|<60, small for 5<|S|<20, and slight for |S|<5.
The results of our spectral mixture analysis of the 16 sites that we hypothesized were undergoing compositional change are consistent with that theory. We see vegetative compositional changes at each of the sites examined, with three of those sites consistent with large compositional changes. Looking at the combined PV trends, for the two sites with statistically significant trends, 13 (S=148) and 14 (S=-105), the results are not always consistent with the vegetative indices and Tasseled Cap transformations. At site 13, the spectral mixture analysis indicates vegetation is increasing while all three vegetative indices are declining. Likewise, the TC B , TC G , and TC W scores all suggest declines in vegetative surface cover. At site 14, the spectral mixture analysis indicates declining vegetative surface cover which is in agreement with two of the vegetative indices (SAVI and MSAVI 2 ) and the tasseled cap transformations. However, the NDVI trend at this site is consistent with a small increase in vegetative surface cover. These results demonstrate the difficulty with assessing vegetative compositional change using only vegetative indices and Tasseled Cap transformations.
Looking more closely at site 13, we see large increases in sage with a moderate increase in stem shrubs and a small decline in leaf shrubs. If we examine how their individual contribution to the composite spectral response has changed, we see that since sage and stem shrubs have higher NIR signatures, age (ρ NIR =0.7153) at 21%, stem (ρ NIR =0.6554) at 40%, and leaf shrubs (ρ NIR =0.5979) at 13%. Higher sage and stem will result in a higher composite NIR spectral response which is consistent with what the trend data shows for this sites NIR response. In the Red band, an increase in sage (ρ NIR =0.2685) at 21%, and stem (ρ NIR =0.1769) at 40%, combined with a decline in leaf shrubs (ρ NIR =0.1667) at 13% will likewise result in a higher cumulative Red band response since although higher vegetative cover is expected to result in increased Red band absorption from higher chlorophyll levels, the composition change to species with higher Red band reflectance result in a higher

Meteorological data
Ecological changes we have identified in the Big Pine Creek watershed are consistent with warming temperatures. In our previous paper, we identified statistically significant increases across the study area in both the maximum temperature (T MAX ) and minimum temperature (T MIN ). Looking at monthly trends, we found that for the maximum temperatures, the largest increases are taking place in the summer with smaller increases in the winter. Higher temperatures are an important factor in driving ecological changes since all biological processes are at their essence chemical reactions and increased temperatures will increase reaction rates. Increased biological activity can alter vegetative composition by changing the availability essential nutrients. Some nutrients will be more available through faster litter breakdown while some nutrients will be consumed at faster rates. This change in resource availability will drive changes in species composition to those species that are better suited to the new environmental conditions and resource make-up [28].
Precipitation and Big Pine Creek stream flow trends were also examined to determine if the moisture deficit conditions were impacting ecological responses over the study period. Here we found that although there was a slight decline in precipitation in the watershed, stream flow was slightly increasing which is consistent with warmer temperatures increasing the melt water contribution to the stream flow from the Palisade glacier [28]. In addition to driving biological activity, higher temperatures will also increase evapotranspiration. This will reduce moisture availability and stress vegetative species, especially those that are not drought tolerant. For the monthly minimum temperature trends, we see the largest increases are taking place in the summer and fall with smaller increases in the winter and spring. What these data demonstrate is that the summers are getting warmer and the winters are getting milder. This is an important finding since as discussed

14
Moderate declines in sage (Big Sagebrush) and stem shrub (Juniper) and a large decline in leaf shrub (Manzanita) 15 Moderate increase in sage (Big Sagebrush) with small decline in stem shrub (Juniper) and a large decline in leaf shrub (Manzanita) 16 Slight increase in narrow leaf tree (Willow) with a large increase in needle (Red Fir) and small declines in stem shrub (Juniper) earlier, warmer summers will increase evapotranspiration during the dry season, increasing potential water stress in the vegetation. Milder winters will also result in reduced water storage capacity as less precipitation will fall as snow, which also results in reduced water supplies in the warmest time of the year [28].
The seasonal trends found in this analysis closely align with future climate regimes predicted by general circulation models showing milder wetter winters and hotter drier summers [37]. Lenihan et al. [37] show that these future climate scenarios can produce shifts in the vegetative composition. In particular, their biological distribution model simulations suggest a shift from shrubs to grasslands under these conditions [37]. The temperature trends demonstrate that the Big Pine Creek watershed is at heightened risk from climate change and highlight the need to develop strategies to adapt to the new climate paradigm.

Confidence levels
Multitemporal satellite imagery is impacted by several factors including changes in sensor response, sensor stability, atmospheric effects, and illumination effects [38]. Geometric pixel registration errors are generally less than ½ pixel [39]. Radiometric uncertainty for the TM data are approximately 5% [40]. The USGS surface reflectance data set has been assessed against MODIS surface reflectance data and found to be highly correlated with discrepancies between 2.2 to 3.5 percent [41].
The spectral mixture analysis performed for this study generated RMS errors of less than 1% on average. The average RMSE error values for each site are shown in Table 5. However, the interpretation of what each ENVI generated endmember represents is somewhat subjective. Accuracy of this interpretation is dependent on the ability of the analyst to correctly match the modeled spectra to actual ground truth spectra.

Summary and Conclusions
This study examined the changes in the ecosystem of the Big Pine Creek watershed as measured by trends in sample site composition at 16 locations over a 30 year time span from 1984 through 2013. Analysis of trends in surface reflectance and vegetation indices do not provide reliable information when the visible and NIR bands are trending in the same direction. We hypothesized that compositional change was taking place at these sites which can account for the same directional trends in the visible and NIR regions. Spectral unmixing can help us elicit information not discernable from simple trend analysis of Landsat spectral bands. In this analysis, we see evidence that compositional changes are taking place at each of the sites. Although none of the individual species components demonstrated statistically significant changes, the cumulative change in photosynthetic vegetation at two sites were statistically significant.
In remote areas that do not have historical surface cover composition data, quantitative analysis of how each individual component affected the composite spectral signal is not possible. Analysis of vegetative indices and Tasseled Cap transformations break down when the Red and NIR bands change in the same direction. However, by performing a spectral unmixing of the sample sites, we can look at the qualitative results of trend data to infer the impact compositional changes are having on the composite spectral response in each band. This analysis demonstrates a way to elicit a plausible explanation for the spectral responses recorded by the Landsat imager. We observed at site 13 how a compositional change with an overall increase in vegetative surface cover can result in both higher Red and NIR signatures even though the vegetative indices suggest decreasing vegetative surface cover. This suggests that previous studies that examined trends in vegetative indices where the visible and NIR trended in the same direction cannot be relied upon to definitively state that increased or decreased vegetation index trends mean higher or lower levels of photosynthetic vegetation. What may actually exist are vegetative compositional change or varying amounts of NPV or litter surface cover. Analysis of vegetative change that relies solely on changes in the composite spectral response may miss important compositional changes that may contradict vegetative indices and Tasseled Cap transformations.