Assessing the long-term hydrologic response to wildfires in mountainous regions

This study aims to understand the long-term hydrologic responses to wildfires in mountainous regions at various spatial scales. The Soil and Water Assessment Tool (SWAT) was used to evaluate hydrologic response of the upper Cache la Poudre watershed in Colorado to the 2012 High Park and Hewlett wildfire events. A baseline SWAT model was established 10 to simulate the hydrology of the study area between the years 2000 and 2014. The effects of wildfires on land cover were accounted for in the model using the SWAT land use update module. The wildfire effects on curve numbers were determined comparing the probability distribution of curve numbers after calibrating the model for pre and post wildfire conditions. Daily calibration and testing of the model produced “very good” results. No-wildfire and wildfire scenarios were created and compared to quantify changes in average annual total runoff volume, water budgets, and full streamflow statistics at different 15 spatial scales. At the watershed scale, wildfire conditions showed little impact on the hydrologic responses. However, a runoff increase up to 75 percent was observed between the scenarios in sub-watersheds with high burn intensity. Generally, higher surface runoff and decreased subsurface flow were observed under post-wildfire conditions. Flow-duration curves developed for burned sub-basins using full streamflow statistics showed that less frequent streamflows become greater in magnitude. A strong (R > 0.8) and significant (p < 0.001) positive correlation was determined between runoff increase and 20 percentage of burned area upstream. This study showed that the effects of wildfires on hydrology of a watershed are scaledependent. Results also revealed that the wildfires had a higher effect on peak flows, which may increase the risk of flash floods in post-wildfire conditions.

Pages 2-3, Lines 30-32 and Lines 1-8: I'm not sure this content achieves the intended (assume to provide justification for the current study).The text here suggests the requisite approaches use static variables to represent dynamic properties.While this approach does have limitations, utilizing a dynamic approach that may not fully represent the dynamics at hand also has limitations.Anyway, the text doesn't necessary make a compelling case for one approach over another.Also, the word "components" need some reference/definition within this text.
That is a valid point that lacking the ability to capture the full dynamic of the problem has limitations.However, incorporating a dynamic approach even with its limitations can enhance the physical base of the simulation and help with more realistic representation of the processes compared to a static approach.We added more references and some more description of the method to better show the merits for its application.By component we mean the modules in the model that provide the capability to implement a specific change in model.Here, LULC component is the land use change module.The explanation is added to the manuscript.
Page 5, Lines 7-9: This is a very broad statement and assumes the continuous model is accurately depicting the dynamic events.Perhaps additional citations would better support this statement as the norm.I'm not entirely convinced the approach used in this study demonstrates a better representation or just a different one.Both approaches can be effective, useful, and yield good results.Any general commentary on one versus the other should be clearly justified and include substantial citations in support.
Both event-based and continuous-time models are useful and one can better represent processes compared to the other depending on the study problem.We did not intend to imply that continuous-time models are better than event-based models or they represent the processes accurately.For this study we had to use a continuous time model and in this paragraph we explain why a continuous time model better represents the processes involved in this specific problem.Some more references were added to the text for this argument.
Pages 5-6, Lines 24-30 and Lines 1-4: More specifics needed here and explanation of Range-Grasses approach is merited.
Thanks.The reason for updating the LULC to Range-Grasses was explained in more detail in the text.
Page 9, Line 2: I'm not sure I agree that the error statistics tell how accurately SWAT is representing processes exactly.
We used the error statistics to assess how good the model simulations conform to real observations.While the good error statistics does not guaranty that the hydrologic processes are represented completely accurate, it gives a good insight into the performance of the model.The text was reworded to avoid confusion.
Page 12, Line 13: The inference that it may be reasonable to use total burn area percentage as a predictor requires some qualification here given the single study.
The idea that increase in percent burn area results in increase in runoff has been reported in several studies (Benavides-Solorio and MacDonald, 2001;Inbar et al., 1998;Lavabre et al., 1993;Robichaud et al., 2000;Scott, 1993).However, we were not able to identify studies where the relationship between runoff increase and percent total burned area was quantified similar to our study.Given the regression analysis done in this study and checking the metrics indicating the strength and significance of the regression model, we expressed that percentage of total burned area may be useful as an indicator of increase in runoff.Anyway to be more accurate we added "in the Cache la Poudre watershed" since as you have mentioned this is a single study on a single watershed.We also mentioned that more studies are needed to generalize the findings of the study.

