Interactive comment on “ Coupled prediction of flood response and debris flow initiation during warm and cold season events in the Southern Appalachians , USA ”

Debris flows associated with rainstorms are a frequent and devastating hazard in the Southern Appalachians in the United States. Whereas warm-season events are clearly associated with heavy rainfall intensity, the same cannot be said for the cold-season events. Instead, there is a relationship between large (cumulative) rainfall events independently of season, and thus hydrometeorological regime, and debris flows. This suggests that the dynamics of subsurface hydrologic processes play an important role as a trigger mechanism, specifically through soil moisture redistribution by interflow. We further hypothesize that the transient mass fluxes associated with the temporal-spatial dynamics of interflow govern the timing of shallow landslide initiation, and subsequent debris flow mobilization. The first objective of this study is to investigate this relationship. The second objective is to assess the physical basis for a regional coupled flood prediction and debris flow warning system. For this purpose, uncalibrated model simulations of well-documented debris flows in headwater catchments of the Southern Appalachians using a 3-D surface–groundwater hydrologic model coupled with slope stability models are examined in detail. Specifically, we focus on two vulnerable headwater catchments that experience frequent debris flows, the Big Creek and the Jonathan Creek in the Upper Pigeon River Basin, North Carolina, and three distinct weather systems: an extremely heavy summertime convective storm in 2011; a persistent winter storm lasting several days; and a severe winter storm in 2009. These events were selected due to the optimal availability of rainfall observations; availability of detailed field surveys of the landslides shortly after they occurred, which can be used to evaluate model predictions; and because they are representative of events that cause major economic losses in the region. The model results substantiate that interflow is a useful prognostic of conditions necessary for the initiation of slope instability, and should therefore be considered explicitly in landslide hazard assessments. Moreover, the relationships between slope stability and interflow are strongly modulated by the topography and catchment-specific geomorphologic features that determine subsurface flow convergence zones. The three case studies demonstrate the value of coupled prediction of flood response and debris flow initiation potential in the context of developing a regional hazard warning system.


Introduction
The Southern Appalachians have been prone historically to devastating landslides, due to the combination of steep terrain, poorly consolidated colluvium soil mantle, and regional climate (Wieczorek et al., 2009;Radbruch-Hall et al., 1982).The most common and dangerous type of landslide in this region is debris flow (Witt, 2005b), which causes frequent damage to critical infrastructure, in particular roads and private property, and has caused numerous fatalities over the years (Wieczorek and Morgan, 2008).For example, landslide hazard risk assessments indicate that up to 50 % of the area of the Pigeon River Basin in the Southern Appalachians is highly unstable (Witt, 2005a, b).Past climatological attribution studies have established that widespread landslides in the Southern Appalachian Mountains are primarily induced by heavy rainfall (Wieczorek et al., 2009) associated with tropical storms: e.g. the remnants of hurricanes Frances and Ivan in 2004 triggered at least 155 landslides and caused ten fatalities (Wooten et al., 2008).The region is considered a landslide hazard area of high potential (United States Geological Survey, USGS, Fact Sheet 2005-3156), and the USGS has been operating a warning system in the region since 2004 when major hurricanes threaten the area (Baum and Godt, 2010).Note that landslides in remote or uninhabited areas remain undetected until a systematic ground survey or a survey flight is undertaken, thus hindering direct attribution.There is therefore an implicit bias in the interpretation of rainfalldebris flow statistics toward the widespread events associated with summertime tropical systems, which remain short of explaining the over 5000 events mapped so far in the Southern Appalachians.
Historical inventories of landslides and susceptibility maps for the Southern Appalachian Mountains are well documented (Witt, 2005b;Wooten et al., 2008;Wieczorek et al., 2009Wieczorek et al., , 2004;;Wieczorek and Morgan, 2008;Clark, 1987;Fuhrmann et al., 2008).Note there are three modes contributing to debris flow mobilization, namely Coulomb failure, liquefaction and transient/mixed modes of the two (Iverson et al., 1997).The Coulomb failure mode initiates shallow landslide activity which can develop into debris flows.This is the key initiation mechanism in the region of study.Debris flow propagation (post-failure) is not addressed in this study.Forensic surveys and maps of historical events provide critical baseline data for qualitatively assessing and predicting potential debris-flow hazards because there is a higher potential for isolated landslides during heavy rainfall events in the areas along the path of previous debris flows.In addition, deterministic or probabilistic empirical approaches based on rainfall thresholds for predicting landslides through analyzing rainfall intensity-duration characteristics and/or calculating a simplified landslide susceptibility index based on terrain topography have been developed also for the Southern Appalachians (Hong et al., 2007;Kirschbaum et al., 2012;Berti et al., 2012;Guzzetti et al., 2007).However, limited rainfall observations in the past have handicapped the effectiveness of rainfall threshold methods, and the triggering mechanisms inducing slope instability and failure are also controlled by many other factors such as aquifer structure (and water pathways at the soilregolith-bedrock interface), soil characteristics (e.g., soil cohesion, friction angle, particle size distribution), vegetation (e.g., root distribution and cohesion), bioactivity (e.g., worms and burrowers), antecedent soil moisture, and subsurface water movement.Simplified steady-state hydrological models, such as SHALSTAB (Montgomery and Dietrich, 1994) and SINMAP (Tarolli and Tarboton, 2006), take most of the static factors and climate information (e.g., soil properties, slope, vegetation characteristics, soil wetness, etc.) into consideration, and thus can provide climatologically meaningful susceptibility or risk assessments based on high-resolution DEM (digital elevation model) and derived geomorphologic characteristics, but cannot predict the dynamic occurrence of debris flow including the effects of antecedent soil moisture during specific events.Baum and Godt (2010) reviewed early warning systems for shallow rainfall-induced landslides in the USA, which consist of evaluating the likelihood of landslide activity in terms of alert levels (null, outlook, watch and warning) by comparing quantitative precipitation forecasts (QPFs) against rainfall intensity-duration thresholds and antecedent precipitation conditions (soil wetness).The challenges in these early warning systems are the accuracy and lead time of the QPF, the uncertainty in the characterization of geotechnical conditions including land use and land cover, and the relationships among hydrological, hydrogeological and slope stability during individual events.Due to the small areas and steep slopes of headwater catchments in mountainous regions, large rainfall events tend to produce flash flood response and multiple debris flows within the same watershed.From the point of view of public safety and warning systems, predicting the flash flood peak and the location of debris flow initiation is essential, though ultimately the utility of the forecast very strongly hinges on the lead time.Nevertheless, in the case of landslides, such diagnostics can be extremely useful to provide guidance to after-event forensic surveys.In order to better simulate the landslide initiation zone, sensor networks monitoring water levels and soil moisture and/or pressure head in soils can be integrated with threshold warning systems, a strategy that holds great potential to manage clustered hazards in urban centers such as Seattle or San Francisco.However, in remote locations the use of distributed sensor networks for near-real-time assessments of hillslope conditions is not economically and even technically feasible at times.Therefore, predictive models are highly desirable.
The need for coupling dynamically distributed hydrologic models with slope stability models required to quantitatively model or predict the debris flow occurrence both in space and time has been articulated earlier (Safaei et al., 2011;Montgomery and Dietrich, 1994;Baum et al., 2010;Simoni et al., 2008;Iverson, 2000).A widely used modeling strategy consists of using some analytic approximations of Richards' equation coupled with the infinite slope stability model (Taylor, 1948).For instance, the Transient Rainfall Infiltration and Grid-based Regional Slope-stability (TRIGRS) model was developed by Baum et al. (2002) based on a transient rainfall-infiltration model coupled with an infinite slope stability model after Iverson (2000).TRIGRS has been widely applied to study landslides triggered during different types of hydrometeorological regimes (Baum et al., 2005(Baum et al., , 2010;;Godt et al., 2008;Morrissey et al., 2008;Liao et al., 2011;Salciarini et al., 2006).Key limitations of TRI-GRS include the assumption that near-surface soils are saturated or nearly saturated, and are homogeneous and isotropic, and the model is not able to simulate space-time flood response.The latter requires a distributed hydrological model with routing capability.
Several distributed models have been coupled with the infinite-slope stability model including stochastic uncertainty analysis to account for heterogeneity and errors in specified soil properties (e.g., thickness, cohesion, friction angle).For example, GEOtop (Rigon et al., 2006) was combined with an infinite-slope geotechnical model (GEOtop-FS) to simulate the probability of shallow landslide occurrence for saturated conditions, using Gaussian distributions to describe the range of independent parameters and linear uncertainty analysis to estimate their combined effect on the factor of safety (FS; Simoni et al., 2008).Similarly, the HIRESSS (HIgh REsolution Slope Stability Simulator) integrates a hydrological and a geotechnical model, computing pressure head and then calculating the FS, to provide the probability of slope failure given an uniform probability distribution for input parameters using a Monte Carlo technique (Rossi et al., 2013).The Connectivity Index-based Shallow LAndslide Model (CI-SLAM) combines a dynamic topographic index-based hydrological model and an infinite slope stability model (Lu and Godt, 2008) to model shallow landslides (Lanni et al., 2012).Lu and Godt (2008) showed that soil texture heterogeneity and hydraulic properties had a large impact on the timing and depth of the landslides initiation for variably saturated conditions.Similarly, Arnone et al. (2011) used the TIN-based Real-Time Integrated Basin Simulator (tRIBS) with an embeded slope failure method to estimate landslide initiation and performed sensitivity analysis of the model to geotechnical parameters (e.g., soil thickness, cohesion and friction angle) for different rainfall events.For long timescales and from the perspective of landscape management (e.g., timber-harvesting impacts, road construction), a distributed slope stability model (dSLAM), based on a surface-subsurface kinematic wave model including vegetation impacts in terms of root strength and vegetation surcharge, was coupled to an infinite slope stability model to analyze rapid, shallow landslides and the spatial distribution of FS in steep forested basins (Sidle and Wu, 2001;Wu and Sidle, 1995).One common trait of these studies is that the simulated hydrologic response to rainfall forcing (e.g., the flood hydrograph) is not evaluated, and the focus is on the landslide initiation indices or prognostics independently of the underlying hydrologic states.However, Mirus et al. (2007) investigated the role of subsurface flow using a three-dimensional numerical solution of Richards' equation based on the control volume finite-element method combined with an infinite-slope equation (Dutton et al., 2005), and demonstrated that pore-water pressures, and thus slope stability, are underestimated without taking into account convergent subsurface flow.In this study, we will further investigate the critical role of subsurface flow (especially interflow) in triggering the debris flow occurrence.Both the flood response and the debris flow initiation produced by a coupled hydrological-stability model are validated against stream gauge observations and the survey report on the debris flow events provided by NCGS geologists (R. Wooten, personal communication, 2012), respectively.
The record of debris flow events for warm-and coldseason events in the Southern Appalachians reinforces the proposition that heavy rainfall alone and local topography are not sufficient conditions to determine the locations at which debris flows initiate.In this study, we investigate the hypothesis that subsurface flow is closely associated with landslide hazards in mountainous regions through altering the water pore pressure, and thus reducing the shear strength of shallow soils at high elevations and on steep slopes.Previous research has demonstrated that the contribution of interflow to total discharge is dominant (50 % ∼ 70 %) for headwater catchments in the Pigeon River Basin (Tao and Barros, 2013).For this purpose, a dynamical uncalibrated hydrological model (3D-LSHM) was coupled to slope stability models to produce spatially and temporally variable depth-dependent (profiles) of the FS estimates over the Big Creek Basin (BCB) and the Jonathan Creek Basin (JCB) (as shown in Fig. 1), two headwater catchments with a long documented history of landslide activity.Three debris flow events of interest are examined in detail: a prolonged wintertime event and a severe short-duration winter storm that took place around 6-7 January and 8-9 December in 2009, respectively, in the JCB; and a summertime event around 14-15 July 2011 in the BCB, about 15 yr after a similar event at roughly the same location that also caused a flash flood in a neighboring basin.The specific objectives of this study are two-fold: (1) to characterize the hydrology mechanisms leading to rainfall-induced debris flow independently of hydrometeorological regime; and (2) to evaluate and assess the potential utility of prediction the spatial distribution of regional slope instability by coupling a 3-D distributed hydrologic model with slope stability models.The latter should be particularly valuable in the Upper Pigeon and French Broad river basins and in the Southern Appalachians generally, which are undergoing very fast urbanization trends, among the highest in the eastern US.
The manuscript is organized as follows.Section 2 describes the methodology used in this study for detection of debris flow occurrence, including the coupled hydrologic model and slope stability models.Section 3 describes the study area and the landslide events of interest, and the meteorological forcing datasets and ancillary parameters.Analysis and interpretation of results are provided in Sect. 4. In particular, the relationship between interflow and debris flow initiation is discussed in Sect.4.2, and model sensitivity associated with uncertainty in soil internal friction angle and cohesion is investigated in Sect.4.3.Section 5 provides a summary and conclusions.

