Investigating the Effects of Land Use and Land Cover on the Relationship between Moisture and Reflectance Using Landsat Time Series

To better understand the Earth system, it is important to investigate the interactions between precipitation, land use/land cover (LULC), and the land surface, especially vegetation. An improved understanding of these land-atmosphere interactions can aid understanding of the climate system and modeling of time series satellite data. Here, we investigate the effect of precipitation and LULC on the reflectance of the land surface in the northern U.S. Great Plains. We utilize time series satellite data from the 45 year Landsat archive. The length of the Landsat record allows for analysis of multiple periods of drought and wet conditions (reflecting climate, as well as weather), such that the precipitation-reflectance relationship can be investigated robustly for every individual pixel in the study area. The high spatial resolution of Landsat (30 m) allows for investigation of spatial patterns in weather (i.e., precipitation extremes) interactions with land surface reflectance at the scale of individual fields. Weather history is represented by a drought index that describes effective moisture availability, the Standardized Precipitation and Evaporation Index (SPEI). We find that effective moisture has a robust and consistent effect on reflectance over many types of land cover, with ∼90% of all pixels having significantly ( p < 0.01 ) higher visible reflectance during dry periods than during wet, occurring in nearly all regional, temporal, and LULC categories investigated. In grassland, the relationship is especially strong; there is an average reflectance increase of more than a third between very wet and very dry conditions (red band), and ∼99% of pixels have a significant relationship. In cropland, the effective moisture-reflectance relationship is more variable, suggesting that management decisions are an important factor in cropland-reflectance relationships.


Introduction
The relationship between land use/land cover (LULC) and climate or weather is of profound importance in the Earth system. Changes in climate and LULC have the potential to produce substantial disruption for individuals and societies [1,2]. Interactions between LULC and weather are of particular importance for understanding the climate system, as they can amplify the impacts of changes in either, through feedbacks. Understanding these land-atmosphere interactions is critical to predicting the effects of changes in either LULC or climate on both, and on society [3].
An important instance of land-atmosphere interactions is the relationship between precipitation and surface reflectance (at a variety of wavelengths). Precipitation and the albedo of the land surface have the potential to produce feedbacks; for example, during drought, albedo usually increases, altering atmospheric energy transfer and often further inhibiting precipitation [4,5]. To improve understanding of these feedbacks, it is important to have evidence for the magnitude of both the effects of albedo on precipitation and the effects of precipitation on albedo. Estimates of the magnitude of the effect of moisture on surface reflectance at a variety of wavelengths help to improve estimates of precipitation-albedo effects and have the potential to inform modeling studies of topics such as land-atmosphere interactions and analysis of time series remote sensing data.
A number of studies have observed a change in remote sensing observations during drought or low soil moisture, e.g., [5][6][7][8][9][10], and this relationship has been investigated at length to produce drought indices [11,12]. However, these studies have generally covered short time periods (less than the 30 years commonly used in creating climate baselines), implicitly relying on a space-for-time-type approach rather than investigating every pixel location across multiple drought episodes. This approach potentially over-emphasizes the attributes of particular drought events and reduces the potential for investigating relationships across space (including between LULC classes) since pixels do not each have a robust estimate independent of other pixels. Furthermore, studies have generally investigated coarser spatial scales than 30 m Landsat data. In addition, much previous work has been focused on vegetation health rather than reflectance. Changes in vegetation are a critical part of reflectance variation, but soil moisture also affects the reflectance of the soil directly [13], so a study of the effect of weather including dry and wet extremes on a direct measurement of reflectance is valuable.
Since different land cover types (e.g., grassland, barren land, forest, and cropland) are all likely to behave differently during periods of high or low moisture, it is also critical to take LULC into consideration when investigating land-atmosphere interactions [3,14]. Grasslands tend to respond to drought quickly, with reduced evapotranspiration and leaf senescence, while evapotranspiration and gross primary productivity in some forest types can change more slowly after the initiation of drought [15,16]. The drought response of croplands can be highly variable since there is a wide variety of management regimes and irrigation practices [17]. Given the differences in LULC types, changes in LULC have the potential to affect reflectance, albedo, and climate. Land cover has a direct effect on reflectance, while land use may affect land cover. Despite the importance of LULC, there has been little work investigating the LULC dependence on large-scale drought responses.
Remote sensing time series investigations over long (multi-decadal) timescales and large (thousands of km 2 ) areas can help to distinguish underlying land-atmosphere relationships as data are gathered synoptically under a wide range of conditions. Landsat is well suited for this type of analysis, as it has extensive coverage across a long period of record (1972-current) including multiple drought episodes [18,19]. Analysis of multiple drought episodes establishes the range of observed responses to drought, allowing a focus on the average drought response of a location instead of the particularities of a specific drought. The 30 m resolution of Landsat is an advantage over studies using 250-1000 m resolution or coarser, because it enables distinguishing smaller features, including field/management boundaries. Furthermore, there are more unmixed pixels that contain only one LULC type [20]. Early data from Multispectral Scanner System instruments (MSS; first launched 1972) are lower quality and resolution (60 m) than data from more recent missions, but are still useful in this kind of study, where a long baseline is potentially more important than high radiometric accuracy [21].
Drought data are necessary to investigate effective moisture-reflectance relationships. The Standardized Precipitation-Evapotranspiration Index (SPEI) is a relatively new drought index and has been utilized by a number of studies, e.g., [22]. The SPEI indicates the water deficit compared to normal conditions for a location [23,24]. One advantage of the SPEI calculation is that it includes the effect of temperature on water availability by subtracting potential evapotranspiration (PET) from precipitation [25]. Furthermore, as the SPEI uses the same statistical approach for all environments and can be applied globally without local correction factors, it provides for the assessment of moisture conditions at any desired timescale and, thus, can describe the water stress on vegetation. Here, we understand the SPEI as a measure of effective moisture.
The goal of this paper is to evaluate the effect of LULC on the relationship between effective moisture and reflectance and to investigate methods for examining geographic patterns in this relationship across the landscape. For each 30 m pixel in the study area, we create a time series of Landsat reflectance (visible, near-IR, and shortwave-IR bands) and compare to a corresponding time series of SPEI data. Summary statistics and linear regression fits of SPEI versus reflectance are averaged by the LULC category to investigate differences among the LULC categories. Individual droughts (and individual wet periods) each have their own unique and interesting features, and so, it is important to sample as many drought events as possible to avoid interpreting patterns seen in one or two droughts as a more general result. To maximize the temporal period in this study, we include MSS data (begin in 1972) in the analysis. The results of this work are intended to contribute to models that might benefit from an improved understanding of the relationship between effective moisture and reflectance, such as time series modeling and land-atmosphere interaction modeling.