Introduction
Assessing the hydrologic and water quality effects of wildfires is becoming more important as the frequency and severity of wildfires in the United States (U.S.) and other regions around the world have shown increasing trends (Westerling, 2016;Doerr and Santin, 2016;Ebel et al., 2012).Wildfires may have undesired consequences for water quality, carbon storage, and ecosystem disturbance (Gould et al, 2016;Holden et al., 2012;Moody and Martin, 2009).Wildfire prone regions, often with increasing populations, are susceptible to loss of life and catastrophic destruction from floods and debris flows as a result of higher runoff and erosion under post-wildfire conditions (Moody et al., 2013).Assessing the hydrologic effects of wildfires in mountainous regions is particularly interesting as these areas often contain the headwater sub-watersheds supplying water that is relied upon downstream (Viviroli et al., 2007).Wildfires could alter the timing and magnitude of runoff, subsurface flow, along with other hydrologic fluxes which eventually will reduce the reliability of water supply in these areas (Ebel et al., 2012).
Characterization of complex responses to wildfires is difficult due to the spatial variability of post-wildfire conditions (Moody et al., 2013).Wildfires can substantially change land use-land cover (LULC) and vegetation within watersheds, which may subsequently result in altering hydrologic regimes including: (1) increased availability of rainfall for runoff by decreasing canopy interception (Moody and Martin, 2009;Robichaud et al., 2000), (2) increased base flow through the decrease of water normally lost through evapotranspiration (Neary et al., 2003), and (3) increased runoff velocities and reduced interception/storage through loss of ground cover, litter, duff, and debris (Moody and Martin, 2001).These alterations can cause increased hillslope erosion and may significantly alter terrestrial habitat.They may also increase channel flooding, decrease channel stability, fill the streambed with fine sediment, and modify temperature regimes (Ryan et al., 2011).
Mathematical modelingmodelling is a useful and well accepted approach for improving our understanding of complex watershed processes (Kiesel et al., 2013).For example, watershed models have been used for simulating streamflow in mountainous regions to identify important hydrologic interactions and processes (Sanadhya et al., 2014).Various approaches are adopted in the literature for modelling the hydrologic and water quality processes in watersheds ranging from purely empirical to fully process-based models (Beven, 2001;Famiglietti and Wood, 1995).Process-based hydrologic models have been used in the literature extensively to simulate the effects of natural and anthropogenic perturbations on the hydrologic regime (Clark et al., 2017).Application of these models to assess the effects of wildfires entails availability of mechanisms for updating model parameters and inputs such as land use which may not be readily available in most models.The Soil and Water Assessment Tool (SWAT) is a process-based distributed parameter watershed model that has been used to characterize and quantify the effects of LULC change, climate change, and mitigation strategies on runoff, evapotranspiration, streamflow, groundwater and other hydrologic responses showing very good results in terms of model performance (Tasdighi et al., 2017;Motallebi et al., 2017;El-Khoury et al., 2015;Fan and Shibata, 2015;Foy et al., 2015).
More specifically, numerous studies involving SWAT model development and calibration have been conducted to evaluate the hydrology in mountainous and snow-driven regions throughout the world, including this study watershed (Foy et al., 2015;Sanadhya et al., 2014); Cannonsville Reservoir watershed, New York (Tolson and Shoemaker, 2007); the Little River watershed, Tennessee (Zhu and Li, 2014); two Himalayan drainages of Nepal (Neupane et al., 2015); and the Yingluoxia watershed of northwest China (Lu et al., 2015).These studies document promising results for application of SWAT for hydrologic simulations in mountainous regions.
One of the challenges in using models for evaluating the hydrologic response of a system to wildfire is developing the mechanism through which the hydrologic effects of wildfires are simulated.Studies have used alteration of model parameters and LULC to represent pre and post wildfire conditions (Parson et al., 2010;Robichaud, 2000).The majority of these studies have used field measurements or implemented a deterministic approach using fixed values for specific parameters during the pre and post wildfire conditions to represent the change in hydrology as a result of wildfire (Parson et al., 2010;Robichaud, 2007).This approach has a direct impact on the results obtained as selection of parameters to change and the values assigned to them during pre and post wildfire conditions are subjective and may not necessarily represent the changes in the real world.Additionally, lack of components required components (modules in the model) for representing LULC change adds to the complexity of the procedure (Batelis and Nalbantis, 2014;Goodrich et al., 2005;McLin et al., 2001).A proper LULC change component in the model is required for continuous simulation and is particularly important when assessing effects of wildfires.The LULC change module within SWAT has been shown useful for evaluating hydrologic condition where LULC has changed as the result of urbanization (Pai and Saraswat, 2001).
The 2012 Hewlett and High Park wildfires have provided a unique opportunity for examining hydrologic response to wildfires, specifically, in a mountainous region.This unique opportunity stems from the fact that a relatively significant proportion of the gaged Cache la Poudre (Poudre) headwaters (approximately 14 percent from the Mouth of Canyon) has been burned as the result of wildfire.The pre and post-wildfire streamflow data availability allows for the development, calibration, and testing of a hydrologic model that accounts for spatial variability in LULC to continuously simulate the hydrology from pre-wildfire conditions through post-wildfire conditions.Due to the magnitude of the 2012 wildfire incident, burn severity mapping is available for the area.This mapping data allows for a land use change module to be implemented during calibration efforts which adjusts hydrologic parameters impacted by wildfire seamlessly during simulation.
The overall goal of this study is to characterize and quantify the long-term (a two-year period after the wildfire) hydrologic responses to wildfires in mountainous regions at various spatial scales (smaller high burn intensity sub-watersheds to watershed scale).To accomplish this goal, the Poudre headwaters located in northern Colorado, USA which experienced the 2012 Hewlett and High Park wildfires was analyzed with the SWAT model.This analysis includes simulation of no-wildfire and wildfire scenarios over a 15 year (2000 to 2014) period.Specific objectives of this study are to: (1) quantify changes in average annual total runoff volume and explore how these changes fluctuate with the percent of the area burned, (2) quantify annual changes in various components of hydrologic budget, and (3) highlight potential implications of these changes using full streamflow statistics through application of flow duration curves at both sub-watershed and watershed scales.While a number of previous studies have examined the hydrologic effects of wildfires, application of a probabilistic approach for characterizing the change in key hydrologic parameters between the pre and post wildfire scenarios and a dynamic LULC updating through the analysis period is novel.The results of this study have important implications for hydrologic effects of wildfires and methods used to assess them.