Methodology
In this study, the focus is on the coupled simulation of flood response and debris flow initiation, with an emphasis on the role of hydrologic processes, and interflow in particular, in the redistribution of infiltrated rainfall in the landscape.For this purpose, a distributed hydrological model (3D-LSHM) was coupled with two different approaches to detect slope instability: (1) an infinite slope stability index (SSI) method modified after a widely used deterministic model (SHALSTAB; Montgomery and Dietrich, 1994), and (2) a dynamic FS model derived using the limit equilibrium method accounting for both unsaturated and saturated soil moisture conditions and including interflow.The objective is to simulate spatio-temporal distributions of SSI and FS at high spatial and temporal resolutions to detect potential locations for rainfall-induced debris flow initiation in headwater basins in the Southern Appalachians, which could be integrated with a quantitative flash flood forecasting framework to improve the effectiveness of regional early warning systems.

Land surface hydrology model (3D-LSHM)
A fully distributed and physically based three-dimensional land surface hydrologic model (3D-LSHM) (Yildiz andBarros, 2009, 2007;Tao and Barros, 2013) is used to solve the coupled water and energy balance equations including coupled surface-subsurface interactions.The temporal and spatial resolution of model simulations is 5 min and 250 m, respectively, which meets numerical stability requirements, and reflects a compromise among the coarse spatial resolution of the atmospheric forcing datasets (1-32 km), the spatial scale of terrestrial ancillary data such as soils properties and vegetation cover (∼ 1 km), and the spatial resolution adequate to capture the governing hydrologic processes (e.g., Tao and Barros, 2013).Each grid element in the modeling domain represents a vertical soil column initially consisting of both an unsaturated and a saturated zone.The unsaturated zone is discretized into three soil layers that also serve as root layers (the number of discrete soil layers used to solve the equations numerically can be significantly larger).The first soil layer in the unsaturated zone functions as the landatmosphere interface.At each grid element, overland flow is first estimated either from infiltration excess or saturation excess mechanisms at each time step, and then routed downslope by a surface flow routing module that relies on a 1-D kinematic wave approximation, assuming a linear flow surface across grid cells (Yildiz and Barros, 2007).The Green-Ampt method is used to describe infiltration.Although the model is equipped with a Richards' equation solver, it is not utilized here.The hydraulic conductivity that governs the gravitational mass flux when the soil moisture is above field capacity follows Campbell (1974).Subsurface flow, comprising interflow and baseflow, is estimated at each grid element in each soil layer, and then is routed to channel segments by a lateral subsurface flow routing module using a modified multi-cell approach (Bear, 1979).The Muskingum-Cunge method of variable parameters (Ponce and Yevjevich, 1978) is utilized for the channel flow routing without significant backwater effects.Sensible and latent heat fluxes are estimated based on the Monin-Obukhov similarity theory, which provides dimensionless variables expressing the buoyancy effects resulting from the vertical density gradients in the stable atmosphere with modifications for unstable boundary layer conditions, and are calculated using the input air temperature, air pressure, wind velocity and specific humidity.Radiative forcing is calculated based on the input downward shortwave/longwave radiation from the atmospheric forcing dataset, and landscape attributes such as albedo and emissivity.Further details on the representation of land-atmosphere interactions in the model are described in Devonec and Barros (2002).Detailed description and applications of the model can be found in Yildiz andBarros (2007, 2009) and Tao and Barros (2013).Note that, in principle, the higher the spatial resolution the more rigorous the coupling between the hydrological and slope stability models, and the more accurate representation of governing processes and spatial gradients.Therefore, a scale effect is expected with simulated hydrologic variables displaying smoother spatial distributions at coarser model resolutions (see, for example, Yildiz and Barros, 2009).