Study Area
This study focuses on a study area in the northern U.S. Great Plains (Figure 1). The Great Plains are a suitable location for this study because the area has experienced frequent droughts and the vegetation (which is mainly herbaceous) is responsive to precipitation [26,27]. Furthermore, research on semi-arid grasslands (including the Great Plains) is of global importance as grasslands compose >30% of the Earth's land surface and are especially vulnerable to anthropogenic LULC conversion [28].
The study area has a semi-arid climate, with average annual precipitation ranging from ∼500 mm/y in the southeast to <400 mm/y in the northwest [29]. Precipitation is highest in the spring, with as much as half of annual precipitation falling in the May-July period. Average annual temperatures are ∼6-9 • C, with average January minimums of −15 to −11 • C and average July maximums of 28-32 • C (based on meteorological stations within the study area, 1981-2010; Figure 2 [29]). On average, there is snow depth of an inch or more on 70 days per year [29]. Vegetation in the study area can be limited by either water or energy, depending on weather conditions. As measured by the SPEI, large parts of the study area were in drought in 1975, 1980, 1988, 2002, 2006, and 2012. The droughts of 1980, 1988, and 2012 are noteworthy as being part of the three most expensive droughts that have occurred in the U.S. since 1980 [30]. The 2002 and 2006 droughts were less widespread, but within the study area, their intensity was second only to 1980. Wetter periods occurred in 1978, 1982, 1986, 1995, 1999, 2010-2011, and 2014. Since both wet and dry events were distributed throughout the study period, the Landsat record provided a baseline for investigating relationships that were linked to precipitation rather than to a particular drought event.
The study area is dominated by grassland (∼72%), with cropland (∼17%) and pasture/hay (∼6%) the next most common (National Land Cover Database 2001; [31]). Livestock grazing is common on grassland throughout the study area; the U.S. Department of Agriculture 2012 Census of Agriculture found 585,704 cattle and calves in the eight counties dominating the study area (4,237,738 ha total county area; [32]). Grassland in the study area is mixed-grass prairie, with species including western wheatgrass (Pascopyrum smithii), blue grama (Bouteloua gracilis), needle and thread (Stipa comata), and green needle-grass (Stipa viridula) [33]. C 3 species are most common, but C 4 grasses compose as much as 50% of production in some parts of the study area [34].
Cropland in the study area is often divided into fields that are at least 65 ha (a quarter section), but can be substantially larger. Hay and wheat are the most common crops grown in the study area. Specifically, the crops covering the most area in 2014 were non-alfalfa hay (∼31% of crop area), spring wheat (26%), winter wheat (12%), corn (9%), sunflower (5%), fallow (5%), and alfalfa (4%) (data from the Cropland Data Layer, [35]). In South Dakota, spring wheat is usually planted April 8-May 12 and harvested July 27-August 20 (this schedule takes advantage of early summer rains), winter wheat is planted September 9-October 5 and harvested July 14-August 5, and non-alfalfa hay is harvested (more than once) June 20-August 14 [36]. These dates are dependent on precipitation. Most cropland is rainfed (>99% in 2007; [37]), although early in the study period, there were small areas of irrigation distributed across the region [38]. Croplands are typically located on table lands parallel to streams. LULC change has been moderate in the region during the past ∼40 years, with LULC change in 7.4 ± 2.0% of the Northwestern Great Plains ecoregion in 1973-2000 [39]. Conversions between agriculture and grassland were the most frequent changes, driven by economic forces, technology, and farm policies including the Conservation Reserve Program (CRP). The CRP removed more than 200,000 ha from production between 1985 and 1992 in counties roughly encompassing the study area [40].

Methods
For this study, Landsat data were compared to SPEI and LULC. A Landsat reflectance time series was constructed for each 30 m pixel in the study area for visible, near-IR, and shortwave-IR bands. Each Landsat time series was then compared to a corresponding SPEI time series. Histograms of reflectance by SPEI and LULC were computed to display the general relationship between effective moisture and reflectance, and linear regression of SPEI versus reflectance was used to investigate this relationship in more detail. We included MSS to extend the temporal range to 1972, in order to incorporate as many dry and wet periods as possible in the analysis. Inclusion of MSS led us to utilize top of atmosphere (TOA) reflectance instead of surface reflectance, because atmospherically correcting MSS was more difficult than later instruments due to the lack of shortwave-IR bands [41]. We also investigated the relationship between the Landsat atmospheric correction and the drought index to evaluate the impact of the utilization of TOA in this analysis.

Landsat Data Processing
Scene-based pre-collection Landsat data from 1972-2014 were acquired from the USGS Landsat archive [20,42]. In general, all available observations from the four-scene study region shown in Figure 1 with <50% cloud cover were utilized (1430 observations, listed in Supplementary Materials S1). This included observations from MSS, Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and Operational Land Imager (OLI) instruments (Table 1). For Landsat 1-3, this area approximately covered six scenes in the Worldwide Reference System-1. Since vegetation is dormant and the ground is frequently snow covered during the winter in the study area, only dates between April 15-October 15 were included in the analysis. One priority of this study was the utilization of data back through the MSS period of the Landsat program, in order to include 1970s era droughts in the analysis ( Figure 3). Data from MSS are the earliest of this type of Earth observation, and the data archive for MSS is extensive in the study area [43]. This time period allowed the incorporation of three noteworthy droughts (1975, 1977, and 1980) in addition to four noteworthy droughts in the later record (1988, 2002, 2006, and 2012). Compared to TM, ETM+, and OLI, the MSS data have lower spatial resolution (60 m versus 30 m), fewer spectral bands (4 versus ≥7), a lower frequency of observation (18 day vs. 16 day repeat cycle; <20 observations/y versus ∼60/y in the 2000s for this dataset), reduced radiometric resolution (6 bit versus 8 bit and 12 bit quantization), higher calibration uncertainty (8-11% for MSS versus 7-9% for TM and 5% for ETM+/OLI), and more data artifact issues [21]. The issues and the lack of ancillary data make the calculation of surface reflectance for MSS more difficult. Here, we utilized TOA reflectance for all instruments. To evaluate the impact of the use of TOA reflectance instead of surface reflectance, we downloaded Landsat Collection 1 TOA and surface reflectance products (available for TM, ETM+, and OLI) for a random selection of 40 TM, 30 ETM+, and 20 OLI observations out of the 1430 Landsat observations utilized in the full dataset. Landsat surface reflectance was generated using the Landsat Ecosystem Disturbance Adaptive Processing System (TM and ETM+; [44]) and Landsat Surface Reflectance Code (OLI; [45]). After applying the cloud masks described below, we calculated the average difference between TOA and surface reflectance across the Landsat observation for each band. We then compared the mean difference for each Landsat observation to the drought index value for the scene center of that observation, to investigate the impact of the atmospheric correction on the relationship between weather and reflectance. In addition, we investigated atmospheric correction of all red band observations in the study using the dark object subtraction method [46]. This method can be applied to MSS in a straightforward way, but was potentially affected by a lack of dark objects in the study area (Appendix B).
To filter out clouds and data artifacts (particularly noticeable in MSS data), we created cloud masks for every observation semi-manually [47]. A simple cloud masking algorithm based on the Fmask algorithm [48] and optimized for the study area was used as a starting point, with Fmask overlaid for reference. More details are described in Appendix A (see also [49]). Any obvious issues with the potential impact on the analysis (including instrument artifacts, visible aerosols, cloud shadow, and snow) seen during cloud masking were masked along with the clouds. Landsat 5 observations were observed to have bright artifacts along the swath edges after activation of bumper-mode in 2002 [50]. These artifacts were located north of gaps in coverage due to scan line offsets, so to remove them automatically, we expanded the initial no-data mask 5 pixels further north along the edges of the image. All valid pixels that were not masked were utilized in our analysis. The pixel-by-pixel approach taken here minimized the influence of missing data of any kind, including Landsat 7 ETM+ scan line corrector artifacts [51].
We combined the scene-based Landsat data to produce time series data stacks. Data stacks were produced as half-degree tiles, to keep the data volume manageable within each tile and to match the spatial coverage of the drought index data. Landsat data were downloaded in Universal Transverse Mercator (UTM) and stacks produced in UTM; as the study area crosses the boundary between UTM Zones 13 and 14 (at 102°W), Landsat scenes that crossed the zone boundary were reprojected to the other UTM zone when producing tiles on the other side of the boundary from the scene center. In order for MSS data to be added to the data stack, the 60 m MSS pixels were subdivided into four identical 30 m pixels and transformed to match the rest of the data stack by nearest neighbor resampling. Images were manually inspected for multi-pixel geolocation offsets by visual comparison with other images from the same area, and problematic scenes were removed from the dataset. However, some one or two pixel geolocation errors remained, especially in MSS data and in areas of scene overlap, meaning that there was likely to be some degree of contamination between neighboring pixels that must be taken into consideration during analysis.