Study area
The Poudre Watershed, with an area of approximately 5,230 km2 above its confluence with the South Platte River on the Great Plains, is situated mostly in northern Colorado, USA with a portion reaching into southern Wyoming, USA (Wohl, 2010).The Poudre River (Figure 1) is supplied by two major tributaries within its headwaters, the South and North Forks, the latter being the longer of the two joining the main-stem farther downstream.After streamflow retreats from the Poudre's headwaters in the Rocky Mountain Range, the river passes through the cities of Fort Collins and Greeley.Eventually, the river joins the South Platte River and winds downstream to join the Platte River and then to the Missouri River.The Poudre River, with its minimally-developed mountainous headwaters, is widely utilized as a source of drinking water for several cities and communities located along its banks (Richer, 2009).
During May and June of 2012 the Hewlett and High Park wildfires burned approximately 384 km2 of primarily forested landscape within the Poudre Watershed.The burned area includes numerous drainage tributaries to the main-stem of Poudre River.Widespread loss of vegetation and burned soils from the wildfires created areas susceptible to severe erosion and flooding.Localized summertime thunderstorms immediately following the wildfire worsened the effects by washing sediment and debris into the river channel posing a threat to the safety of people and homes in the area (Oropeza and Heath, 2013).The affected area extends along the Poudre River from the mountain front upstream to several kilometers south of the community of Rustic, Colorado.Therefore, a study watershed outlet was defined near the mountain front at Colorado Division of Water Recourses' (CDWR) surface water gauge CLAFTCCO18 (formally USGS Gage 06752000), commonly referred to as the Mouth of Canyon (Figure 1).
The resulting study watershed is approximately 2,732 km2.At higher elevations, streamflow is dominated by snowmelt runoff and at lower elevations rainfall runoff from summer convective storms greatly affect streamflow.The storms combined with the upstream snowmelt runoff, can produce high-magnitude, short-lived floods at times (Wohl, 2010).The resulting hydrograph is snowmelt dominated with a rise typically beginning in April and a recession lasting into August.
Generally, peak streamflow occurs at the end of May or early June and base flow levels occur in September or October (Richer, 2009).

SWAT model
SWAT is a continuous-time, distributed-parameter, process-based watershed model, which has been used extensively for hydrologic and water quality assessments under varying climatic, land use, and management conditions in small watersheds to large river basins (Gassman et al., 2007;Neitsch et al., 2011;Arnold et al., 2012, CARD Staff, 2016).
SWAT model allows for numerous physical processes to be simulated in a watershed.These processes may be separated into two coarse divisions of the hydrologic cycle: the land phase and the routing phase.The main processes include precipitation, infiltration, surface runoff, evapotranspiration, groundwater flow, snowmelt, and flood routing.SWAT is driven by a water balance equation which relates individual components of the hydrologic cycle.Additional details including specific equations associated with the water balance and the individual hydrologic processes may be found in the SWAT Theoretical Documentation, Version 2009 (Neitsch et al., 2011).
In order to assess the hydrologic effects of wildfires, application of a continuous time model is necessary as it will provide the capability to explore the long-term hydrologic effects (Jeong et al., 2010).Compared to event-based models, continuous time models better represent watersheds where channel storage may be significant and/or where significant variability exists in land use (e.g., urbanization), soil types, and/or topography (Arnold et al., 1995;Arnold et al., 1998;Nicklow et al., 2006;Jeong et al., 2010).Being a distributed-parameter model SWAT divides a watershed into sub-basins, which are further divided into hydrologic response units (HRUs).HRUs are the smallest spatial units in SWAT, and are defined as areas within each subwatershed with unique combinations of land use, soil, and slope class.Sub-basins can be assigned unique climate and hydrologic properties which in combination with unique land use characteristics of HRUs provides the capability to investigate the effects of land use change scenarios under varying climatic conditions both spatially and temporally.

Updating LULC
The hydrologic modelingmodelling process was initiated by first collecting and preparing the necessary data, summarized in Table 1.A detailed description of each data type is presented in Appendix A.1.The NLCD 2011 Land Cover spatial dataset was preprocessed to allow the High Park and Hewlett wildfires to be simulated by SWAT.The NLCD 2011 Land Cover was overlaid with the Thematic Burn Severity dataset (Monitoring Trends in Burn Severity Project, 2014).Then the NLCD 2011 Land Cover was reclassified to incorporate low, medium, and high burn severity categories.The reclassification was accomplished using spatial analyst toolset within ArcMap from ESRI's ArcGIS software package.The preprocessing retained the pre-wildfire classification, but added a burn severity identifier.For example, portions of the NLCD 2011 Land Cover that consist of Evergreen Forest and overlap with a low burn area were reclassified to a newly created Evergreen Forest Low Burn classification.
The SWAT Model Database contains various pre-defined model parameters for different LULC types and the SWAT LULC lookup table relates NLCD classifications to the LULC types found in the SWAT Model Database.In order to seamlessly represent wildfire during the simulation the SWAT Model Database and SWAT LULC lookup table were altered to reflect the conditions after the wildfire.
For the pre-wildfire database, the newly added LULC types consisted of attributes identical to the original classification, but with a new description and identification code.Thus, the SWAT model created using this database will represent prewildfire condition, but areas influenced by wildfire will be delineated from non-burned areas.
For the post-wildfire database, the newly added LULC types included a new description and identification code similar to the pre-wildfire database; however, for the post-wildfire case, attributes were also altered from their original classifications.
Wildfires result in loss of canopy.Most of the canopy in the study watershed was evergreen forest (NLCD 2011).This canopy is mainly replaced by shrubs and range grasses after the wildfire.In order to mimic this alteration in canopy, fFor all burned areas, LULC attributes in the database were changed to match those of the Range-Grasses LULC.This land use change results in change in several land use-specific predefined parameters available in model databases.This change was implemented to aid with appropriately representing loss ofof the alteration in canopy in burned areas.

Updating Curve Numbers
Curve Numbers (CNs) were adjusted to account for expected increases in runoff.The change in CNs was based on a pre and post wildfire calibration of the model explained in section 2.2.7.Comparing the probability distribution of CNs before and after the wildfire, an average CN increase of 5, 10, and 15 for areas with low, moderate, and high burn intensity were considered respectively.The original and edited SWAT LULC lookup tables as well as curve numbers for both pre and postwildfire conditions can be found in Appendix A.2, tables A3 through A5.