Stability index mapping
Shear strength testing of soils in the Southern Appalachians indicates that soils in debris flow initiation zones in the region of study either are cohesionless or have very low cohesion (Witt, 2005b).Based on the assumption that the water table follows topography at small scales, and thus is parallel to the slope, and that the soil material is cohesionless, Dietrich et al. (1993) proposed a simplified infinite slope stability model using the conventional limit equilibrium method (i.e., at equilibrium, driving forces are equal to resisting forces): where ρ s and ρ w are soil bulk density and water density, respectively; θ is local slope angle; ϕ is soil internal frictional angle; h is the saturated soil depth; and z is the total soil depth to bedrock.Equation (1) then is used to map SSI (x, y, z) , based on the basin topographic characteristics as per Montgomery and Dietrich (1994).Here, Eq. ( 1) was modified to incorporate simulated soil wetness by the 3D-LSHM since the term [h/z] is equivalent to the saturation degree [W = w/φ], where w is the simulated volumetric soil moisture and φ is soil porosity: Thus, the slope stability classification will be performed using spatio-temporal soil moisture distributions.The SSI classes range from unconditionally unstable to unstable, stable and unconditionally stable (as shown in Fig. 2).The obvious merit of this effective and simple approach is that it requires a small number of parameters and state variables, such as soil internal friction angle, slope angle and wetness, as input fields.This is a significant advantage compared to methods requiring many parameters that are not easy to measure and thus induce uncertainties, especially over topographic complex regions.Even though the SSI method is based on quantitative information relating soil moisture and slope static properties, this information is aggregated into broad qualitative categories using a threshold-base classification, which creates ambiguity as many different slope states belong to the same category.That is, locations classified as unstable or unconditionally unstable pixels are highly susceptible to debris flow, but it is not necessary that debris flow will initiate.

Factor of safety
Although the SSI can provide spatio-temporal instability information, it neglects the soil cohesion and suction effects, as well as the relational position of a grid element with respect to its neighbors.In order to quantitatively analyze the debris flow triggering mechanisms accounting for all the dominant factors, the spatio-temporal distribution of a factor that can represent the dynamic net forces acting on the slope should be determined explicitly.The factor of safety (FS), defined as the ratio between resisting forces and the driving forces, is a widely used factor for slope instability analysis.
In this section, a dynamic form of the FS equation using the limit equilibrium method for an infinite slope was derived as described next.
Figure 3 shows the diagram of forces acting at a generic location x and depth z on a cross section of the conceptual infinite slope model.The z axis is normal to the surface and positive in the downward direction.The x axis is parallel to the slope surface, and thus normal to z. F N and F P are the normal component and parallel components of the weight G along the slope, respectively.The normal component of the weight is counteracted by the normal resisting force N. In the along slope direction, F P is counteracted by the friction F f , suction F s and cohesion F c forces.The Coulomb failure mode occurs when the shear stress at failure on the failure plane equals or exceeds the resultant of friction, suction and cohesion stresses, that is, FS ≤ 1.The infinite slope model has two critical assumptions.First, it assumes that the slope failure occurs within a thin soil layer, and second it assumes that the failure plane is of infinite length, i.e., H L, in Fig. 3.In this study, the "effective" L is the spatial resolution of the model (250 m), whereas H is the soil mantle thickness, which is a spatially variable ranging 10s of centimeters at higher elevations to 100s of centimeters at lower elevations and in the valleys; thus L/H is always larger than the critical ratio of 25 above which the infinite slope assumption is valid (Milledge et al., 2012).
For unsaturated conditions, as shown in Fig. 3, the limit equilibrium equation should be written as where θ is the slope angle, F N is the normal component of gravity G(z, t) = γ s (z, t) Az with γ s being the depth-averaged soil-specific weight, z is the depth below the surface, and A is the nominal area where the force is applied (i.e., the area of the slice shown in Fig. 3), that is, the spatial resolution in our model; F P is the parallel force to the surface due to gravity; F f , F s and F c are resisting forces due to soil friction, soil suction pressure P w (z, t) and cohesion due to both soil and vegetation, expressed as follows (Rossi et al., 2013): As stated earlier, ϕ is the soil internal frictional angle.Here the tanϕ is the friction coefficient.γ w is the specific weight of water; c is the combined measure of soil and vegetation cohesion; and (z, t) is the pore suction head distribution in space and time.Instead of using an analytical approximation (e.g., Lu and Godt, 2008), the dynamic pressure head profile is simulated by the physically based 3D-LSHM according to the dynamic soil moisture characteristic curve as described by the soil water retention equation (Campbell, 1974): where λ is the pore-size index and b is the bubbling capillary pressure head.The parameters λ and b are assigned according to soil texture (Rawls et al., 1982).Other relevant model parameters are discussed in Sect.3. At the parallel direction to the surface, the resisting forces include the friction force, suction force and cohesion force, whereas the driving force is the gravitational force.The FS is equal to one when the slope is at equilibrium.Rearranging Eq. ( 4) to separate the resisting forces from the driving forces in the direction parallel to the slope, and then dividing the two sides of the equation by the driving forces, we can obtain the final form of the FS equation for unsaturated conditions: Typically, the shear stress induced by the water flow is neglected due to the small equivalent kinetic energy head (V 2 /2 g) caused by subsurface flow in each layer.However, here it is explicitly incorporated in the F s term in Eq. ( 6).
For saturated conditions, the suction force vanishes; instead a hydrostatic force F H (z, t) will act on the slope, and Eq. ( 3) can be rewritten as where F H (z, t) = ρ w ghAcosθ and h is the depth of fully saturated soil at soil depth z.Rearranging Eq. ( 7) yields the equation of FS for saturated conditions (note the kinetic energy head is included in the second term): and γ sat is the specific weight of saturated soil.
During rainfall-triggered landslide events, the soil pore water pressure on steep slopes increases towards positive suction head, reducing the suction force and then shear strength.Meanwhile the shear stresses increase and cohesive resistance decreases as the soil becomes wet, causing the slopes to become unstable.When shear strength exceeds shear stress, i.e., resisting force is larger than driving force, FS > 1 and the slope remains stable.When FS < 1, the slope fails.Equations ( 6) and ( 8) account for the essential processes that play interactive roles in the initiation of debris flows.In this study, values of basic soil properties were extracted and compared against previous studies in or near this region (Witt, 2005b;Liao et al., 2011), and are summarized in Table 1.