Drought Index: SPEI
We selected the Standardized Precipitation-Evapotranspiration Index (SPEI; [23]) to quantify drought in this study. The SPEI is a drought index that is calculated from PET, as well as precipitation; as we used it to investigate wet conditions in addition to dry conditions, we describe what SPEI was measuring as effective moisture. The SPEI is broadly applicable (it uses the same statistical approach in all locations), estimates drought condition based on water balance and PET, and provides for the assessment of drought conditions at any desired timescale. We utilized data from SPEIbase (Version 2.4; [52]). This dataset contains global drought data with half-degree spatial resolution and monthly temporal resolution, in the period 1901-2014 (based on Climatic Research Unit time series Version 3.23 meteorological data; [53]). This dataset is based solely on meteorological station data, so we were able to perform a comparison of two independent measurements (satellite and drought), but the spatial resolution was coarse. Negative SPEI values indicate dry (<50th percentile) conditions; −1 is one standard deviation below the average of the full time series. SPEI values during the study period covered a wide range of drought conditions ( Figure 4).  One strength of the SPEI is that it can be computed over any desired timescale. The SPEI is calculated by comparing the difference between precipitation and PET over a given period of time to what would be expected based on the entire record. For example, the 9 month SPEI for October 1, 2000, would be calculated based on precipitation-PET for the period January 1, 2000 to September 30, 2000, as compared to the distribution of January-September precipitation-PET for all years 1901-2014. To simplify the current study, we focused on 9 month SPEI (i.e., based on a 9 month time period). For the summer months in this study, the 9 month timescale included the current growing season and the previous winter snow season. A shorter timescale would begin to exclude precipitation data from the previous winter and the critical spring precipitation maximum, while a longer timescale would include data from the previous growing season. Furthermore, this timescale is similar to the timescale of the Palmer Drought Severity Index [23]. Preliminary testing demonstrated that other timescales, both shorter and longer, also showed a reflectance response.
Since the SPEI dataset (available monthly, at half-degree spatial resolution) does not have the same characteristics as the Landsat observations, it was necessary to determine how to compare them. Spatially, we utilized the SPEI data for a particular half-degree area for all Landsat pixels contained within that half-degree. This meant that the 20 SPEI pixels present in the study area were each applied to ∼2 million pixels of stacked Landsat data. Temporally, the SPEI is available at monthly intervals, while Landsat observations occur throughout the month. Here, for each Landsat observation that intersected the current SPEI pixel, we linearly interpolated between the drought index at the end of the preceding month and the index at the end of the current month based on the day of the month of the Landsat observation. This method did not account for when in the month precipitation occurred, but potential issues were minimized since in the 9 month SPEI, 8 of the 9 months of input for the SPEI calculation were the same for 2 consecutive months.

Land Use/Land Cover: NLCD
Land use and land cover data at Landsat spatial resolution (30 m) are available from the National Land Cover Database (NLCD; [54]), which generally has over 80% agreement with reference data [55]. We acquired all available years of NLCD data (1992, 2001, 2006, and 2011 [31,56-58]). In order to focus on areas that had stable LULC throughout the study period, we selected only pixels that were in the same LULC category in all NLCD datasets (the 1992 NLCD had somewhat different categories; see Table 2). This was intended as a conservative approach to removing pixels with possible land cover change rather than as a measurement of the amount of change occurring in the study area, since the 1992 NLCD was not directly comparable to later data, and the pixels that were not used (∼30% of the total) should not be interpreted as a measure of change in the study area. A limitation of this approach is that pixels that had LULC change before 1992 would be included in the stable category. Major farm policy changes during this period were of potential importance, and some degree of LULC miscategorization in this period was plausible.