Initial model development
Two models representing pre and post-wildfire conditions were developed and then later merged to create a unified model.
Two sets of initial SWAT model input files for the study watershed were created using ArcSWAT 2012 (U.S.Department of Agriculture Agricultural Research Service, 2014).ArcSWAT is an ArcMap extension that provides a graphical user interface for creating a SWAT model.The interface was used to process the previously described model data to generate initial SWAT input files.This process is summarized in Figure 2.
The ArcSWAT Automatic Watershed Delineation tool was used to create a stream network, define sub-basin outlet locations, delineate the watershed, and calculate the sub-basin parameters.Additional outlets were manually placed at locations where a large tributary entered the study reach and the whole watershed outlet was defined at the Mouth of Canyon.ArcSWAT was then used to determine LULC/soil/slope combinations within each sub-basin which is then used to determine the distribution of HRUs for the entire watershed.

Model options
Options for both models are identical and were selected based on previous modelingmodelling studies using SWAT in mountainous regions (Foy et al., 2015;Lu et al., 2015;Neupane et al., 2015).A modified version of the commonly applied United States Soil Conservation Service (now the NRCS) curve number (CN) procedure was adopted to simulate surface runoff in the watershed.The CN depends on the soil type, LULC, and hydrologic condition (Lu et al., 2015).Penman-Monteith method based on energy balance components was selected to estimate potential evapotranspiration.Lastly, channel routing was represented using the Muskingum River Routing Method.Other model options were left as default configurations.

Accounting for Mountainous Terrain
The study watershed is located within the rainshadow of the Rocky Mountains and overall experiences a topographicallydriven climate.Significant difference in elevation within the study watershed yields large variability in the quantity and form of precipitation.Thus, lapse rates as well as elevation band parameters were assigned to each sub-basin to account for orographic effects.The precipitation lapse rate (i.e., increase in mean annual precipitation with an increase in elevation) of 658.4 mm/km obtained from Foy (2015) was incorporated into the model.Additionally, the temperature lapse rate (i.e., decrease in mean annual temperature with an increase in elevation) of -5.5 ⁰C/km reported by Foy (2015) was used.
SWAT is capable of integrating up to 10 elevation bands in each sub-basin.These bands were derived by topographically discretizing each sub-basin within the watershed.SWAT requires the input of the elevation at the center of each band and the fraction of sub-basin area within the elevation band.Data from the Topography Report generated by AcrSWAT was used to discretize each sub-basin into 10 elevation bands.The minimum elevation was subtracted from the maximum elevation and divided by ten, which creates ten equal-interval elevation bands.Next, the elevation at the center of each band and the fraction of sub-basin area within the elevation band is calculated.Lastly, the previously generated SWAT input files were modified to contain these parameters.These parameters allow SWAT to use the elevation band equations described in the SWAT Theoretical Documentation, Version 2009 (Neitsch et al., 2011) to simulate orographic effects.
The curve numbers provided in the SWAT Model Database are appropriate for slopes up to 5 percent (Neitsch et al., 2011).
An analysis of the HRUs revealed that many of the average HRU slopes in the study watershed exceed 5 percent.We adjusted the curve numbers for different slopes at the HRU level using the following equation (Neitsch et al., 2011) where CN 2S is the moisture condition II CN adjusted for slope, CN 3 is the moisture condition III CN for the default 5 percent slope, CN2 is the moisture condition II CN for the default 5 percent slope, and S is the average fraction slope of the subbasin.Note that upon simulation SWAT caps CN values at 98.

Land use update module
The pre and post wildfire models were used to create a unified model by activating the land use update module in SWAT.A

Model calibration and testing
Calibration of process-based distributed hydrologic models is often challenging due to high data demand to support the complexities (Beven, 2001).Many studies use single source of data (e.g.streamflow or nutrient concentrations) to calibrate such models (Ahmadi et al., 2014;Foy et al., 2015).While this is a common approach in calibrating the models due to difficulties in data acquisition and limited resources, it also has limitations.Application of several data sources for calibration of process-based models can improve the performance of the model by assuring that various modules within the model have a good performance.The SWAT model was calibrated and tested for the daily naturalized streamflows at the Mouth of Canyon.The naturalized flows were determined using the daily historical flow records of diversions or reservoirs.
Separate calibrations were conducted for pre and post wildfire conditions with the purpose of characterizing the change in CN as a result of wildfire.Once the change in CN was characterized, it was implemented in the unified model.The unified model was then calibrated and tested for the whole period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014).Calibration, pre-wildfire testing, and post-wildfire testing periods were 2005-2013, 2000-2004, and 2014, respectively.These simulation periods were selected based on data availability.Initial calibration parameters were identified from previous modelingmodelling and sensitivity analysis efforts for the study watershed (Sanadhya et al., 2014;Foy et al., 2015).They used eFAST (Saltelli et al., 1999) which is a variancebased global sensitivity analysis method to determine the sensitivity of streamflow simulations to SWAT parameters.Both first and total order sensitivity indices were determined and used to rank the parameters.These parameters were supplemented with additional parameters identified from a previous sensitivity analysis study utilizing using the method of Sobol with the SWAT model (Ahmadi et al., 2014).A total of 38 modal parameters were used for calibration (Table A6).
Parameters related to subsurface processes, specifically the soil hydraulic conductivity (SOL_K) and tributary channel hydraulic conductivity (CH_KI) had the highest effects on simulations in the watershed (Sanadhya et al., 2014).A Matlab code was developed for auto-calibration of SWAT model using a global optimization algorithm named dynamically dimensioned search (DDS; Tolson and Shoemaker, 2007).DDS is designed to arrive at good solutions within a maximum number of user-defined function evaluations for use in model calibration with many parameters (Tolson and Shoemaker, 2007).This auto-calibration tool was used to generate 498 model runs.Each model run consisted of a unique combination of the 38 model calibration parameters.The tool works towards minimizing an objective function.In this case, we based this objective function on two primary error statistics, relative error (RE) and Nash-Sutcliffe efficiency coefficient (E NS ).E NS is a normalized statics that indicates how well observed versus simulated plot fits a 1:1 line.E NS is computed as: where for this study Y obs is the observed streamflow, Y sim is the simulated streamflow, and Y mean is the mean of observed streamflows (Nash and Sutcliffe, 1970).The optimal value for E NS is 1.E NS can range between -∞ and 1.0, with values between 0.0 and 1.0 generally regarded as satisfactory levels of performance.Values equal to or smaller than 0.0 indicate that the simulated values are a same or worse predictor than the mean of observed values (Moriasi et al., 2007) respectively.
The RE gives an indication of how good simulated values are relative to the magnitude of corresponding observed values.
RE in percentage is computed as: where for this study Y obs is the observed streamflow and Y sim is the simulated streamflow.These error statistics are used to determine how well the model simulations match the observations accurately SWAT is representing hydrologic processes through comparison of observed and simulated streamflows at the Mouth of Canyon.While they may not necessarily mean that all hydrologic processes are accurately represented, they are a good indicator of the performance of the model with regard to fitting the observations.Model calibration parameter starting values and ranges are displayed in Appendix A2, table A6.

