## 2 Materials and Methods
### Satellite Data Processing  
In accordance with [PERSON] and [PERSON]'s [32] validation of STARFM for the same Landsat Path/Row, Landsat TM data were used as the high spatial resolution images. The choice of high temporal resolution was based on previous findings [39] that favored the MODIS Bidirectional Reflectance Distribution Function (BRDF) adjusted 16-day reflectance (MCD43A4) product over MODIS 500 m daily surface reflectance (MOD9 GA) data. Therefore, we chose to use the nadir BRDF-adjusted reflectance (NBAR) data with an eight day overlapping temporal resolution from the MODIS collection 5 data products. USGS's EarthExplorer website provided the download for the Landsat TM reflectance images from the Climate Data Record (CDR) collection. Two different methods were used to select the Landsat images used in the base pairs. The first used an informal screening process, similar to [PERSON] and [PERSON] [32], where we visually inspected Landsat scenes to ensure that there was relatively minimal cloud, shadow, snow, or water in our study area and immediately adjacent to it. To do this, we downloaded 80 images for path 38/row 31 with less than 10% cloud cover spanning the period from the start of the MODIS data production until the end of Landsat 5 (Day 49, 2000-Day 361, 2011). The CDR data also included identification for clouds, cloud shadows, snow, and water [40]. Python scripts were written to automate the process of preparing Landsat data for STARFM by looping through the images at each step, relying on the ArcPy library [41] for many steps. The first step was to extract the files and rename the Red, Near Infrared (NIR), and Quality Assurance (QA) data for sorting and pairing when defining the STARFM inputs. The data were then reprojected to the WGS84 UTM_Zone12N coordinate system and clipped to the same extent. Based on the CFmask QA band that comes as a part of the CDR data, the scripts masked out all areas that represent fill, water, shadow, snow, and cloud, removing them from the analysis. Last, the scripts converted the images to a generic 16-bit signed binary file. Figure 2 shows a schematic that details the flow of the Python algorithm. Processed images were visually inspected under the informal screening criteria to remove those with high levels of no data values in and immediately adjacent to the study area, which left 54 Landsat images available for analyses. Removing a processed Landsat scene does not reduce the number of images used in determining seasonality parameters, instead it adds in an additional prediction day by removing a base pair. Figure 2: The steps contained in the python script to process the MODIS and Landsat data for input into STARFM. The second Landsat selection criteria used a statistical approach similar to [PERSON] et al. (2017) where the selection was based on the percent of unclear pixels. To do this, for all available Landsat surface reflectance images we calculated the proportion of clear pixels within a 1500 m buffer of our study area, which represents the entire area that STARFM used for the search window. To determine the total number of clear pixels, we reclassified the clear flag in the CFMask file from the Landsat Surface Reflectance product so that clear pixels are set to a value of 1 and snow, water, cloud, and cloud shadows are set to 0. We then preformed zonal statistics on each reclassified CFMask that summed up all pixel values within our buffer area, considering no data pixels as not clear. To determine the total number of possible pixels, we created a raster image with 30 m resolution that encompassed the study area plus the 1500 m buffer and obtained a count of the pixels (504,714 in total). We then divided the total number of clear pixels by the total count and multiplied by 100. We set two statistical thresholds of the amount of clear pixels needed to include the image, one at 90% and the other at 95%. The same Python script then prepared all images that met either of these criteria for STARFM.