Data Analysis Methods
For each NLCD class in Table 2, we computed histograms of reflectance as a function of SPEI (reflectance bin size 0.001, selected as a number that was well below the expected error in Landsat measurements; SPEI bin size 0.5, selected to produce a reasonable number of categories for viewing). Histograms were computed for the six Landsat bands, including MSS in green and red. Since reflectance varies based on the time of the year, changes through time were investigated by including Landsat observations based on a 40 day moving window. Inclusion was based on day of year rather than calendar date. The 40 day interval was selected because it was the smallest multiple of 10 that included two full 16 day cycles of Landsat observations, and while it served to dampen out short-term variability, it should not otherwise affect the analysis.
Since it is in common use, we computed the Normalized Difference Vegetation Index (NDVI, (near-IR − red)/(near-IR + red)) to provide context and compare reflectance to the vegetation response to drought. As a ratio calculation, NDVI is not directly comparable to reflectance, and it was shown here to aid in the interpretation of the reflectance-effective moisture relationship. While the Landsat sensors have been broadly consistent through time, the spectral response functions are somewhat different for each instrument, which has the potential to affect NDVI calculations. Since the near-IR bands on the MSS instruments have a somewhat different spectral range, we restricted NDVI data analysis to TM/ETM+/OLI. The spectral response of the TM/ETM+ instruments was modeled to produce only a few percent difference in NDVI [59], but OLI NDVI was measured to be ∼10% higher than for ETM+ [60]. Since OLI was only a small fraction of the observations included, here this should not greatly affect the results.
Correlations were calculated using linear regression between SPEI and reflectance. Since the number of Landsat observations per year is heavily weighted toward more recent times, we weighted each year in the regression equally (Appendix C). Simple linear regression assumes that all error is on the y variable (reflectance) and the x variable (SPEI) has no error. The fit slope and significance for simple linear regression tend to be reduced in absolute value when there is error on the x variable [61]. Here, we used simple linear regression values since it is problematic to quantify error due to the calculation of SPEI and since the variability of pixel reflectance changes between land cover types. Furthermore, no adjustment was made for temporal autocorrelation, which might inflate fit significance due to clustering of multiple observations from one particular drought year.
The Interactive Data Language 8.4 (IDL) programming language and ENvironment for Visualizing Images 5.2 (ENVI) geospatial analysis software from Harris Geospatial Solutions were used for spatial analysis, including reprojection. Some statistical work was performed using the R programming language [62], and the geographic information system application QGIS [63] was utilized for visualization.

Results
A correlation between SPEI and TOA reflectance was clearly apparent throughout the study area. In the example pixels shown in Figure 5, much of the variability between observations could be explained by SPEI. Other sources of variability included instrument noise, atmospheric effects, changes in vegetation type, differences in land use/grazing intensity, and other types of weather impacts. A summary of the results from the SPEI versus reflectance regression showed that this correlation was apparent across the study area (Table 3). In Table 3, the slope, intercept, and root mean squared error (RMSE) columns give the per-band results of the regression averaged across all pixels in the study area, and the significant pixels column indicates the percent of pixels that have regression results that are significant at the 0.01 level. The visible and shortwave IR bands were negatively correlated with SPEI (i.e., during high SPEI periods, reflectance was low), while the near-IR band was positively correlated with SPEI. This result was reasonable because the loss of vegetation and decreased soil moisture increased visible/shortwave reflectances, but produced decreases in near-IR (and NDVI, which is driven 60% by near-IR; [64]). Most of the pixels in the study area had a relationship with SPEI. For example, regression of SPEI versus red reflectance produced a significant slope in ∼90% of the pixels in the study area.
Not only was the relationship between effective moisture and reflectance detectable, but its magnitude was large enough to be important. Based on the linear regression of the red band, the expected difference of reflectance between very wet (SPEI = 2) and very dry (SPEI = −2) conditions was 0.05, almost 50% of the average reflectance of 0.11. The estimated difference based on regression was consistent with the mean reflectance of all pixels with SPEI < −1.75 or > 1.75 (0.14 and 0.09, respectively).
There was a small magnitude, but significant relationship between the mean impact of the Landsat atmospheric correction and SPEI ( Figure 6; the last two columns of Table 3). As expected, the atmospheric correction produced a significant (p < 2 × 10 −16 ) change in average reflectance for all bands, decreasing reflectance in the visible bands and increasing reflectance in the IR bands (last column in Table 3); the largest impact was on the blue band, for which the atmospheric correction decreased by ∼0.050 on average. The relationship between the atmospheric correction and SPEI was small in magnitude (the largest was SWIR2 with 0.003 as compared to −0.021 for the SPEI reflectance relationship), but statistically significant (p < 0.005 for all bands, except near-IR, with p ≈ 0.08). In visible bands, the atmospheric correction tended to be larger when conditions were wet (positive SPEI). It was plausible that the humidity was higher and clouds more common, requiring more atmospheric correction during the wetter periods. This suggested that the slope of an SPEI versus surface reflectance regression would be slightly steeper than that of SPEI versus TOA reflectance. In IR bands, the atmospheric correction also had the potential to increase the slope of the SPEI reflectance relationship, due to a larger decrease at dry SPEI values. Use of the dark object subtraction atmospheric correction method also had a small magnitude impact on the slope of the SPEI-red reflectance relationship, but in the opposite direction, perhaps due to the influence of dust or to finding a more appropriate dark object within cloud shadows during wetter periods (Appendix B).    While the relationship between the atmospheric correction and SPEI was statistically significant, its magnitude was small enough to limit its importance for the interpretation of the effective moisture-reflectance slope. For example, the atmospheric correction produced a potential change in the slope of the SPEI-red reflectance of 0.001, but the slope was substantially larger in magnitude (−0.013). The magnitude of the change in the intercept was large enough to be important for interpretation (especially in the blue band, with an intercept of 0.116 and a potential change of 0.050), but here, we focused on the slope of the relationship.
The impact of effective moisture on reflectance varied over the course of the growing season. As others have previously observed [6,65], we found that the land surface was less responsive during early spring and fall, when vegetation activity was low at any moisture level (Figure 7). Effective moisture had a substantial impact on reflectance throughout the summer, as shown by the large differences between the solid red and dashed blue lines in Figure 7. These temporally varying impacts produced differences between dry and wet years in the temporal pattern of reflectance (in visible wavelengths); during dry years, red reflectance decreased only slightly in spring and then had a late summer peak, while in wet years, red reflectance decreased substantially through May into June and was followed by increasing brightness through to the end of the summer (Figure 7). Since chlorophyll absorption reduces visible reflectance, the decrease of red reflectance through the summer was presumably due to a combination of vegetation senescence/seed production and changes in the reflectance of the underlying soils. For context, the pattern of NDVI is shown in Figure 8; it was expected to be dominated by vegetation and was similar to red reflectance, but in the opposite direction.  Figure 7. Average red reflectance by day of year, drought intensity, and LULC. Each point on the lines shows the average red reflectance for all pixels in the study area/period that are within a 40 day window of the day of year, the specified SPEI range, and the specified LULC category. The solid red line is the driest category (SPEI < −1); the short dashed maroon line is drier than average (SPEI between −1 and zero); the dashed dotted indigo line is wetter than average (SPEI between zero and one); and the dashed blue line the wettest (SPEI > 1). Note the different scale for the barren category. The observation count in each point of this plot is available in Figure A4. It makes sense that vegetation activity continues later into the season during wet periods in a semi-arid environment, where vegetation is limited by water availability. Even in drought years, there is often an appreciable amount of precipitation during the winter and spring, which is then available for vegetation growth early in the growing season, but vegetation is only able to continue increasing productivity into the hot summer when there is still a substantial amount of moisture available given increased evaporative demands at higher temperatures.