Case studies
Three case studies are conducted in this work: two coldseason events, on 7 January and 9 December 2009 in the JCB, and a warm-season event, on 15 July 2011 in the BCB.However, neither the BCB nor the JCB is equipped with stream gauges.Consequently, streamflow simulations for the same three events were also conducted for the Cataloochee Creek Basin (CCB), a USGS hydrologic benchmark watershed and the closest watershed to the BCB and the JCB (Fig. 1), for hydrological verification and to demonstrate the robustness of the estimated rainfall fields.Verification of the location of landslide initiation is based on the survey data provided by North Carolina Geodetic Survey (NCGS; R. Wooten, personal communication, 2012). of the 2nd and 3rd layers are from 0.5 to 1.5 m varying with elevation and slope (Fig. 6) Ksat (m s −1 ) Spatially varying (Fig. 7) Scaling factors for Kv None Scaling factors for Kh 1000-300-1-0.1 Porosity (m 3 m −3 ) Spatially varying (Fig. 7) Field capacity (m 3 m −3 ) Spatially varying (Fig. 7) Wilting point (m 3 m −3 ) Spatially varying (Fig. 7) Channel cross section Rectangular, channel width ranging from 1 to 30 m Channel threshold (pixels) 5 Parameters for slope stability models Values Soil density (kg m −3 ) 1922 (Witt, 2005b) Soil friction angle (degree) 26 (Witt, 2005b) Soil and vegetation cohesion (Pa) 2000 (Witt, 2005b)

Study area
The Big Creek Basin (BCB), the Cataloochee Creek Basin (CCB) and the Jonathan Creek Basin (JCB) are three headwater catchments in the Pigeon River Basin, in the Southern Appalachians in North Carolina, USA.The Cataloochee Creek is a small tributary to the Pigeon River and has a drainage area of 128 km 2 .The BCB and JCB have drainage areas of about 95 km 2 and 172 km 2 , respectively.The three headwater catchments are heavily forested and are characterized by steep slopes.In recent years, the JCB has witnessed significant land-use and land-cover (LULC) change due to increased urbanization.The Pigeon River Basin is underlain by crystalline-rock aquifers comprising crystalline metamorphic and igneous rocks covered by an extensive mantle of unconsolidated material consisting of saprolite, colluvium, alluvium, and soil (Trapp Jr. and Horn, 1997;Miller, 1999).The colluvial deposits are mainly found on the hillsides due to rock weathering and are highly susceptible to landslides.Substantial alluvial deposits appear along streams and are built over time due to sediment transport in the streams.The dominant soil types are Edneyville-Chestnut complex soil, Plott fine sandy loam, Wayah sandy loam, and eroded Wayah loam soil (Allison et al., 1997).The climate for the study area is subject to moisture-rich winds from the Gulf of Mexico and westerly mesoscale convective systems in the warm season, whereas westerly and northwesterly flows govern most of winter weather activities.Previous research has shown that the orographic rainfall enhancement is very strong, on the order of 60 % at ridge compared to valley locations (Prat and Barros, 2010).The rainfall threshold for debris flows based on the historical record is 125 mm over a 24 h period (Witt, 2005a).However, recent observations such as during the July event studied here indicate that such rainfall can accumulate in periods of less than 90 min (Prat and Barros, 2010;Tao and Barros, 2013).Existing landslide hazard risk assessments indicate that most of the area of the Pigeon River Basin is highly unstable, especially the headwater catchments (Witt, 2005b, a).

Landslide events
Basic geologic and geomorphic conditions at the debris flow sites for the three landslide events investigated here are summarized in Table 2.

Warm season
In the middle of July in 2011, an extremely heavy storm event triggered debris flows in the Big Creek Basin (Fig. 1), which also caused flash flooding that damaged the Cherokee fish hatchery on the evening of 14 July (Lee et al., 2011).Observational field data of soil and rock materials were collected by NCGS geologists at three debris flow locations.The debris flow tracks were scoured to bare soil and, in some locations, down to the underlying bedrock for most of their lengths.Most vegetation was stripped or downed along the tracks, including large trees.It was determined that there is a potential for further slope movements originating from the steep slopes in the head scarp regions of all three debris flows.

Cold season
A debris flow in the JCB was caused by a severe winter storm on 8-9 December 2009.Although the slope failure narrowly missed some infrastructure, it destroyed a portion of the Rich Cove Road, and one lane had to be removed.Another debris flow near Bear Trail in the JCB was triggered by a persistently heavy rainfall system during 5-8 January 2009, cleared all vegetation in its path, eroded away a large tract of a local road, destroyed a private home and caused personal injuries 1 .NCGS geologists visited the initiation site several times and established that the debris flow initiated at a colluvium catchment with localized residual deposits filling in-between the colluvium and overlying bedrock located on a north-facing steep slope (Fig. 1), and that bedrock controlled locally the geometry of the initiation zone, a common characteristic in the Southern Appalachians (Sas Jr. and Eaton, 2008;Wooten et al., 2008Wooten et al., , 2009)).

Forcing data and model parameters
One factor that limits landslide and hydrologic studies in mountainous regions, independently of the modeling approach, is that to predict dynamically the initiation of debris flow on an event basis, availability of good-quality spatiotemporal rainfall distributions at the resolutions required to capture the subsurface physics of soil wetting and water flow processes is critical.A spatially dense, high-elevation rain gauge network has been recording observations in the Upper Pigeon River Basin in the Great Smoky Mountains since 2007 to investigate the 4-D distribution of precipitation in the region (Prat and Barros, 2010).These rain gauge observations have been used to characterize the spatial-temporal error structure of radar-based quantitative precipitation estimates (QPEs) and to improve QPEs for hydrological modeling with success (Tao and Barros, 2013).The recent record of landslide activity in populated areas indicates that debris 1 http://www.ncdc.noaa.gov/stormevents/eventdetails.jsp?id= 151707 flows are all-season events in the region, and that mesoscale convective systems and isolated thunderstorms play an important role in concurrent flash flooding and debris flows in the warm season.Therefore, assessing the quality of rainfall data and bias correction or adjustment procedures to improve the accuracy of rainfall input to the hydrological system is necessary.This is discussed in detail in Sect.3.3.1.Other meteorological forcing datasets and model parameters for analyzing slope stability are discussed in Sect.3.3.2.

Rain gauge observations
A spatially dense, high-elevation rain gauge network has been recording observations in the Pigeon River Basin to investigate the 4-D distribution of precipitation in the Great Smoky Mountains (GSMRGN) since the summer of 2007.The network comprises 35 stations at elevations ranging from 1150 to 1920 m along exposed ridges in the Southern Appalachians (purple circles in Fig. 1 (Prat and Barros, 2010;Tao and Barros, 2013)).More detailed information about the network can be found at http://iphex.pratt.duke.edu/.
Similar to Tao and Barros (2013), GSMRGN observations are used here to assess and improve existing radar-based QPEs required by the distributed hydrologic model, specifically the Q2 product described below.Only rain gauges along the topographic divide and within individual basins are used for assessing and correcting Q2 over each particular basin.For example, rain gauges numbered 3## were used for the Big Creek Basin, as shown in Fig. 1.Overall, 12 rain gauges in total were used for the CCB and for the BCB during the 12-17 July 2011 event, and 9 rain gauges were used for the JCB during 5-10 January 2009.There were 7 rain gauges for the CCB during the event in January of 2009 but 12 rain gauges during the event in December of 2009, because 3## rain gauges were not installed until the summer in 2009.Availability of rainfall observations is one major reason why these three most recent events were selected for this modeling study.