Scenario analysis
With the SWAT model calibrated and tested, two scenario models were created.First, a no-wildfire scenario model was created.This was achieved by simply removing the land use update files thus representing no wildfire activity throughout the entire simulation period.Second, a wildfire scenario model was created.This was achieved by adjusting the land use update files to reflect a wildfire occurring at the beginning of the simulation.Thus, wildfire is simulated throughout the entire simulation period.Note the simulation period for each scenario was between 2000 and 2014 (15 years).

Output Data Post-Processing
SWAT outputs were post-processed in Matlab.Simple summing functions were used to calculate total runoff volumes and water budgets throughout the study watershed.Full streamflow statistics were used to develop flow-duration curves for burned sub-basins.These represent the percentage of time that streamflow is likely to equal or exceed a given streamflow value for both scenarios.The code used sorts, ranks, and plots the input streamflow data to generate flow-duration curves.
Flow-duration curves are a widely accepted method for characterizing streamflow regime.They are commonly used for hydropower, water resource management, water quality management, habitat suitability, and flood control applications (Fan and Li, 2004).However, they have not been frequently used in evaluating response to wildfires (Newtson, 2013).Next, the ecodeficit and ecosurplus metrics introduced by Vogel (2007) were computed for each flow-duration curve.These metrics provide a simplified representation of hydrologic impacts (Vogel et al., 2007).For this study, ecodeficit is defined as the ratio of the area below the no-wildfire scenario flow-duration curve and above the wildfire scenario flow-duration curve divided by the total area under the no-wildfire scenario flow-duration curve.Conversely, ecosurplus is defined as the ratio of the area above the no-wildfire scenario flow-duration curve and below the wildfire scenario flow-duration curve divided by the total area under the no-wildfire scenario flow-duration.Thus, these values represent the overall loss (ecodeficit) and gain (ecosurplus) in streamflow (Vogel et al., 2007) between scenarios.

Model performance
The performance of the model during calibration and testing at various temporal scales was assessed using the common criteria in the literature (Motovilov et al., 1999;Moriasi et al., 2007).Model performance at a daily timestep was considered good if E NS ≥ 0.75 and was considered satisfactory for values of E NS between 0.75 and 0.36 (Motovilov et al., 1999).At monthly timestep, the performance of the was categorized as very good (0.75 < E NS ≤ 1.00), good (0.65 < E NS ≤ 0.75), satisfactory (0.5 < E NS ≤ 0.65), and unsatisfactory (E NS ≤ 0.5) (Moriasi et al., 2007).
During the separate pre and post wildfire calibration the model had a good performance giving E NS = 0.92; RE = -0.32%,and E NS = 0.81; RE = 2.16% respectively.This separate calibration step was done to characterize the change in CN between pre and post wildfire conditions.Figure 3 illustrates boxplots of CNs for the pre and post wildfire conditions.The values of CN for the best solution (highest E NS ) and mean of the CNs are also marked on the boxplots.Based on these results, an average CN increase of 5, 10, and 15 for areas with low, moderate, and high burn intensity were considered respectively.
The optimal parameter set found during the calibration effort generally yielded good results.The model performed best during the post-wildfire testing period, but still performed well during the calibration period and pre-wildfire testing period.
Final values for the 38 calibration parameters are displayed in Appendix A2, table A6.Model performance was evaluated based on primary statistical results (at both the daily and monthly timesteps) and visual inspection of the graphical results.
The best calibration achieved for the Mouth of Canyon naturalized streamflow at the daily timestep was E NS of 0.82 and RE of 1.68.The testing E NS values for the pre-wildfire and post-wildfire periods were 0.71 and 0.88, with RE values of -19.52% and 9.31%, respectively.Table 2 presents a summary of the model performance at the daily timestep.
All simulation periods earned a performance rating of very good at the monthly timestep.Monthly results were generally comparable to those from other SWAT modelingmodelling studies involving mountainous watersheds (Foy et al., 2015;Lu et al., 2015;Neupane et al., 2015).Table 3 presents a summary of the model performance at the monthly timestep.
Generally, simulations yielded good visual agreement between observed and simulated daily streamflows and total runoff volume, as shown in Figure 4.A slight discrepancy between the observed and simulated total runoff volume exists for the no-wildfire testing period.This difference propagates to the statistical results, most notably, the RE value of -19.52%.A negative relative error shows that the model overestimates runoff volume compared to observations.Based on visual examination of the hydrographs, the calibration period may be slightly "wetter" relative to the pre-wildfire testing period, which may be the cause of the noted discrepancy.
Also, the simulated and observed flow-duration curves for the entire simulation period yielded good visual agreement, as shown in Figure 5.The simulated flow-duration curve generally follows the observed flow-duration curve with the exception of a slight deviation for less frequent flows.For the less frequent streamflows the model is underestimating streamflows.A deviation is expected as less frequent streamflows correspond to larger streamflows which are less predictable and less understood.
Previous studies have used SWAT along with similar calibration techniques throughout this region for hydrologic analysis.
However, use of the SWAT land use change module to investigate hydrologic response to wildfires has not been well documented.Moreover, characterization of change in CNs as a result of wildfires using a probabilistic approach has not been performed previously.The performance results above indicate that the comprehensive methodology of using the SWAT land use change module along with multi-variable parameter calibration was an effective technique to represent the hydrology of an area which has been exposed to wildfire