Effects of LULC and Effective Moisture on Spectral Response
Grassland, the most common LULC class in the study area (∼72% in 2011 NLCD), was observed to be highly responsive to effective moisture ( Table 2). The SPEI versus red reflectance relationship had the steepest average slope (−0.014) of the classes considered, and ∼99% of grassland pixels had a significant relationship. Compared to other cover types in the study area, grassland response was also quite uniform in time and space, with an average RMSE of 0.016 for individual pixel fits and relatively small variability across the study area (Figure 9). Cultivated crops also occurred frequently across the study area (∼17% in 2011 NLCD), but were less responsive to effective moisture and more variable than grassland. The red reflectance average slope of crops (−0.011) was lower than grassland; the mid-summer reflectance difference between dry and wet conditions was less than that found in grassland ( Figure 7); and ∼62% of pixels with a significant relationship were the lowest among the LULC classes, highlighting the variability between pixels. Cropland was also variable within individual pixels, with the highest average RMSE for pixel fits (0.030). The SPEI provided less predictive power for cropland than for other categories (R 2 of 0.19); agricultural field management activities like plowing greatly affected the land surface, but did not occur at exactly the same time in every field in every year, so management decisions could have large effects on individual pixels. During the wettest periods (SPEI > 1), the red reflectance minimum was delayed from mid-June to the beginning of July (Figure 9), perhaps due to late planting (for example, in years with SPEI > 1, the median of spring wheat was ∼12% complete in the last week of April versus ∼50% in wetter years; [66]). This was one of the few times where the average reflectance during a wetter period crossed above that of a drier period.
The next most common LULC class was pasture/hay (∼6% of study area in 2011 NLCD), with results similar to crops, but with generally lower red reflectance (average intercept 0.09 versus 0.11 for crops) and higher NDVI ( Figure 10). Furthermore, average RMSE (0.022) was lower than in cropland, but higher than in grassland. These observations were reasonable if well-watered land with good soil tended to become pasture where there were steeper slopes that were difficult to cultivate and tended to become cropland in large upland areas with a flatter topography. The other LULC classes were less common in the study area (combined ∼5%, none individually >2%). Although all still had a large number of pixels (forest, with the least, had >3000), pixels of these classes tended to be more isolated (>90% of pixels in each class had a neighboring pixel of a different class) and so were more likely to be affected by sub-pixel variability and geolocation errors. Since grassland is the most common LULC class, these less common classes potentially contain substantial numbers of pixels that are partially grassland. Furthermore, these isolated pixels might not be representative of larger areas of that LULC class; for example, the riparian forest in the study area might behave differently than more extensive forested areas. However, it is still interesting to analyze these LULC categories, as they display some different patterns than the three more common classes.
Open water was the only class with a non-significant SPEI versus NDVI fit; NDVI values were near zero throughout the growing season, as expected. However, red reflectance had a shallow but significant fit (slope of −0.008). This result could be attributed to drying and water turbidity during drought, because bare soil and turbid water are more reflective than clear water. The Shadehill Reservoir (the largest water body in the study area, approximately 45°45 N, 102°15 W)) could lose more that 20% of its area during droughts [67,68], and smaller reservoirs and stockponds were likely to be even more affected by drought. Water bodies that dried out during drought were not necessarily filtered out of the dataset with the NLCD mask since none of the NLCD mapping eras (1992,2001,2006, and 2011) featured a serious drought in the study area (nine month SPEI showed drought in late 2006, but this late 2006 drought was not apparent in longer, e.g., two year, SPEI, and the Shadehill water levels were not especially low in 2006; [68]).
Developed land in the study area had a phenological and effective moisture response pattern similar to that of grassland, although fewer pixels had a significant relationship between effective moisture and reflectance (∼80% versus >98% for grassland). This was reasonable, since developed land is a mixture of structures, trees, and grass [54].
Barren land had by far the highest average red reflectance among the LULC types (0.17; Table 2), due to limited vegetation coverage. Unlike open water, barren land pixels tended to have low, but non-zero NDVI values; inspection of images in Google Earth showed that these pixels often contained small patches of grass or shrub. The slope of the SPEI versus red reflectance fit was fairly steep (−0.012) compared to the SPEI versus NDVI slope, suggesting that some of the reflectance relationship was due to darkening of the bare soils by moisture. This was also supported by the fairly constant value of reflectance through the summer. The high average RMSE (0.026) and low fraction of significant pixels (∼75%) demonstrated that this relationship was less well predicted for bare soil than, e.g., grassland, perhaps due partly to soil color responding more quickly than vegetation condition to precipitation.
Wetlands have a shallow red reflectance average slope (−0.008) compared to other LULC types, despite high NDVI. Wetland vegetation might be buffered due to more stable moisture availability, even during drought. A stronger response might be seen during a decadal drought (such as the 1930s). Due to their lush vegetation and the presence of water, wetlands had the lowest average red reflectance (0.077). As in the open water class, the absence of surface water during droughts might increase red reflectance. The phenology of wetlands includes a steep increase in vegetation in the spring (especially May), which is expected for the almost exclusively C 3 wetlands species, during which time, effective moisture has no apparent effect on NDVI.
The phenology curve of scrub/shrub was similar to that of grassland, with a slightly less prominent minimum in June, but with lower reflectance continuing into the fall (particularly during drought), which produced a slightly less steep red reflectance slope in the late summer period tested (−0.012; grassland: −0.014). In both grassland and scrub/shrub, the driest NDVI curve was nearly flat, with little spring green-up. In the study area, the relatively small number of scrub/shrub pixels was commonly isolated in small clusters interspersed with grassland; if these pixels often contained some grass, it was not surprising that they were similar to grassland. Shrubs tend to be deeper rooted than grass, explaining the continuation of low reflectance into the fall.
Forest was not well represented in the study area, since effective moisture was not usually high enough to support tree cover. The pixels of forest that did occur were mainly in river valleys and near watercourses. The average slope for SPEI versus red reflectance was shallow (−0.010) compared to grassland; however, pixels had low average RMSE (0.014), and most (∼99%) were significant.

