IMPACTS OF BEETLE KILL ON MODELED STREAMFLOW RESPONSE IN THE NORTH PLATTE RIVER BASIN OF BEETLE KILL ON MODELED STREAMFLOW RESPONSE IN THE NORTH PLATTE RIVER BASIN.”

: A beetle epidemic across the western United States has resulted in the death of millions of acres of forests. This beetle outbreak, referred to as “beetle kill”, has caused many to believe that such dramatic changes in land cover could potentially alter the hydrology of the impacted regions. One of the most important hydrological processes that beetle kill has the potential to impact is streamflow. This research evaluates the hydrologic impacts on streamflow from land cover change due to beetle kill in the North Platte River Basin (NPRB) (Colorado and Wyoming, USA) by utilizing the Variable Infiltration Capacity (VIC) hydrologic model. Utilizing the National Agricultural Imagery Program (NAIP) dataset from 2005 / 2006 (onset of “beetle kill”) to more current conditions (2009), a decrease in tree cover of 16% to 40% was estimated. This decrease in tree cover was applied to VIC modeled streamflow from 1950 to 2000. The VIC model predicted a minimal increase in streamflow of approximately 5% which was not statistically significant .


Introduction
Many coniferous forests across western North America are experiencing a vast epidemic of bark beetles that has caused widespread tree mortality to peak to surprisingly unprecedented levels. The current epidemic, which began in North America roughly around the mid to late 1990s, has claimed billions of trees ranging from parts of northern Mexico all the way to the southern Canadian provinces of Alberta and British Columbia [1] [2] [3]. Since its beginning in the United States, this outbreak of bark beetles has impacted roughly 42 million acres of affected forests, which compares closely in magnitude to the entire state of Oregon [4]. The widespread pandemic has become so extensive throughout much of the West, that many locals have coined the phrase of "beetle kill" to describe the process the bark beetles employ to kill off the trees and eventually, given the appropriate amount of time, entire stands of forests. Beetle kill impacts include falling trees [5] [6] [4], wild fires [7], shifting ecosystems [8], reduction in tourism and skiing [9] [6] [4] [8], and homeowners costs associated with clear cutting [9].
Many researchers have come to realize that such large-scale tree mortality could have implications on the water resources of many beetle infested forest regions [10] [11] [12] [5]. The more forests that are lost due to beetle kill dramatically increases the chances for flash floods throughout the West because there are no longer any trees in place to catch the snow and attenuate the snow melt [5]. Flash floods create the possibility of more stream bank erosion, which decreases water quality by introducing more sediments into streams [10] [6] [8].
Of all the regions severely affected by beetle kill across western North America, no area has been more greatly affected than northern Colorado and southern Wyoming [13] [6]. Since the late 1990s, the combined amount of coniferous forests affected by bark beetle in the states of Colorado and Wyoming have totaled approximately 10 million acres [4]. The particular species of bark beetle that is causing so much devastation throughout much of the Rocky Mountain region is the Mountain Pine Beetle (MPB) (Dendroctonus ponderosae) [7] [14] [8]. It has been estimated that since 2008, in the states of Colorado and Wyoming, the MPB has affected nearly 5 million acres of lodgepole pine forest alone; this is an increase of roughly 25 million trees from the baseline year of 2006, in which only 2.6 million trees were infected [15]. Although the MPB has caused a great deal of devastation to forests throughout much of the West, MPBs have always been a native species to the pine forests of western North America. According to historical records from the past century, outbreaks of MPB have always occurred throughout the West but never to the extent to which levels are currently being observed [17]. Many scientists have speculated that the underlying cause for the massive beetle outbreak can be largely linked to climate change and other associated abiotic factors [2] [14] [5] [17] [18].
The current research utilized high resolution remote sensed data to assess the percentage of tree loss (with error) due to beetle kill. Land cover was then changed in a hydrologic model to reflect the loss of tree cover and the change in yearly streamflow volume (with error) was observed.