QPE adjustment
The experimental Next-Generation Multi-sensor QPE was obtained from the National Mosaic and Multi-sensor QPE (NMQ) project at the National Severe Storms Laboratory (NSSL).The local gauge-corrected hourly radar-based QPE product (Q2RAD_HSR_GC, Q2 in short) at high spatial resolution (0.01 • × 0.01 • ) (Vasiloff et al., 2007;Zhang et al., 2011) was used in this study.We first evaluate the Q2 datasets using the GSMRGN rain gauge observations to characterize the spatial-temporal error structures in Q2, and then apply bias correction to improve the accuracy of Q2 based on the error structures identified.
Hourly Q2 accumulations were spatially interpolated using a nearest-neighbor method to downscale rainfall fields from 1 km resolution to higher resolution at 250 m, and the downscaled values were subsequently compared to rain gauge observations at the grid scale.Figure 4 shows that the original Q2 fields generally underestimate rainfall for the summer storms (Fig. 4a and b) and winter storms (Fig. 4cf), despite the large difference (on the order of one order of magnitude) in rainfall intensity between the events.The inaccuracies in Q2 are attributed mainly to radar-terrain configuration issues (e.g., radar beam blockage or overshooting) and the radar rainfall retrieval algorithm (Young et al., 1999;Smith et al., 1996;Fulton et al., 1998;Prat and Barros, 2009).
A simple bias-correction adjustment method based on linear regression relationships between hourly rain gauge observations and Q2 data was developed and was successfully applied previously to adjust Q2 for a tropical storm (Tao and Barros, 2013).We employed the same procedure to improve the Q2 accuracy at small basin scale in this study, taking advantage of the very high-density GSMRGN observations in the Pigeon River Basin.The adjusted Q2 product demonstrates significant improvement compared to the original Q2 (Fig. 4).The root mean square error (RMSE) (mm h −1 ) between observed rainfall rate (mm h −1 ) and the Q2 product before and after adjustment is provided in Table 3.The adjusted Q2 outperformed the original data over all the basins for the three rainfall events with RMSEs reduced significantly, resulting in large improvement in storm cumulative rainfall amounts (Fig. 5).The limitation is that no rain gauges are installed in the inner basin to characterize the error structure of Q2 in the valley.This data void may cause uncertainty in the areas at lower elevation.The downscaled and adjusted hourly Q2 fields were interpolated to five-minute temporal resolution using the methodology described by Tao and Barros (2013), where further discussion and a detailed description of the Q2 spatial error structure in Q2 can also be found.