Fine-Scale Spatial Variability in the Impact of Effective Moisture
The slope of the effective moisture-reflectance relationship displayed detailed patterns when shown for each 30 m pixel across space (Figure 11). Impacts of LULC were visible as expected, with grassland pixels having generally steep slopes compared to cropland, water, or forest. However, within-class variability was also visible. For example, while many bodies of water had a low slope, smaller water bodies often had a steep slope. Small water bodies in the study area could be greatly reduced in size and depth during drought, increasing their reflectance and producing a strong relationship between effective moisture and reflectance. Similarly, the Cheyenne River was notable both for carrying a substantial sediment load [69] and for having a strong effective moisture-reflectance relationship. 46  Patterns within grassland are of particular interest since grassland often has a strong effective moisture-reflectance relationship. For example, the Moreau Prairie ecoregion (cyan outline in Figure 11a) had lower productivity than surrounding areas because of the alkaline soils derived from the underlying Hell Creek Formation [70]. While the difference between the grassland in the Moreau Prairie and in surrounding areas was not clearly apparent in visible reflectance (Figure 11c), the central areas of the ecoregion appeared to have a less strong relationship than surrounding areas of grassland (Figure 11a). This might be attributable to lower productivity (and higher red reflectance) even in wet periods, leading to less change between wet and dry periods. Patterns in the effective moisture-reflectance relationship also followed management boundaries, topography, and soil type in some places, demonstrating that there was information available to be drawn from the slope.
The image of RMSE also displayed patterns across the land surface. Croplands stood out as locations of substantial variability, as did some bodies of water, including the Cheyenne River (visible as a dark blue line across the southern portion of the study region in Figure 11b). RMSE contained interesting within-class variability, for instance croplands commonly had high RMSE, but varied between management areas. In grasslands, one source of variability was locations having prairie dog activity at some point in the study period, which showed larger RMSE than surrounding locations.
One other advantage of this image approach was that it had the potential to highlight data inconsistencies. Although there was a sharp boundary in the SPEI values used at the half-degree latitude/longitude lines, there was little difference in slope across the boundaries (the vertical line at 102°W was due to reprojection across the UTM boundary). This suggested that the gridding of the SPEI values was not producing inconsistencies between the SPEI pixels. There was a small amount of striping due to Landsat 7 scan line corrector artifacts, and there were features at the edge of Landsat scene boundaries in the southwestern portion of the study area. However, the slope images were generally consistent between scene center and overlap zones, despite approximately double the number of observations in overlap zones and no bidirectional reflectance distribution function (BRDF) correction between the paths. In the study area, Worldwide Reference System-1 (WRS-1) overlap zones broadly covered the area of the WRS-2 scene center, so the relative observation frequency also changed through time, but was not visible in the slope or RMSE images.

Discussion
The relationship between effective moisture driven by weather and LULC (i.e., high visible/shortwave-IR reflectance during drought) was consistent with expectations and with previous work, e.g., [5,10,71]. During a drought, the vegetation senesces and loses its leaves, while the soil dries. These processes combine to increase visible/shortwave-IR reflectance. Near-IR reflectance is high for actively growing green vegetation and decreases with the loss of active vegetation; however, similar to visible and shortwave-IR, near-IR reflectance of dry soils tends to be high [13], suggesting that vegetation was dominating the near-IR reflectance signal. The temporal pattern was also as expected, with a larger impact from effective moisture during the growing season. Since the visible and shortwave-IR reflectance bands were observed to increase during drought, albedo could be expected to increase also, although the near-IR wavelengths would counteract this relationship.
The robustness of the relationship between effective moisture and reflectance presented here was striking; it was apparent in nearly every regional, temporal, or LULC subset that was investigated. The study area was selected as one that would likely respond to effective moisture, but even so, the wide applicability of the relationship was noteworthy. The 30 m resolution was detailed enough to show that this relationship was not just a result of averaging over a large area, but was relatively consistent across a variety of soil, vegetation, and management conditions. The magnitude of the relationship was not greatly affected by atmospheric correction, although the average reflectance did change.
The uniformity across the study area in the response of grassland to effective moisture (grassland/herbaceous in Figure 9) was quite interesting. Most of the 30 m grassland pixels in the study area responded to effective moisture in a similar way, indicating a robustness to different levels of grazing. The higher variability seen in barren land suggested that uniformity in soil reflectance was not the main driver of this uniformity and that vegetation was an important factor. Plant communities were not uniform across the study area, with variability in the fraction of C 3 and C 4 grasses [34]. This variability in vegetation suggested that different communities of grass species across the study area responded to effective moisture in a way that was broadly similar, which was potentially related to evolution in a highly variable environment with grazing pressure.
Not using time as a variable in the analysis highlighted the effective moisture-reflectance relationship by land cover type and reduced the influence of temporal biases on the analysis. For example, the full Landsat record was cross-calibrated, but if the calibration of one of the early Landsat instruments produced data that were slightly lower than that from later instruments, a time-focused analysis might find that there was an increase in reflectance even if nothing had changed on the land surface. However, the measurements from that earlier instrument would cover a range of drought index values and thus should have a more limited impact on the slope of the relationship between effective moisture and reflectance.
It is useful to note that large changes in reflectance during drought do not necessarily imply that drought is producing changes that will be long lasting or difficult to reverse once precipitation returns. For example, when grass cover senesces during drought, it can return to normal by the next growing season if precipitation returns, e.g., [72], whereas forest cover that has more stable reflectance might have substantial tree mortality due to drought, which continues to affect the character of the landscape even after the return of rain [73,74]. The strong effective moisture-reflectance relationship observed in this study suggested that much of the observed change in reflectance was reversible; if the land surface was undergoing (relatively) permanent changes during drought, then reflectance values might not return to their previous levels, and there would be more noise in the effective moisture-reflectance relationship. In a landscape that was out of equilibrium, the strong relationship found here between effective moisture and reflectance might break down, providing a potential method for locating such areas.

Potential Applications for Spatial Variability in Effective Moisture Impacts
A spatial perspective on the variability in the effective moisture-reflectance relationship potentially provides interesting tools to better delineate the characteristics of the land surface. For example, impacts related to grazing intensity might be investigated by comparing the slope between areas with different management regimes. One potentially interesting focus of study is variability between soil types, as highlighted by the Moreau Prairie ecoregion. This type of approach might provide additional information to help map small areas of differing soils, using freely available remote sensing data. Residual differences from expected values also might be used for interpreting differences through time, especially if higher resolution moisture data were available. The slope data from this study provided a time-averaged result, so that small regions of variability in the short-term (e.g., opposite sides of a field receiving different amounts of precipitation in a thunderstorm) had limited influence on the slope results.
The 30 m spatial scale allowed the possibility of detecting important patterns that might be obscured within coarser resolution satellite observations such as from Moderate Resolution Imaging Spectroradiometer (MODIS), Visible Infrared Imaging Radiometer Suite (VIIRS), or Advanced Very High Resolution Radiometer (AVHRR). For example, the difference we observed between small water bodies and their surroundings would be averaged out in imagery with coarse spatial resolution, losing the pattern of the water and reducing the slope seen in the surrounding grassland. The 30 m scale was also helpful for distinguishing patterns within classes, for example drainage patterns or irrigation features that varied within an agricultural field, or the difference in moisture response between a north-versus a south-facing slope. Even differences that were small relative to the range of variability could be observed in the slope images.