Model performance
The daily simulation outputs from both no-wildfire and wildfire scenarios were analyzed and compared in order to characterize an average hydrologic response to wildfire during the simulation period of 15 years (2000 to 2014).Total runoff values, represented as both depth and volume for each burned sub-basin as well as for the entire study watershed are shown in Table 4. Also, Figure 6 displays the burn severity distribution and average annual total runoff percent increase (based on the values presented in Table 4) for each burned sub-basin and for the entire study watershed.The average annual total runoff includes surface runoff, lateral flow, and base flow.
Figure 6 shows that in the case of sub-basins 28, 30, 26, and 32, more than 50 percent of the area experienced burning as a result of the High Park and Hewlett wildfires.Sub-basins 28 and 30 were the most severely burned with large high burn severity percentages.The remaining sub-basins had smaller burned area percentages.
The total runoff percent increase between scenarios was greatest on average for sub-basins 28, 30, 26, and 32.For these subbasins, increases in runoff between the no-wildfire and wildfire scenarios ranged from approximately 66 to 75 percent.For the remaining sub-basins, as well as the entire study watershed, runoff percent increases are found to be considerably less.
This is likely because those sub-basins were not as heavily burned.Nevertheless, the results indicate wildfire effects at larger scales are still substantial, but only in terms of the magnitude rather than percent change of total runoff volume increase.
Larger areas (i.e., sub-basin 35 and the entire study watershed) appear to experience much greater absolute increases in total runoff volume between scenarios, despite having smaller total burn area percentages.This is what we might expect given that each sub-basin is nested within the study watershed, resulting in a cumulative effect.
Other studies have documented total runoff increases under post-wildfire conditions (Benavides-Solorio and MacDonald, 2001; Inbar et al., 1998;Lavabre et al., 1993;Robichaud et al., 2000;Scott, 1993).For example, Lavabre et al. (1993) used a lumped conceptual hydrological model to evaluate a small Mediterranean basin which experienced a burn covering 85 percent of its surface area in 1990.They suggested a 30 percent increase in the annual runoff yield.Scott (1993) showed total streamflow volume increases of 15.3 and 9.4 percent in response to burning in two small mountainous catchments using a paired catchment method.In contrast, Mahat et al (2015) reported no significant change between the modeled streamflow from burned and unburned models.They suggested that this outcome may be the result of using a conceptual modelingmodelling approach instead of using a physically based model.The amount of total runoff volume increase following wildfire disturbance varies greatly between locations depending on wildfire intensity, proportion of the forest vegetation burned, climate, precipitation, geology, soils, watershed aspect, and tree species (Neary et al., 2003).Thus, it is not surprising that results vary.Also, comparison between studies is difficult because of changes in size of disturbance (i.e., wildfire) in relation to the size of the catchment (Robichaud et al., 2000).This emphasizes the need to examine increases based on percent burn area upstream.
Figure 6 is arranged in descending order of percent burned area from left to right.Generally, we see an increase in total runoff as percentage of total burn area increases.This observation is consistent with reports in the literature indicating total runoff volume increase following wildfire disturbance is in part a function of the proportion of the contributing area burned (Neary et al., 2003;Robichaud et al., 2000).This relationship is further explored by applying linear regression to the data.
Figure 7 shows a linear regression model fitted between the total runoff volume increase and total burned area percentage in the Cache la Poudre Watershed.Note that the entire study watershed results were not included in this regression.Also, subbasin average slope was categorized as low (slope < 0.30), moderate (0.30 ≤ slope < 0.40), and steep (slope ≥ 0.40) for each sub-basin.
An F-test was performed using Matlab to determine if this particular model fits the data well.The regression generally yields a good fit, with a p-value < 0.001 for the F-test.No previous study was found documenting this relationship with linear regression.Thus, tThis study suggests it may be reasonable to use total burn area percentage as a predictor for increase in total runoff volume in Cache la Poudre Watershed.Also, the figure indicates that generally for the High Park and Hewlett wildfires the sub-basins with moderate to steep slopes experienced wildfire in a larger percentage of their area relative to low slope sub-basins.It should be noted that the linear regression model was developed for Cache la Poudre Watershed and it may not necessarily be applicable to other watersheds.Further studies are required to draw a potential generalized relationship between percent burned areas and increase in runoff.