Ancillary data and model parameters
The soil internal friction angle and cohesion are two important parameters required by slope stability models.Uncertainty in these parameters can induce very large uncertainty in the resultant FS values.The present study benefited from previous research conducted by Witt (2005b) in the same area.Witt examined the same two parameters using SIN-MAP and SHALSTAB, and reported representative values.We adopted those representative values as well as the ranges reported by this study for slope stability analysis and sensitivity analysis (shown in Table 1).
The 3D-LSHM needs basin geomorphic information (such as elevation, slope and flow direction), soil hydraulic and geotechnical parameters, landscape attributes and atmospheric forcing datasets to perform the hydrological simulations.The DEM over the Pigeon River Basin was obtained from the National Elevation Dataset (NED) provided by the US Geological Survey at 3arcsec resolution, and was reprojected and spatially resampled to the model grid at 250 m resolution.The spatially varying soil depth was estimated using two alternative approaches, the Z-and S-methods (Saulnier et al., 1997): where h(x, y) is the estimated total depth of the middle and deep soil layers at pixel (x, y); h max and h min are the maximum and minimum depths, respectively; and z max and z min , and θ max and θ min are the corresponding maximum and minimum elevations and slope angles.The Z-method assumes that soil depths increase as elevation decreases.The S-method assumes that soil depths increase as topographic slope decreases, because soil cannot accumulated on steeper slopes due to erosion and landslides.However, the Z-method tends to result in too-thin soil depth at very high elevations, while the soil depth in the valleys calculated by the S-method tends to be too thick (Fig. 6) based on the authors' observations in the field.Consequently, the mean of the soil depth estimated by both methods is adopted in this study.The ratio of the thickness of the second to the third layer is 2 : 3, to roughly represent root density distribution.
The top soil layer in the model is fixed as 10 cm all over the basin.The total depth of the second and third layer is the mean of soil depths estimated by the Z-method and Smethod, with h max = 1.5 m and h min = 0.5 m, and θ max = 40.96degree and θ min = 0.01 degree.The base layer is 1m deep at elevations above 1300 m, and 4 m deep below 1300 m to represent thicker alluvial deposits in the valleys.
Soil hydraulic properties including saturated hydraulic conductivity, porosity, field capacity and wilting point were extracted from the State Soil Geographic (STATSGO) database provided by the US Geological Survey (Schwarz and Alexander, 1995).Standard soil layers defined in STATSGO were selected according to soil depth for modeling layers at each pixel first.Then, soil parameters for each model layer were extracted from the STATSGO layers taking into consideration the depth of the soil column for each grid element.The minimum value of vertical saturated hydraulic conductivity from STATSGO was used as a representative for each soil layer since the minimum hydraulic conductivity controls the hydrological response.For other soil properties, such as porosity, field capacity and wilting point, average values were used.It must be stressed that all soil hydraulic parameters are spatially varying across the basin as shown in Fig. 7, which displays large heterogeneity in 3-D space.Although there is certainly uncertainty in these soil properties, which in turn affect the calculation of both the hydrological response and also the slope instability analysis, we did not perform further investigation addressing these uncertainties, and assumed that the values extracted from the STATSGO are physically based and are representative of the actual soil properties in the region as in Tao and Barros (2013).
Space-time varying landscape attributes such as broadband albedo, broadband emissivity, fractional vegetation coverage, and leaf area index were derived from NASA's MODIS (Moderate Resolution Imaging Spectroradiometer) products (MCD43B1, MOD11C2, and MCD15A2, respectively).The original products were first reprojected, bilinearly interpolated to the model grid, and then linearly interpolated into five-minute temporal resolution.Missing data gaps are addressed using physically meaningful constraints based on ancillary data.Lastly, quality-control adaptive temporal filtering for the landscape attribute data was performed using TIMESAT software to reduce the discontinuity caused by cloud contamination following the adaptive Savitzky-Golay filtering method (Eklundha and Jönssonb, 2012).
The meteorological forcing data required by the hydrological model were extracted from NCEP (National Centers for Environmental Prediction) North American Regional Reanalysis (NARR) products originally at 32 km spatial resolution and 3 h temporal resolution (Mesinger et al., 2006), including air temperature, air pressure, wind velocity, downward shortwave and longwave radiation and specific humidity.The bilinear interpolation method was utilized to interpolate NARR fields to finer spatial resolution at 250 m, and linear interpolation was applied in time.Elevation adjustments and corrections to near-surface variables were applied between NARR terrain and local terrain at high resolution for each time step based on predicted atmospheric conditions (e.g., using dynamic lapse rates).Special bias corrections for downward shortwave radiation were applied through dynamical adjustment, accounting for cloudiness and topographic effects.The atmospheric forcing and landscape property datasets are subsets from the high-resolution datasets developed to provide the hydrologic modeling/forecasting for the southeastern US, in support of the Integrated Precipitation and Hydrology Experiment (IPHEx, http://iphex.pratt.duke.edu/).Other ancillary parameters were specified according to prior studies (Dickinson et al., 1993;Chow, 1959;Campbell, 1974;Clapp and Hornberger, 1978;Jackson, 1981;Yildiz and Barros, 2005, 2007, 2009;Price et al., 2010Price et al., , 2011;;Tao and Barros, 2013).

Hydrological verification over the Cataloochee Creek Basin (CCB)
The 3D-LSHM was used first to simulate hydrological response to the storm events in January of 2009, December of 2009 and July of 2011 over the CCB.Model-simulated streamflows were compared against stream gauge observations to evaluate the model's hydrologic performance.Initial conditions and essential model parameters are provided in Table 1.In order to allow the model state variables to reach internal consistency, model spin-up simulations for the same duration of the simulation were conducted before the event simulation proper.The end of the spin-up period is the beginning of the event simulation.The basin soil moisture conditions for the spin-up simulations were initialized by specifying soil wetness based on seasonal climatology modeled to be consistent with the streamflow at the beginning of the simulation period (shown as in Table 1).
The comparison between the discharge observations and simulated hydrographs over the CCB generated from the 3D-LSHM driven by Q2 rainfall datasets before and after adjustment for the three events are presented in Fig. 8.The original Q2 fields significantly underestimate rainfall, yielding much lower streamflow and completely missing the storm response (Fig. 8a.1, b.1, and c.1).By contrast, the simulated streamflows using the adjusted Q2 forcing show very good agreement with the stream gauge observations with regard to the peak flow and peak time of the hydrographs, as well as the shape of the rising limb of the hydrograph, for the summer and the severe winter storm simulations (Fig. 8b.1 and c.1).For an extreme event with high rainfall intensity such as the summer storm in July of 2011 or the severe winter storm in December of 2009, large overland flow is produced concurrently with the heavy rainfall as illustrated in Fig. 8a.2 and c.2.Nevertheless, despite the strong and fast response of overland flow, the interflow produced by subsurface soils is the governing contribution to the basin's hydrological response by water volume.For a prolonged and persistent rainfall event such as the winter storm in January in 2009, interflow plays a governing role in the hydrological response with regard to flow rate and flow volume determining the peak time and overall shape of hydrograph, as illustrated in Fig. 8b.2.Compared to the large interflow produced from the top soil layer in the summer and the severe winter storm, the interflow in the second soil layer is dominant in the prolonged winter storm event, consistent with persistent rainfall lasting for several days.Overall, interflow dominates the flow processes and determines the water redistribution in the basin, which is of vital importance for the initiation of shallow landslides and debris flow mobilization.

Landslide event analysis
Hydrological evaluation of simulated streamflow against observations for the CCB indicates that the estimated rainfall forcing is robust and the specified model parameters are representative for the region.Thus, the same parameterization and datasets are used to implement the model for the BCB and the JCB watersheds.

Warm-season events -debris flow in Big Creek
Basin (BCB) The 14-15 July 2011 debris flows in the BCB were associated with the passage of a nocturnal convective rainfall system, and one particular convective cell that remained stationary for nearly two hours at the ridge above the location where the debris flows initiated.The short-duration but severe rainstorm produced large amounts of precipitation on the border between the BCB and CCB (as shown in Fig. 5), with rainfall rates as high as 60 mm h −1 .The three debris flow initiation zones mapped by the NCGS are located in three nearby pixels in the modeling domain.The time series of volumetric soil moisture and interflow produced at each soil layer for one of the pixels in which debris flow occurred are shown in Fig. 9a.Due to similarity, the plots for the other two pixels will not be shown here.Note that the negative flow rates indicate the flow is leaving the pixel; in other words, the combined infiltrated rainfall at the pixel and incoming flow received from upstream locations is smaller than the outflow.As can be seen from the figure, the top two layers respond promptly to rainfall infiltration and produce large interflow.The dash line indicates the time when the magnitude of total interflow reaches its peak, which is concurrent with the time when the debris flow occurred.The spatial distribution of soil moisture, absolute interflow magnitude at each soil layer as well as the total interflow are shown in Fig. 10a.The debris flow locations marked by circles show very large interflow compared to other locations with steeper slopes and nearly the same rainfall where debris flows did not initiate.The histograms of soil moisture, interflow, slope and rainfall rate for the three events are shown in Fig. 11.The solid red lines indicate the local conditions in the unstable grid element selected for analysis (corresponding to the gray solid line in the upper interflow time series).The histograms of these variables provide an alternative view of the same data that illustrates the concurrency of slope steepness, high rainfall intensity, and large and fast interflow response especially from the top two layers at the unstable locations.As can be seen from Fig. 11a, the histograms of  rainfall and total interflow show uneven distributions, skewed to the left and with very long tails on the right.For the conditions when and where the shallow landslide initiated the debris flow (marked by the vertical solid lines), rainfall, total interflow and the slope as well show intermediate high values on the right of the distribution.Figure 12a highlights the relationship between interflow and the initiation of slope failure.The temporal evolution of vertical profiles of soil moisture, pressure head, interflow and the FS are presented in Fig. 13a.When soil moisture increases, the (negative) suction pressure head increases, leading to a decrease in FS as the slope becomes less stable.Note that the kinetic energy head included in the estimation of FS is potentially grossly underestimated here because the relationship between particle size and pore size distributions, and the distribution of soil pipe networks are not accounted for in determining the effective hydraulic area and interflow pathway system.The inset shows that at around 02:30 UTC on 15 July the value of FS crosses the theoretical stability-to-instability line (FS = 1) in the base layer.For this warm-season event, the debris flow coincidently occurred at the location and time where and when the heaviest rainfall occurred.This is different from the case in the JCB for the cold-season event presented next.
Figure 14a depicts the spatial distributions of SSI and FS at the time of debris flow initiation (dash line in Fig. 9).The three pixels where debris flow took place are classified as unconditionally unstable by the slope stability index mapping method, meaning they are highly susceptible to debris flow.The slope instability simulated by FS is also below unity at the three pixels, indicating unstable conditions toward slope failure, consistent with NCGS field surveys.The number of total unstable pixels identified in the basin by the SSI and the FS metrics varies with time (Fig. 15), mimicking closely the spatial distribution of interflow and the space-time evolution of the storm system.Note however that the number of unstable pixels is almost one order of magnitude larger using the SSI method, which suggests that it overestimates the extent of unstable areas, mainly because it neglects the cohesion and suction effects in the soil.This is consistent with Witt (2005b), who reported that 80-90 % of the region is highly unstable and susceptible to the debris flow occurrence using a static SSI (Dietrich et al., 1993), a clearly excessive estimate based on the historical record.Finally, debris flow proper is not simulated in this study, and therefore the simulation is not representative of realistic conditions after initiation.It is expected that, following mass movement, the shear stresses at locations surrounding the initiation points will decrease, and thus there should be a strong decrease in the number of unstable pixels.

Cold-season events -debris flow in Jonathan
Creek Basin (JCB) The first debris flow examined here for the JCB was caused by a persistent rainfall system from 5 to 10 January 2009.This winter storm produced persistent rainfall of moderate intensity (< 10 mm h −1 ) typical of stratiform orographic systems for about two days continuously.The second debris flow was triggered by a severe winter storm on 8-9 December 2009, presenting relatively larger rainfall intensity but lasting for about one day.The time series of soil moisture and interflow for each soil layer at the pixel where debris flow occurred are presented in Fig. 9b and c.As the persistent rainfall infiltrates, it is stored at first in the top two layers, which produced relatively small interflow.When the second soil layer finally reaches saturation, interflow increases rapidly and reaches the peak value as indicated by the dash line in the bottom panel (Fig. 9b.2).The same situation is found for the severe winter storm (Fig. 9c).It should be noted that the timing of peak interflow for the severe winter storm in December is not concurrent with the rainfall.It occurs about two hours after the rainfall ends, indicating that the subsurface flows redistribute water and take some time to concentrate at this point.Note the interflow in the top layer is positive for the persistent winter storm (Fig. 9b.2), meaning the top layer overall is receiving more water from incoming interflow from upslope areas and rainfall input than the interflow it releases as outflow at this pixel.The second soil layer has negative interflow, meaning that the net interflow is leaving the pixel.The opposite flow directions in the soil column contribute to a more complex shear stress profile than just the gravitational stresses.These results reinforce the premise that interflow plays an important role in destabilizing slopes.
The spatial distributions of soil moisture and interflow at the interflow peak time are shown in Fig. 10b and c for the two winter storms.Note that the interflow shown in the spatial map is the absolute magnitude, emphasizing the impact of interflow rates on the slope stability.Both pixels show relatively large interflow compared to the surrounding pixels, indicating a concave area concentrating subsurface flow.For the event in January, the basin received some rainfall at the end period of the prolonged storm system at the interflow peak time.However, for the severe winter storm in December, the rain ceased two hours before the interflow reached the peak as mentioned earlier.This fact illustrates that, to determine the initiation time for debris flow, considering rainfall alone is not enough because the most important mechanism controlling the process is subsurface flow, particularly the interflow.
Figure 11b shows the histograms of soil moisture and interflow, with the conditions at two different times marked by the red dash and solid lines corresponding to the vertical gray dash and solid lines shown in the interflow time series (top panel) for the persistent winter storm in January of 2009.Both the red and gray dash lines refer to the condition when the largest rainfall rate took place at the pixel; both the red and gray solid lines refer to the condition when the interflow reaches the maximum at the pixel.This case is representative of conditions when rainfall thresholds are not a necessary condition for debris flow initiation.Rather it is the interactions among antecedent soil moisture and interflow that differentiate this condition (red solid line, which coincided with the debris flow initiation), from the heaviest rainfall condition (red dash line).This is also clearly shown for the severe storm in December (Fig. 11c).It is the nonlinear interactions among steep slope, antecedent soil moisture conditions, and basin-received rainfall that lead to the production of large and fast interflow that destabilized the slope, which is again consistent with the hypothesis articulated in Sect. 1. Figure 12b and c illustrate how the geomorphology of the JCB such as the concave landform (as shown in Table 2) with modest slopes at intermediate elevations favors interflow concentration in the debris flow initiation zone, similar to the findings reported by Mirus et al. (2007).This is in contrast with the BCB (Fig. 12a).Nevertheless, as Fig. 13b and c show, even though the third and base soil layers reach saturation, the FS remains slightly above unity (FS = 1.04 for the event in January and FS = 1.10 for the event in December), and thus the soil column would be classified as stable.Given the uncertainty in specifying soil properties and in capturing soil structural heterogeneity, it is important to recognize throughout  Note that in this work, the debris flow proper is not simulated, and therefore the simulation is not representative of realistic conditions after debris flow initiation.For example, it is expected that, with mass movement, the shear stresses at locations surrounding the initiation points will decrease, and thus there should be a strong decrease in the number of unstable pixels.
our discussion that FS estimates are also uncertain.This is addressed in part through the sensitivity analysis presented in Sect.4.3 below.Unlike the warm-season event in which the debris flow occurred in the wetting period, the profiles for the event in December (Fig. 13c) demonstrate that the top layer was already in a drying period, while the bottom layers received drained water from the upper layer and upstream pixels, not from infiltrated rainfall.This is very important for issuing debris flow warnings.Debris flow could still occur, especially at these concave areas in the basin, after rainfall has stopped.
Figure 14b and c display the spatial distributions of SSI (top panel) and FS (bottom panel) at the time when the debris flow was initiated.The initiation location is unambiguously identified as unstable using the SSI method, but not for the FS as expected based on the FS profile presented in Fig. 13b and c.Nevertheless, immediate neighbors at higher elevation do exhibit FS values below unity, and thus become unstable at the critical time.This begs the question of spatial uncertainty, which can be associated with the rainfall forcing, soil depth, soil hydraulic properties, root and soil cohesion, etc.As in the summer case, the number of unstable pixels using the SSI metric is larger by almost two to ten times than that for FS (Fig. 15b and c).The trends of unstable pixels identified by SSI closely follow the change in soil wetness, due to its intrinsic dependence on instantaneous soil moisture in the basin.The FS metric tends to be a more conservative (and realistic) approach to detect slope failure.The number of unstable pixels still increased after the rainfall ceased, illustrating that subsurface flow continued redistributing water to concave areas in the basin.Yet, there is still a large number of unstable or nearly unstable locations at each time, which is an indication of spatial ambiguity.On the other hand, Figs.13b and 12b show that interflow peaks locally at the time of initiation, which can be used as an additional constraint in assessing local instability.Overall, the results highlight the role of interflow in slope mobilization for these cold-season events.

Sensitivity analysis
Whereas a full-fledged uncertainty analysis including Monte Carlo simulations encompassing ancillary parameters and rainfall forcing such as that described by Tao and Barros ( 2013) is out of the scope of this manuscript, it is important to characterize the elasticity of FS with regard to changes in key soil properties.Here the focus is on the physical basis of the initiation process, and thus we conduct a targeted sensitivity analysis focusing only on two critical parameters including soil internal friction angle and cohesion.Specifically, motivated by the cold-season events' results for the JCB and by previous work (e.g., Arnone et al., 2011;Lu and Godt, 2008, among others), the uncertainty in FS caused by the specification of the soil internal friction angle and cohesion parameters is examined in detail.Other parameters, such as soil properties, are not tested here.Recall that, due to the lack of site-specific measurements or estimates, the soil internal friction angle and cohesion used in Sect.4.2 and summarized in Table 1 were defined based on the representative values reported by Witt (2005b).
Physically reasonable ranges of 10  respectively.The FS profiles as a function of friction angle (Fig. 16) show that instability takes hold for friction angles below 20 • and for relatively shallow soil depths.The uncertainty in FS associated with the soil friction angle in the top layer is relatively smaller than for bottom layers.Changes of ±20 % in friction angle lead to similar change in FS magnitude for the 2nd to 4th soil layers.The FS is more sensitive in the bottom soil layers as indicated by the steeper slopes in the rightmost panels.In reality, the soil internal friction angle should vary horizontally and vertically to capture changes in soil texture and soil structure with depth, which are not considered explicitly in this study.In addition, soil heterogeneity, LULC change, bioactivity, and prior landslides can play an important role in determining effective soil internal friction angles locally.
Compared to the large uncertainty caused by changes in the soil friction angle, the changes in cohesion have a smaller impact on the magnitude of FS (Fig. 17) than in previous studies using other models (e.g., Wu and Sidle, 1995).Contrary to the sensitivity behavior with respect to friction angle discussed above, FS is more sensitive to changes in cohesion in the top layer.For instance, 50 % change in cohesion causes 22 % change in the magnitude of FS for the top layer, but just about 5 to 10 % change in the bottom layers for the case in the BCB.The uncertainty caused by cohesion for the coldseason event in the JCB is much smaller (Fig. 17b and c), within ranges of −15 to 10 %.Differentiating Eqs. ( 6) and ( 8) with respect to cohesion implies that the changes of FS actually depend on wet-soil-specific weight and depth.Regarding the role of root systems in forested catchments, note that the density decreases with depth, and more so between the third and base layers as specified in the model (Sect.3).

Conclusion and discussion
A fully distributed hydrologic model coupled with slope stability models was used to investigate the mechanisms triggering initiation of debris flow in two headwater catchments of the Pigeon River Basin in the Southern Appalachian Mountains, USA, for both warm-and cold-season debris flow events.The summer event took place during the passage of an intense (heavy rainfall intensity) nocturnal convective system.The winter events took place during a long-lasting stratiform (light to moderate rainfall intensity) orographic storm system, and a severe short-duration winter storm.Two slope stability models derived from the infinite slope model using the limit equilibrium method were utilized in this study: one is the modified slope stability index (SSI) calculated from soil wetness and slope neglecting cohesion and suction effects; the other is the factor of safety (FS) accounting for most of the dominant factors controlling slope instability.The SSI is based on temporal-spatial quantitative information about instability, but the subsequent aggregation of this information into threshold-based classes introduces ambiguity.For instance, pixels classified as unconditionally unstable are not necessarily always highly susceptible to slope failure.Sensitivity analysis of the FS estimates to soil strength parameters at rest, specifically the soil internal friction angle and cohesion due to soil and vegetation, was conducted.The results indicate that the FS exhibits strong sensitivity to friction angle, which increases with soil depth.The opposite occurs with respect to cohesion: sensitivity is modest, and it is significant only in the top soil layers reflecting the model's implicit representation of root systems.This is an important result, as anthropogenic activity, bioactivity, as well as prior slope movements, in addition to heterogeneities in soil structure and composition, can have a strong impact on the effective value of soil friction angle at small scales.The SSI approach tends to strongly overestimate the spatial distribution of slope instability due to its high sensitivity to instantaneous soil moisture at local scales, while neglecting soil cohesion and suction effects.Nevertheless, there is still ambiguity in the FS method in that it yields a large number of pixels with FS ≤ 1 + ε (ε is a measure of uncertainty for model estimates).Whereas the coupled modeling framework presented here does capture the locations of known debris flows, there are a number of locations where no debris flow was initiated and yet are nominally unstable.Clearly, not all factors determining initiation are included here, such as previous history of landslide activity, which should impact locally soil depth and structure.In addition, we hypothesize that there should be a scale effect associated with the spatial resolution of the model itself, and thus there should be practical utility in investigating the dependence of simulated soil moisture and interflow conditions at the time of landslide initiation on model resolution.Furthermore, the ability to represent heterogeneity and subgrid scale variability in subsurface flow dynamics should have a strong impact on the magnitude of interflow at small scales.
Although the debris flow initiation time with respect to the beginning of the storm differs for warm and cold seasons and from basin to basin, interflow magnitude controls the flow responses for all the events and is closely related to the trigger mechanism of shallow landslide initiation followed by debris flow mobilization.We demonstrated that for all the three case studies the interflow reaches the peak magnitude around the time when debris flows occurred at the initiation locations.That is, the timing of debris flow initiation is that when the interflow peaks independently of watershed, or storm type.We propose that the spatial ambiguity in FS prognostics can be addressed, at least in part, by monitoring the temporal evolution of interflow virtually using a modeling system such as described here.Finally, the methodology employed in this study fits in the same general framework for operational QFFs (quantitative flash flood forecasts) using quantitative precipitation estimates (QPEs) and quantitative precipitation forecasts (QPFs) described by Tao and Barros (2013).Thus, the prediction of debris flows could be made concurrently with QFF.The prediction of debris flow is necessary in order to issue timely warnings that can prevent or decrease loss of life and property especially downslope of debris flow initiation points, as well as to identify locations for forensic surveys in uninhabited areas, and where observing systems are not available.

Fig. 1 .
Fig. 1.Topography, major rivers and rain gauges over the Pigeon River Basin in North Carolina, USA.The Big Creek Basin (BCB) and the Jonathan Creek Basin (JCB) are marked, and the Cataloochee Creek Basin (CCB) used for hydrological verification is illustrated by shaded area.Simulated events of interest are marked using circles.The debris flow that occurred in the JCB on 7 Januray 2009 destroyed a house completely (shown in the picture below, courtesy goes to Dr. Richard Wooten).Land cover and soil texture are also provided, indicating the spatially varying vegetation and soil types over the basins.

Fig. 2 .
Fig. 2. Slope stability index classified by relationship between degree of soil saturation and slope, modified from Fig. 2 in Montgomery and Dietrich (1994).

Fig. 3 .
Fig. 3. Conceptual schema of the geotechnical system, explicitly showing the essential forces acted on a slope.

Fig. 4 .
Fig. 4. Comparison of hourly precipitation rate (mm h −1 ) between rain gauge observations and Q2 estimations before and after adjustment during 12 to 17 July 2011 for CCB (a) and BCB (b), during 5 to 10 January 2009 for CCB (c) and JCB (d), and during 6 to 11 December 2009 for CCB (e) and JCB (f).

Fig. 6 .
Fig. 6.The spatially varying soil depth estimated by two simple methods, and the ultimate soil depth used in this study averaging the two estimated soil depth.

Fig. 7 .
Fig. 7.The saturated hydraulic conductivity, soil porosity, field capacity, and wilting point extracted from the State Soil Geographic (STATSGO) database for four soil layers from the left to right, according to spatially varying soil depth.

Fig. 8 .
Fig. 8.The comparison between simulated streamflow at the outlet of the CCB, generated from the 3D-LSHM driven by Q2 rainfall datasets before and after adjustment for the event in July of 2011 (a), and in January (b) and December (c) of 2009; the flow components of estimated streamflow by adjusted Q2 including overland flow, interflow and baseflow are shown in the middle row (a.2, b.2 and c.2); and the interflow produced from three soil layers are shown in the bottom.The upper and right axis in figures indicate basin areal averaged storm hyetograph.

Fig. 9 .
Fig. 9.The time series of soil moisture (top) and interflow (bottom) produced at each soil layer at the pixel in which debris flow occurred.The x axis is zoomed into the rainfall period to show details more clear.The dash lines indicate the time when the magnitude of total interflow reaches its peak.

Fig. 10 .
Fig. 10.The spatial distribution of soil moisture, interflow for each soil layer and total interflow in the basins at the time when the debris flow occurred (indicated by dash line in Fig. 9) in July of 2011 (a), and in January (b) and December (c) of 2009.The debris flow locations are marked by circles.Slope and rainfall rate are also shown for reference.Channel pixels are not shown.

Fig. 11 .
Fig. 11.The histograms of soil moisture and interflow in each soil layer, slope, rainfall rate and total interflow generated using data from all over the basin for the entire simulation period, for the event in July of 2011 (a), and in January (b) and December (c) of 2009.The vertical red solid lines mark local conditions in the unstable grid element selected for analysis (corresponding to the gray solid line in upper interflow time series).In (b), both the red and gray dash lines indicate the condition when the largest rainfall rate took place at the pixel but the debris flow did not occur.

Fig. 12 .
Fig. 12. Scatter plots of the factor of safety (FS) versus slope and elevation for each grid element in the basins during the simulation for the event in July of 2011 (a), and in January (b) and December (c) of 2009.The circles are colored according to the magnitude of outgoing interflow.The pixel of interest is highlighted by a black circle.

Fig. 13 .
Fig. 13.Temporal evolution of the vertical profiles of soil moisture, pore pressure head, absolute interflow values and factor of safety at one of the debris flow initiation points in the basins in the event in July of 2011 (a), and in January (b) and December (c) of 2009.Color scheme indicates time.

Fig. 14 .
Fig. 14.The spatial distribution of slope stability characterized by the slope stability index (SSI, top), and the factor of safety (FS, bottom) at the time the debris flow occurred in the basins during the event in July of 2011 (a), and in January (b) and December (c) of 2009.The debris flow locations are marked by circles.The right-hand side panels are spatial zooms into the initiation zone.

Fig. 15 .
Fig. 15.The time series of the number of the total unstable pixels identified in the basins using the SSI (top) and the FS (bottom) metrics.Note that in this work, the debris flow proper is not simulated, and therefore the simulation is not representative of realistic conditions after debris flow initiation.For example, it is expected that, with mass movement, the shear stresses at locations surrounding the initiation points will decrease, and thus there should be a strong decrease in the number of unstable pixels.

Fig. 16 .
Fig. 16.Sensitivity analysis of the vertical profile of FS to soil internal friction angle, for the slope failure cases in the BCB in July of 2011 (a) and JCB in January of 2009 (b) and December of 2009 (c).The dark line is the actual failure case using the representative friction angle, 26 • .

Fig. 17 .
Fig. 17.Sensitivity analysis of the vertical profile of FS to the combined cohesion for soil and vegetation, for the slope failure cases in the BCB in July of 2011 (a) and JCB in January of 2009 (b) and December of 2009 (c).The dark line is the actual failure case using the representative cohesion, 2000Pa.

Table 1 .
Major parameters specified in LSHM for the three basins.

Table 2 .
Summary of debris flow sites' characteristics.The NCDOT Materials Testing Laboratory in Asheville, North Carolina, conducted soil quality tests on the soil samples from the debris flow initiation zones (data provided by NCGS, from Dr. Richard Wooten).
• -35 • for the soil internal friction angle and 500-3000 Pa for soil and vegetation cohesion were tested; the results are shown inFigs.16 and 17,