Using Reflectance and Weather Data in Modeling
Modeling of the relationship between moisture and surface reflectance in time series remote sensing data can be utilized in algorithms for change detection, e.g., [75,76]. For change detection, it is useful to understand and to be able to model processes that cause variability in satellite reflectance measurements. Quantification of the moisture-reflectance relationship has the potential to improve the ability to distinguish the effect of moisture on reflectance from other processes occurring on the land surface (e.g., clearing of land for development).
An advantage of linear regression is that it is a useful approximation that would be tractable to implement in a model of land-atmosphere feedback. In climate and weather models, albedo is often estimated based on biome and is set as a constant through time, e.g., [77,78]. To simulate feedback, albedo can be estimated from vegetation and soil fractional coverages modeled based on climate, e.g., [79]. Another approach is to utilize satellite measurements to estimate albedo, either by measuring albedo across an average year and using that for all years or in retrospective studies of recent time periods, e.g., [5,80,81]. However, it is challenging to measure albedo on near-real timescales [82] and impossible to measure albedo before the satellite record or in the future. With a simple model for the response of albedo to weather, it might become possible to improve albedo estimates beyond a constant value. Although the relationship between reflectance and effective moisture is likely to be different in other study areas and the relationship between reflectance and albedo is not straightforward, the robustness of the effective moisture-reflectance relationship presented here suggested that there was a possibility for using a drought index to model improved estimates of albedo.
Here, we modeled the relationship between effective moisture and reflectance using linear regression. Qualitatively, the plots of SPEI versus reflectance appeared to be somewhat linear ( Figure 5). However, there was no mechanistic reason to assume the relationship had to be linear. In study areas less limited by water, a non-linear response is plausible; a rainforest, for example, might have increased visible reflectance (due to decreased active vegetation) in its driest periods, but not have much reflectance difference between average and wet periods (or even have decreased active vegetation causing increased reflectance during wet periods, due to less available sunlight; [83]). In any study area, reflectance cannot go beyond its theoretical bounds, and so, a linear fit must break down at extreme values of the drought index. Even within this study area, there is reason to believe that a linear relationship is only an approximation; the trend through time of pixels with SPEI < −1 is more offset than the other lines in Figure 7, and the medians of the wetter bins are roughly equal for several land cover types (e.g., cultivated crops) in Figure 9.
In addition, the simple linear regression used in this study assumed that all error was in the y-values (reflectance). This was clearly not the case, as the weather variables used in the SPEI calculation had errors, and even more importantly, converting point meteorological data to gridded data introduced significant errors. Patil [84] estimated a mean absolute error of 0.4 in SPEI for the study area due to spatial averaging. That error was similar to or larger than the 5-10% Landsat instrument errors [21]; however, surface variability will also increase the y-error. Given the extra conceptual and computational complications of a more sophisticated fit, we utilized simple linear regression, but the fit slopes would be steeper if the x-error were taken into account.
Other potential sources of error include change in LULC, the cloud mask, the use of top of atmosphere reflectance, registration, and calibration differences between instruments or over time. These errors will increase the amount of noise in the reflectance measurements, but more concerning would be if a source of error produces bias that affects the effective moisture-reflectance fit. Missed clouds during wet periods would act to decrease the slope measurement, while increased atmospheric water vapor would decrease TOA reflectance and increase the slope. The significance and robustness of the measured fit suggested that there was little cloud contamination, as thick clouds are highly reflective and would be highly significant outliers for the fit. Dust or haze in dry periods would increase the slope. Registration errors could affect pixels with subpixel or neighboring pixels with different LULC, so the uncommon LULC classes largely composed of isolated pixels were likely to be contaminated by other LULC classes (most commonly grassland). A change in LULC might also lead to incorporation of data from one class into another in this analysis, for example a pixel that was cropland early in the study period and then became grassland before 1992 would be included in the grassland data.
Since time is not directly a component in the effective moisture-reflectance relationship, instrument calibration errors and other changes during the study period should not produce a large bias; any effect would be second-order (due to more droughts occurring towards the beginning of the study period). This minimization of time bias was an important advantage of using this method of analysis.

Effect of LULC on Albedo-Drought Interactions
The impact of effective moisture on reflectance can provide insight into the potential for climate-albedo feedback within different LULC categories. Since increases in albedo have been observed to reduce precipitation [5], it is important to understand albedo changes in response to effective moisture. Barren land that has a consistently high albedo, which might inhibit precipitation at all times (e.g., in the Sahara Desert), whereas grassland might be more subject to intensifying drought feedback and have large interannual precipitation variability.
The small variability in grassland visible reflectance for a given drought index shown in this study suggested that drought-albedo feedbacks were relatively predictable in grassland; albedo values would generally be similar after a given amount of effective moisture, and so should produce a similar amount of feedback. Modeling precipitation feedback might be relatively straightforward in grasslands. In contrast, the higher reflectance variability of cropland suggested that the feedback potential in cropland was not consistent throughout the study. Management decisions likely were an important part of this variability, suggesting that management could affect feedback. For example, Vick et al. [85] suggested that discontinuation of summer fallow across a large area of the northern Great Plains in recent years produced increased evapotranspiration during the summer months, leading to increased humidity and precipitation and decreased temperature. Feedback between reflectance and drought might also have importance for forest/grassland and grassland/barren transition locations.

Conclusions
Effective moisture had a substantial effect on reflectance across the northern U.S. Great Plains study area. Most pixels had higher visible reflectance during dry periods than during wet periods, especially in grassland, where the difference was particularly substantial. In cropland, there was more variability in reflectance at a given drought index value than in other LULC categories, although there was still a correlation between effective moisture and reflectance.
Feedbacks between albedo and moisture are an important part of the Earth system, and further research is needed to explore these interactions. The pattern of LULC on the landscape is also potentially important; boundaries between LULC types affect the initiation of convection and precipitation [86]. Using linear regression to estimate reflectance could be a helpful approach in climate models. One important topic for future consideration concerns the boundaries of LULC categories and whether feedback plays any role in maintaining these boundaries.