Wildfire effects on hydrologic budgets
The daily simulation outputs from both no-wildfire and wildfire scenarios were further analyzed and compared in order to quantify changes in average annual hydrologic budgets as a result of wildfire during the simulation period of 15 years (2000 to 2014).Figure 8 shows hydrologic budgets for select sub-basins as well as the entire study watershed.These hydrologic budgets show the fate of average annual precipitation along with the fate of average annual total runoff.The fate of precipitation (rainfall and snowfall) is shown as evapotranspiration, total runoff, and other (deep aquifer contribution and soil water storage).Also, the major hydrologic processes for the fate of runoff were defined as surface and subsurface (lateral flow and base flow) runoff.
It is evident that hydrologic budgets change on the sub-basin scale following wildfire; however, little change is seen at the watershed scale.Batelis and Nalbantis (2014) also documented that wildfire effects are practically indiscernible on a regional scale.Generally, Figure 8 shows under the wildfire scenario an increase in surface runoff and a corresponding decrease in subsurface flow at the sub-basin scale.For example, the hydrologic budget for sub-basin 30 (a heavily burned area) shows a change in surface runoff from 21 to 61 percent under the no-wildfire and wildfire scenarios, respectively.This is consistent with previous studies in which it seems to be generally accepted that infiltration rates decrease after wildfires.
For example, Moody and Martin (2001) showed that infiltration rates were decreased by a factor of two to seven after wildfires.
At the sub-basin scale under the wildfire scenario we also see less evapotranspiration.This connects well with the results from Section 3.2, where generally we see an increase in total runoff for the wildfire scenario.Increased water yields (i.e., total runoff) primarily due to reduced evapotranspiration has been a reported effect on post-wildfire hydrology (Neary et al., 2003;Townsend and Douglas, 2004).

Implications of wildfire effects
Lastly, the daily simulation outputs from both no-wildfire and wildfire scenarios were analyzed and compared in order to determine potential implications of wildfire effects during the simulation period of 15 years (2000 to 2014).Figure 9 shows flow-duration curves for select burned sub-basins as well as for the entire study watershed and Table 5 lists the ecosurplus and ecodeficit values associated with each computed flow-duration curve.Flow-duration curves were generated using total runoff, which includes both surface and subsurface water fluxes leaving the sub-basin or watershed.The ecosurplus and ecodeficit metrics are a dimensionless measure which represent the overall loss (ecodeficit) and gain (ecosurplus) in streamflow (Vogel et al., 2007) between scenarios.
Similar to findings from the hydrologic budgets, it is evident that flow-duration curves change under wildfire conditions on the sub-basin scale.Also, little change is seen at the watershed scale (Figure 9 and Table 5).This is perhaps the result of wildfire effects at the watershed scale being damped by non-burned portions of the contributing area.
Figure 9 also suggests that wildfire has little impact on flow-duration curves for areas with low total burn area percentages, but seems to impact flow-duration curves for area with higher total burn area percentages.For example, in sub-basins 30 we see that less frequent streamflows become greater in magnitude under the wildfire scenario (i.e.we see an ecosurplus).
Whereas, in sub-basin 19 (a less burned area) we see little change in the flow-duration curve.Previous research efforts have involved a paired-catchment analysis to compare flow duration curves for pre and post-wildfire conditions (Liu et al., 2004;Newtson, 2013).Both Newtson (2013) and Liu et al. (2004) found a general increase in percentile streamflow as a result of wildfire.However, Liu et al. (2004) examined precipitation duration curves for the study areas and concluded that changes in precipitation between locations explained the difference in streamflow and not necessarily wildfire.For this study, the two scenarios approach uses an identical precipitation record for both scenarios.Thus, the study eliminates limitations associated with temporal and special variation in precipitation.Table 5 indicates the streamflows for the burned sub-basins appear to be ecosurplus versus ecodeficit when the wildfire scenario is compared with the no-wildfire scenario.The ecosurplus values range from 0.004 to 0.279.Kannan and Jeong (2011) indicate that for high streamflows a large ecosurplus is likely to have moderate to high impacts to stream health.In this case, the ecosurplus values associated with the heavily burned sub-basins (i.e., sub-basins 28, 30, 26, and 32) are much greater in magnitude when compared to the other ecosuplus values.Thus, impacts to stream health are expected to be the greatest in heavily burned areas.specifically, those during the month of June when streamflows are elevated due to mountain snowpack melting.Also, the model appears to slightly underestimate streamflows during late summer into autumn.These systematic errors may be due to SWAT releasing snowmelt too quickly during spring runoff, thus, rising streamflows are simulated earlier than observations during the melting season.Further, perhaps the tendency of the model to simulate earlier snowmelt results in higher simulated streamflow during the latter part of summer and early autumn.This deficiency may be the result of SWAT misrepresenting snowmelt processes or perhaps faulty model parameterization.Thus, it is thought that hydrologic model uncertainty is introduced here and it is recommended that additional research be focused on better representing snowmelt processes in mountainous watersheds.