Study Area
The area of study is the North Platte River Basin (NPRB), which is located in the states of northern Colorado and southern Wyoming (Fig. 1). The location of the NPRB lies between the Sierra Madre range to the west and the Medicine Bow range to the east. United States Geological Survey (USGS) streamflow gage 06630000, which is located above Seminoe Reservoir, the first in a series of impoundments in the state of Wyoming, was identified as the gage of interest. The NPRB has a drainage area of approximately 10,800 km2 (4,175 mi2) and the watershed has an elevation range from roughly 3,650 meters to 1,950 meters (12,000 feet to 6,400 feet) with varying terrain consisting from rugged mountains to gently sloping grasslands. Of all the macro-scale watersheds in the Colorado/Wyoming region, the NPRB is one of the most severely affected by the recent beetle kill infestation. Presently, the NPRB, which was severely hit with the MPB epidemic in 2007/2008, has seen virtually all the lodgepole pine forests completely annihilated by beetle kill (Pendall et al., 2010).

NAIP Imagery Data
Scale appropriateness of evaluation imagery was an important consideration when choosing an imagery dataset. The resolution of imagery selected for assessing should be appropriate for the features to be reviewed [19]. O'Neill et al. (1996) [20] suggests a resolution of one-fifth to onehalf the size of the features of interest. Woodcock and Strahler (1987) [21] recommends an image resolution of one-half to three-fourths the size of target objects. Because of the density of forest cover and the irregular locations of affected trees within the NPRB, in order to identify forest cover that has been affected by beetle kill within a dense area of unaffected trees, a high resolution imagery dataset was needed. In addition to this, the approximately 4175 mi² extent of the NPRB adds another constraint and narrows the choice of datasets even further.
From a macro-scale basin level perspective, assessing forest loss not only requires high resolution imagery, but imagery with extensive spatial and temporal coverage, as well. Since high resolution remote sensing data for land cover classification, such as LiDAR, are not always readily available or cost effective, especially for this particular project, an alternative high resolution dataset was found with the use of the National Agricultural Imagery Program (NAIP). NAIP is well suited to serve all these needs and is regularly used for land cover classification [22] [23] [18] [24]. NAIP provides one meter ground sample distance (GSD) ortho imagery rectified to a horizontal accuracy of within +/-5 meters of reference digital ortho quarter quads (DOQQ's) from the National Digital Ortho Program (NDOP) with coverage spanning the vast majority of the United States. The imagery is collected in natural color during the agricultural growing season and contains no more than 10% cloud cover. Since NAIP imagery is collected on a state by state basis, as a consequence, imagery for one state may not be available for the same year as a neighboring state. This was the case for Colorado and Wyoming. Beginning in 2003, NAIP was acquired on a 5-year cycle, with imagery being collected for Colorado in 2005 and Wyoming in 2006. Since these years were both before the beetle infestation became widespread, any differences in classification that might have occurred before this time period were considered negligible. However, NAIP began a 3-year cycle in 2009; resulting in complete coverage of both Colorado and Wyoming, and thus the entire NPRB.
The NAIP imagery were obtained from the USDA Geospatial Data Gateway website. These imagery are made available as either compressed county mosaics or uncompressed 3.75 minute by 3.75 minute quarter quadrangles with a 300 meter buffer on all sides. For labeling purposes, the period spanning up to and before 2005/2006 for the NPRB can be considered the "pre-beetle kill" era, whereas the period spanning up to and after 2009 for the NPRB can be considered the "postbeetle kill" era. Referencing these two periods was not completely accurate since beetle kill did not truly manifest in the NPRB until roughly 2007/2008, but for the purposes of labeling ( Figure  2a & Figure 2b) it was adequate. The grey ( Fig. 2a) portions of the NPRB represent living unaffected forests, whereas the black (Fig. 2b) represent dead affected forests. It can be seen from the figures that the bark beetle epidemic has severely impacted the forested regions.

Hydrological Model Data
The hydrologic model used in this study is the Variable Infiltration Capacity (VIC) model [25] [26] [27] [28]. In order to run successful model simulations, VIC requires a variety of meteorological forcing data (precipitation, temperature, wind speed) from an array of credible sources [28]. All daily meteorological forcing data were obtained from the Soil and Water Modeling Group, University of Washington (www.hydro.washington.edu/SurfaceWaterGroup/data.php) [29]. One of the most important datasets for this analysis was the land cover classifications as described by Hansen, DeFries, Townshend, and Sohlberg (2000) [30], which were based on and obtained from the Department of Geography, University of Maryland (http://www.geog.umd.edu/landcover). This data contains different land cover types all of which were at a one kilometer spatial resolution and consists of 14 total different land cover classes. For this analysis, the land cover in the NPRB did not consist of all 14 possible land cover classes that VIC specifies. Since this was the case, it was concluded that only 11 of the 14 land cover classifications would apply to this analysis. Table 1 below reveals the 11 land cover classifications that were used to classify the entire vegetative coverage of the NPRB. Wooded Grasslands 8 Closed Shrublands 9 Open Shrublands 10 Grasslands 11 Crop Land (Corn)

Landcover Analysis (Maximum Likelihood Classification)
The accuracy of classification results is determined by which classification method is selected. Generally, landscapes with homogeneous spectral responses are estimated through a multivariate Gaussian distribution and the assignment of pixels to classes is often based on a maximum likelihood classification (MLC) [31]. Of all the parametric classification techniques, MLC is the most commonly used for classification of forest landscapes [32]. Being consistent with previous literature, a MLC technique is used to classify forest cover using the interactive supervised classification tool in ESRI's ArcGIS 10.0. This tool uses an MLC algorithm to classify imagery based on a signature file created from user derived training samples. The MLC algorithm is based on two principles; the cells in each class sample in the multidimensional space follow a Gaussian distribution and Bayes' theorem of decision making [31].
The MLC computes a variance-covariance matrix for each of the class signatures when assigning each pixel to one of the classes represented in the signature file. Under the assumption of normality a class sample can be characterized by the mean vector and the variance-covariance matrix [33]. Given these two characteristics, a statistical probability is computed for each class to determine the membership of each cell to a class within the signature file. An equal a priori probability weighting option is specified, allowing all classes to have the same a priori probability; resulting in each cell being assigned to the class to which it has the highest probability of being a member [34].

Hydrologic Model
Since its inception nearly 20 years ago, VIC has undergone a variety of updates and has been extensively used in topics focused primarily on water resource applications ranging from climate  [32] change to land use change studies [35]. Specifically VIC is a macro-scale, physically based, semi-distributed, land surface hydrologic model. The model operates by using 1/8° spatial resolution using gridded meteorological forcing data (precipitation, max/min temperature, and wind speed) in conjunction with other watershed characteristic data (land cover, soil, elevation bands, snow bands, etc.) in order to estimate surface water runoff and base flow [45].
For this analysis, all simulations were performed using the VIC model (version 4.1.1) downloaded from (http://www.hydro.washington.edu/Lettenmaier/Models/VIC) and were completed by using VIC's full energy and water balance mode of operation at daily time intervals with 1/8° spatial resolution. The VIC model used for the current research has been utilized in two previous research efforts [45]. Similar to previous modeling efforts of Acharya, Piechota, Stephen, et al. (2011) [45], the routing model was not used for this analysis due to the small basin size (only 95 grid cells) and the need to observe annual streamflow volume.

Land Cover Classification Results
Aerial images for the NPRB were taken in 2005/2006 and again in 2009. For classification purposes, 2005 (Colorado) / 2006 (Wyoming) was defined as the "pre-beetle kill" era or in other words, the endemic stage before the epidemic fully manifested itself. The 2009 image was defined as the "post-beetle kill" era or when the epidemic had progressed. In order to quantify the approximate amount of dead forests primarily due from beetle kill, NAIP aerial images from both 2005/2006 and 2009 were imported into ArcGIS and a maximum likelihood classification performed. ArcGIS outputs a classification raster for each aerial image and the difference in land cover change between the two periods showed that a 28% decrease in forested vegetation had occurred. ArcGIS then output a confidence raster based on the classification raster, which determined the approximate amount of uncertainty in this classification procedure. It was determined that an uncertainty of 12% was likely to occur in the classification procedures of the dead forests. From this classification analysis, a 95% confidence interval of 28% +/-12% was established.
Along with each output classification raster, the mean vector and the variance-covariance matrix were used to create an output confidence raster from which confidence intervals were constructed for each classification. Knowing the mean and confidence intervals, a population density function was created and subsequently used to vary the forest cover decrease in the VIC model, capturing the uncertainty in the classification results. In order to mimic beetle kill for the affected areas throughout the NPRB, a variety of simulations that decreased forested land cover over the entire basin were implemented to show if changes in streamflow would become apparent. A total of five different simulations were performed with the VIC to simulate decreases in land cover change due to beetle kill. Each individual simulation was decreased uniformly and not spatially over the entire basin, based on the results from the maximum likelihood classification. The following were the percent decreases in forested land cover change for each simulation in VIC: 1) 16%, 2) 22%, 3) 28%, 4) 34%, and 5) 40%, which are based on a 95% confidence interval.
In order to simulate land cover change by the respective amounts above, the vegetation parameter file that contains each of the 95 grid cells being modeled in VIC along with each land cover classification for each grid cell had to be modified. Each individual grid cell within the VIC model contains a certain number of land cover types. Referring back to Table 1, for the NPRB there are a total of 11 classifications that can reside within each grid cell. As stated previously, the NPRB contains 95 grid cells (1/8°), and it was concluded that 54 of those cells contained forested cover (Fig. 1). For this analysis, "forest coverage" was defined as only land cover classes #1-6 (refer back to Table 1).
For each of the 54 forested grid cells, if a grid cell contained land cover classification #1-6, then the value of the fractional coverage (Cv) for that particular land cover type was decreased by the specified amount for whichever of the five simulations were currently being performed. If a grid cell contained multiple types of forested vegetation, which occurred in most of the grid cells, then the percent decrease was divided by the total number of forest cover type classes #1-6 and then evenly distributed amongst the total number of forested land cover types within the grid cell. Once this was completed, then the total percent decrease that was applied for that particular scenario was then added to a less influential hydraulic land cover class, which for this analysis was assumed to be grasslands. In order to try to rationalize what was occurring in the NPRB, if a forest was initially in place and then was affected by beetle kill, in the process of modeling in VIC, it was concluded that the forest would basically transition from forest to grassland. In order to simulate decreasing forest cover in VIC, because all the total fractional coverage's for each land cover type in a grid cell had to sum to one, the above mentioned procedure was carried out for each of the 54 forested grid cells. In essence, for modeling purposes, it was assumed that a forest affected by beetle kill would become similar to that of grassland, thus enabling a greater potential for surface water runoff to occur and potentially increased streamflow. The above rationale was applied to each of the five previously mentioned simulations by use of VIC and modeled streamflow results produced.