Appendix A. Cloud Masking Methods
The cloud masking goal for this study was to minimize (to the extent possible) statistical bias in the Landsat dataset due to clouds. Because clouds are highly reflective, they have substantial potential to skew data distributions; we utilized an aggressive cloud removal strategy, even at the expense of removing numerous clear pixels. This approach could be expected to remove more pixels during wet periods (when there are more clouds), but since wet periods are already expected to be less represented due to clouds, discarding some clear pixels during these periods does not add a new bias. However, areas that are consistently detected as clouds during clear conditions have the potential to produce statistical issues not already present in the data; for example, if a barren pixel is preferentially flagged during periods of high reflectivity, both its median and the slope of its relationship with the Standardized Precipitation and Evaporation Index (SPEI) would be expected to be artificially reduced. Thus, we pursued the approach of being relatively conservative in flagging pixels as clouds, but then aggressively buffered the flagged pixels to remove neighboring unflagged clouds; this approach also served to remove most cloud shadows without any other detection strategy.
The starting point for cloud masks production for this study was a simplified version of the Fmask algorithm [48]. It was also partially similar to that of Braaten et al. [49], who described a cloud mask for the Multispectral Scanner System (MSS). In the simplified algorithm, pixels are flagged as clouds if all visible bands are above a reflectance cutoff and all visible bands have a similar reflectance (since clouds tend to be uniformly bright). Cutoff values were selected manually, optimized for this study area to minimize extraneous detections seen around Badlands National Park (where highly reflective barren pixels can easily be flagged) and Rapid City, SD (the largest city near the study area). For MSS, pixels were flagged if red (R) and green (G) were both >0.35 and if G − R was within (−0.02,0.025). For Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and Operational Land Imager (OLI), the test was R, G, and blue (B) all > 0.2 and both G − R and B − R within (−0.03,0.025). Alternately, TM, ETM+, and OLI pixels were flagged if R, G, and B were all > 0.35, since some oversaturated cloud pixels (particularly in TM) had significantly different R, G, or B values. To further reduce false detections, pixels were removed if there were not at least 10 other detected pixels in the surrounding 5 × 5 pixel region. Finally, all flagged pixels were surrounded with an aggressive buffer (20 km for MSS, 10 km for OLI/ETM+/TM), to remove undetected clouds, as well as cloud shadows.
To further reduce cloud errors and instrument artifacts (particularly in MSS), we performed a manual filtering step. The starting point was the simple cloud mask and for MSS, a mask filtering out any scan line that averaged at least twice that of surrounding lines or of the whole image. Results of Fmask (Version 3.3) were overlaid for reference during manual masking. Clouds, cloud shadows, snow, or instrument artifacts that were not detected by the simple mask were manually flagged (clear pixels in the vicinity were also flagged in the interest of analyst efficiency). False cloud detections caused by surface features were manually unmasked to reduce statistical bias (clear buffers around valid cloud detections were generally not unmasked).
The resulting dataset [47] contained a few areas with a low number of observations, but fewer than the Fmask algorithm or simple mask alone ( Figure A1). The manual mask also drastically reduced the number of visible clouds in the dataset (as with any manual procedure prioritizing efficiency, there were a few misses). Manual masking also reduced haze and smoke pixels, although these were more difficult to identify and flag than clouds that were more opaque. While not done for this dataset, the simple masking procedure used as manual input could be further improved by running it a second time with the few problematic/consistently flagged pixels required to never produce cloud detections. This also might allow a slight reduction in the cutoffs used, to reduce the number of clouds missed. 46

Appendix B. Atmospheric Correction with Dark Object Subtraction
We investigated the dark object subtraction method of Chavez [46] for atmospheric correction since it can be computed using MSS data. Dark object subtraction assumes that the darkest pixel(s) in an image provides an approximation of the radiance from the atmosphere, and so, subtracting that value from the image provides an atmospheric correction. Here, we focused on the red band. We found the minimum value for the red band in each Landsat observation, subtracted that from the entire image, and then converted to reflectance by multiplying by the same factor as the conversion to top of atmosphere (TOA) reflectance. Then, linear regression was performed with SPEI for all pixels, following the same method as for TOA reflectance. For each pixel in the study area, the dark subtraction linear regression intercept and slope were subtracted from the TOA intercept and slope, and then, histograms were computed across all pixels to determine how much the fit coefficients changed when the dark subtraction method was applied.
As expected, the dark subtraction correction decreased the reflectance and the SPEI-reflectance intercept (by 0.022 on average; Figure A2a). However, the effect on the slope of the SPEI-reflectance fit was small (average increase of 0.0003 compared to the average slope of −0.013; Figure A2b). This suggested that the dark subtraction correction did not have a large impact on the analysis of the SPEI-reflectance relationship.
A challenge for dark object subtraction in this study area was the lack of dark objects. Dense vegetation, deep clear water, and shadowed slopes are examples of dark targets, but none were common in the study area. Shadows from clouds are also potential dark targets, but not all observations included cloud shadows. Qualitatively, we observed that water was often the minimum used by dark object subtraction (at least in observations without cloud shadow). However, water turbidity was a potentially serious complication as it was common in the study area. In addition, we did not use the cloud masks before searching for a dark subtraction value in order to include cloud shadow as a possible source of dark subtraction. This allowed for the possibility of data artifacts producing the dark subtraction value, which we saw in some observations. To further investigate the relationship between dark subtraction and SPEI, we compared SPEI to the magnitude of the dark subtraction correction. For each of the 1430 Landsat observations in this study, we found the SPEI at the scene center and the difference between TOA reflectance and dark subtraction reflectance for the red band. In drier conditions (negative SPEI), the dark subtraction correction was slightly larger than during wet conditions ( Figure A3). This was the opposite direction from the trend seen between SPEI and the surface reflectance correction, where the correction was larger in wet conditions ( Figure 6 in the main text). It was possible that this was an accurate description of the atmospheric correction, i.e., that there was more atmospheric dust during dry periods, and so, a larger correction was needed. Alternately, it was possible that this trend was an artifact of the dark subtraction correction. If cloud shadow were more common in wet periods, the dark subtraction correction might be more accurate during wet periods and tend to produce an overcorrection during dry periods. It was also possible that the dark subtraction correction was higher due to the land surface effects we were trying to measure, or in other words, that the surface reflectance of dark targets was higher during dry periods than wet periods. Dark vegetation is likely to have higher reflectance during drought, and stock ponds, wetlands, and reservoirs are likely to have increased reflectance during drought when they become more shallow or dry completely. Thus, it would be problematic to use dark subtraction to investigate the SPEI-reflectance relationship. While it was possible that more atmospheric correction could be needed in dry conditions due to increased aerosol optical depth, it was difficult to disentangle that from the land surface effects, and the surface reflectance correction-SPEI relationship suggested that the land surface effects might be more important.

Appendix C. Weighted Regression
Since there have been many more Landsat observations in recent years than in the early part of the record, we calculated the linear regression fits so that each year in the record would have equal weight. To do this, we weighted each observation by 1/m i , where m i is the number of Landsat observations in the year of observation i. If the model for regression is: where y is the reflectance, x is the SPEI drought index, is the residual error, and β 0 and β 1 are the fit solutions, then we maximize the likelihood L such that: where N is the total number of observations. Solving for the β values,