Conclusions
Long term simulation scenario analysis at the sub-basin and watershed scales was used to characterize hydrologic response to wildfires in mountainous regions.This was achieved by applying the hydrologic model SWAT simulation period were created from the optimal parameter set achieved during model calibration.These scenarios were used to characterize the hydrologic response to wildfires.
Specific objectives of this study were to investigate changes in average annual total runoff volume, average annual hydrologic budgets, and flow-duration curves across multiple scales as a result of wildfire.At the watershed scale, wildfire conditions appear to have little effect on the hydrologic responses with the exception of total runoff volume.However, at the sub-basin scale, simulations suggest that wildfire effects trend with burned area upstream.A total runoff increase up to approximately 75 percent between scenarios was found.Generally, water budgets showed more surface runoff versus subsurface flow, which suggests infiltration rates decrease under post-wildfire conditions.Flow-duration curves for burned sub-basins showed that less frequent streamflows become greater in magnitude leading to ecosurplus values of up to 0.279.
Results reported in this study show an overall acceptable performance of the SWAT model in simulating daily streamflows under pre and post-wildfire conditions to characterize the hydrologic response to wildfires.However, this method required comprehensive knowledge of the watershed, was time consuming, and was computationally intensive.Further, this study demonstrates the need for improvement in understanding the rainfall-runoff prediction relationship for burned areas.comprehensive LULC change analysis for the study watershed using NLCD 2000, 2006, and 2011 Land Cover is    is accomplished using a multi-tiered approach including a formatting check as well as a quality test looking for a variety of data problems.Based on this, no further quality control beside removal of flagged data was conducted.The stations were selected based on location, type of data provided, length of record, and completeness of record.A complete list of stations may be found in Table A2.Mean annual precipitation ranges from 330 mm at the lower elevations to 1350 mm at the higher elevations and mean annual temperature ranges from approximately 9° C at the lower elevations to -5° C at the higher elevations.Precipitation within the study watershed is greatest during the winter months.Snow accumulates which generates the mountain snowpack that is then released during the spring and early summer months.In an effort to support economic, environmental, and recreational water demands downstream, manmade structures such as diversions, storage reservoirs, and irrigation canals are used to store and distribute the snowmelt runoff during times of the year when the demand of water exceeds its availability.Thus, the Poudre River flow regime is modified.One study of the Poudre watershed described several flow regime modifications including delayed hydrograph rise, decreased peak streamflows, and lower winter base flows (Richer, 2009).In an effort to ensure hydrologic processes are represented appropriately, naturalized streamflows were used for model calibration and testing.Naturalized streamflows remove the influence of afore mentioned features such as diversions and impoundments.Daily naturalized streamflows were collected from Northern Colorado Water Conservancy District at the Mouth of Canyon.Fig. A5 shows the relationship between naturalized daily average streamflow versus observed daily average streamflow.Observed Daily Average Flow, Q obs-avg (cms) A.2. Supplementary Tables Table A3 Matlab code was developed to prepare the land use update files.This code was prepared to add the burned HRUs from the post-wildfire model to the pre-wildfire model and create a land use update file to make a unified model.Land use update files tell SWAT to change the pre-wildfire HRU fractions to nearly zero (The HRU fraction is a HRU level parameter that specifies the fraction of sub-basin area represented by that HRU,Neitsch et al., 2011) and increase the post-wildfire HRU fractions to represent the burn area at the appropriate time during the simulation.In this case, the High Park and Hewlett wildfires occurred during May and June of 2012.Thus, for model calibration the land use update was initiated on July 1, 2012.

Figure 10
Figure 10 displays simulated versus observed monthly streamflows as well as average monthly simulated and observed streamflow for the Mouth of Canyon.This figure suggests that the model slightly overestimates larger monthly streamflows: to a watershed recently exposed to significant wildfire incident located in northern Colorado, USA.The model represents pre-wildfire and postwildfire conditions by implementing the SWAT land use change module as well as curve number updating during simulations to represent burned area as a result of wildfire.Geospatial data representing LULC, soil, terrain, and climate attributes of the study watershed was used to develop the model.An optimal parameter set was obtained for pre-wildfire and post-wildfire conditions through the automated DDS optimization algorithm.Error statistics were calculated to evaluate model performance with regard to daily observed naturalized streamflows.Results indicate a good model performance, with an E NS of 0.82 during calibration as well as 0.71 and 0.88 for the no-wildfire and wildfire testing periods, respectively, for daily streamflows at the Mouth of Canyon.No-wildfire and wildfire scenarios representing a 15 year (2000 to 2014)

Figure 1 :Figure 2 :
Figure 1: Study area map which includes the location of study watershed and the CDWR surface water gauge.

Figure 5 :Figure 6 :
Figure 5: Flow-duration curve at the Mouth of Canyon for the entire simulation period.

Figure 8 :
Figure 8: Hydrologic budgets showing the fate of average annual precipitation (i.e., evapotranspiration, total runoff, and other) with the fate of average annual total runoff (i.e., surface and subsurface) for select sub-basins and the entire study watershed.
Fig. A3.Distribution of burn severity of the Hewlett and High Park wildfires within the study watershed based on MTSB's High Park Fire Assessment.

Fig
Fig. A5.CDWR naturalized daily average streamflow versus observed daily average streamflow with 1 to 1 reference line.

Table 1 .
SWAT model input data.

Table 2 .
Motovilov (1999) between observed and simulated daily streamflows for the calibration period as well as the testing periods.Performance ratings based onMotovilov (1999).5

Table 3 .
Moriasi and Arnold (2007)observed and simulated monthly streamflows for the calibration period as well as the testing periods.Performance ratings based onMoriasi and Arnold (2007).

Table 4 .
Average annual total runoff volumes and depths for both the no-wildfire and fire scenarios, shown for the burned sub-basins as well as for the entire study watershed.Area is also include for reference.

Table 5 .
Ecosurplus and ecodeficit values for the burned sub-basins as well as for the entire study watershed.

Table A1 .
Comprehensive distribution of LULC in study watershed based onNLCD 2001NLCD  , 2006NLCD  , and 2011.5   .5Burnedareas within the watershed were identified using the High Park Wildfire Assessment (Monitoring Trends in Burn Severity Project, 2014) conducted as a part of the Monitoring Trends in Burn Severity (MTBS) project directed by groups within the USGS and United States Forest Service.The MTBS project was introduced to consistently map burn severity and

Table A2 .
Meteorological stations used for this study.

Table A4 .
. Original SWAT database land use / land cover lookup table.Pre-wildfire editied lookup table and corresponding curve numbers.

Table A5 .
Post-wildfire edited lookup table and corresponding curve numbers.TableA6.SWAT calibration parameters.These parameters were varied as a percentage of to maintain spatial variability.