Model Calibration and Validation
The VIC model was calibrated and validated with respect to the historical monthly observed streamflow for the period of 1950-1980 and 1981-2000, respectively [45].

Streamflow Response Due to Land Cover Change
By decreasing tree cover loss in the designated 54 forested grid cells in VIC by the specific amounts (16%, 22%, 28%, 34% and 40%) based on the results from the maximum likelihood classification performed in ArcGIS, annual modeled streamflow results were obtained for the period of 1950-2000. Uncertainty was captured by varying the forested vegetative cover from 16% to 40%. By understanding the effects from the beetle kill infestation in the past an estimate could be made as to how the effects of beetle kill could present itself in the near future.
The impacts of beetle kill on streamflow have been summarized in Table 2. The data presented here shows annual (1950 to 2000) streamflow for modeled (without "beetle kill") and for each of the five simulations (16%, 22%, 28%, 34% and 40%). For each of the five simulations, as tree cover continually decreased, streamflow continues to increase, which what was suspected.
Modeling efforts suggest that an increase in decadal modeled streamflow of approximately 1% to 10% for a decrease in land cover ranging from 16% to 40% will occur. The average change in decadal streamflow volume for the NPRB can be expected to be somewhere in the realm of a 5% increase for the average decrease in land cover of 28%. The t test and Rank Sum test were applied to compare each of the five "beetle kill" simulations to the modeled (no "beetle kill") to determine if there was a significant difference (increase) in streamflow due to "beetle kill" [46] [47]. For all five simulations, p values were large and, thus, none of the simulations showed a significant increase in streamflow (

Uncertainty
The authors recognize the current research and application of the VIC model was a "macroscopic" approach. The variation in tree loss (12% to 40%) and the resulting land cover change in VIC captures the uncertainty in this approach. However, possible sources of uncertainty in this analysis include: 1) VIC has been historically known to not perform as well in arid and semi-arid environments, such as the NPRB, because the model has problems with under or over-estimating peak flows for some years. However, for the VIC model used in this analysis was well calibrated and tested. 2) The NAIP aerial imagery data did initially present itself with uncertainty, but it was realized every dataset has some degree of associated uncertainty. However, because the NAIP imagery was thought to be the best available dataset for a variety of reasons, it was employed. 3) For each simulation, forest cover loss was uniformly decreased across all forested grid cells, but in reality land cover should have varied spatially to limit uncertainty. 4) It was decided that for modeling purposes, if a tree was "dead" then it would act the same as grassland, which was an assumption and has potential associated uncertainty. 5) For this analysis, when forested vegetative cover was [36] decreased, it was assumed to be either "alive" or "dead". According to Pugh and Gordon (2012) [5], this is a common misconception because a tree cannot be considered "alive" one second and then once infected by beetle kill, "dead" the next second. The death of a tree due to beetle kill is a gradual process that can take up to a year and a half depending on the tree species. 6) Lastly, when interpreting aggregated data such as this, it is important to be aware of the Modifiable Areal Unit Problem (MAUP) and its effect; depending on how the data is aggregated, be that different sized and/or positioning of the grid cells, the results may change to some degree. The MAUP is a potential source of bias that can affect spatial study results, which utilize aggregate data sources [47].
According to Pugh and Gordon (2012) [5], there are 6 stages for the life of a tree when associated with beetle kill. Initially, all healthy trees, prior to beetle infestation are considered in Stage 0: undisturbed forest, where natural hydrologic processes are functioning properly. Once beetle kill has commenced within the tree, the tree has entered into Stage 1: green phase, where the tree begins to slowly die due to lack of nutrients. At this stage, the tree retains its green needles but quickly begins reducing ET rates due to lack of absorption from the root system. This ultimately increases the amount of moisture remaining in the soil. Within about a year, the tree will enter into Stage 2: red phase, where all the needles have turned red and then eventually brown. At this point ET and other natural processes within the tree have ceased entirely and the tree can be considered completely dead.
Based on Pugh and Gordon (2012) [5] discussion of the phases that trees go through, this could potentially be a source of uncertainty for this analysis since we do not know how many trees are in what stage of beetle kill. Obviously, the trees, most of which are lodgepole pine, that were affected initially when beetle kill occurred in the NPRB around 2007/2008 are probably entering into Stage 3: grey phase in which they are essentially skeletons with only trunks and braches exposed. However, trees that have been recently, within the last year or two, affected could possibly still have some functioning ecohydrologic processes within them, thus the reason for our uncertainty. It is also important to realize that as the trees enter Stage 4: tree fall phase, where trees have died and begun falling to the ground, this can have a tremendous impact on the saplings and other smaller trees located in the understory of the forest. As the bigger trees die, many of the smaller trees in the understory will become more susceptible to more amounts of precipitation, due to the fact that the larger tree canopies are completely gone, and will commence the process of absorbing more water from the ground. This process that has just been described is the beginning steps of Stage 5: forest regeneration.

Conclusions
This research applied the VIC model and evaluated the impacts of beetle kill on streamflow within the NPRB. The impacts of the current (2009) land use change due to beetle kill infestation have been applied to past  historical meteorological forcing data to attempt to model beetle kill in the past and determine the results on streamflow. By performing this analysis, the potential effects of beetle kill on past streamflow could be applied to the present day conditions and possibly predict future water yields.