Quantification of the Influence of Preferential Flow on Slope Stability Using a Numerical Modelling Approach

The effect of preferential flow on the stability of landslides is studied through numerical simulation of two types of rainfall events on a hypothetical hillslope. A model is developed that consists of two parts. The first part is a model for combined saturated/unsaturated subsurface flow and is used to compute the spatial and temporal water pressure response to rainfall. Preferential flow is simulated with a dual-permeability continuum model consisting of a matrix domain coupled to a preferential flow domain. The second part is a soil mechanics model and is used to compute the spatial and temporal distribution of the local factor of safety based on the water pressure distribution computed with the subsurface flow model. Two types of rainfall events were considered: long-duration, low-intensity rainfall, and short-duration, high-intensity rainfall. The effect of preferential flow on slope stability is assessed through comparison of the failure area when subsurface flow is simulated with the dual-permeability model as compared to a single-permeability model (no preferential flow). For the low-intensity rainfall case, preferential flow has a positive effect on drainage of the hillslope resulting in a smaller failure area. For the high-intensity rainfall case, preferential flow has a negative effect on the slope stability as the majority of rainfall infiltrates into the preferential flow domain when rainfall intensity exceeds the infiltration capacity of the matrix domain, resulting in larger water pressure and a larger failure area.


Introduction
Landslides are commonly triggered by rainfall events.Hydrological models can be integrated with slope stability analysis methods to calculate the factor of safety and predict the time and magnitude of landslides (Crosta and Frattini, 2008;Shuin et al., 2012;Aleotti and Chowdhury, 1999;Westen et al., 2006).Combined hydro-mechanical models can roughly be divided into two types: simplified conceptual models (Montrasio and Valentino, 2008;Dai et al., 2002) and numerical models (Stead et al., 2001;Jing, 2003;Brinkgreve et al., 2010;Pastor et al., 2008), and have different levels of complexity depending on the scale and the research purpose.
The limit equilibrium method or infinite slope stability approach is frequently integrated with Richards' equation (Lanni et al., 2013;Ng and Shi, 1998;Godt et al., 2008;Shuin et al., 2012;Wilkinson et al., 2002;Talebi et al., 2008;Greco et al., 2013) or conceptual models (Arnone et al., 2011;Simoni et al., 2008;Qiu et al., 2007) for landslide hazard evaluation.The limit equilibrium method and infinite slope approach assume or search for a potential failure surface.The factor of safety is defined as the ratio between the maximum retaining force and the driving force (e.g.Lanni et al., 2013;Lu et al., 2012).Although the underlying assumptions of the slope failure mechanism have limitations (Huang and Jia, 2009;Griffiths et al., 2011), the simplified slope stability analysis method has low computational demand and is widely used for geotechnical analyses at the slope scale Published by Copernicus Publications on behalf of the European Geosciences Union.(Talebi et al., 2008;Tsai and Yang, 2006;Abramson, 2002), watershed scale, and catchment scale (Borga et al., 2002a, b;Baum et al., 2010;Wilkinson et al., 2002).
The strength reduction method (Griffiths and Lu, 2005;Huang and Jia, 2009) or local factor of safety method (Lu et al., 2012) can result in similar factor of safety values and locations of the critical slip surface as the limit equilibrium method, while no assumption is needed for the critical failure surface (Griffiths and Lu, 2005;Hammouri et al., 2008;Kim et al., 1999).The location, shape, and magnitude of the plastic deformation area are used to quantify the slip surface and factor of safety (Griffiths and Lane, 1999).Geotechnical engineering software and numerical models -such as FLAC ® (Itasca, 2002) or PLAXIS (Brinkgreve et al., 2010, based on the strength reduction method) -have been widely applied for slope stability analysis under the influence of transient hydrological conditions, such as rainstorms (Mukhlisin et al., 2008;Hamdhan and Schweiger, 2011) and reservoir water level variations (Huang and Jia, 2009;Zhou et al., 2014).
The Darcy-Richards equation is the most widely used approach in current software packages, but cannot effectively simulate preferential flow resulting in rapid infiltration (Nieber and Sidle, 2010;Beven and Germann, 2013).In highly heterogeneous slopes, preferential flow and transport can fundamentally influence subsurface flow (Jarvis, 2007;Hendrickx and Flury, 2001) and contaminant transport (Köhne et al., 2009;Allaire et al., 2009;Debieche et al., 2012;Zehe et al., 2001).A chain of connected macropores is commonly found in various types of soils, including forest soil and semiarid land (Uchida et al., 2001;Jarvis, 2007;Flury et al., 1994).For example, an earthworm burrow can extend from the surface deep into the soil, as can decayed plant roots or soil cracks (Jarvis, 2007;Beven and Germann, 1982;Hendrickx and Flury, 2001).The self-organising preferential flow network will become active and hydraulically connected with an increase in soil saturation (Nieber and Sidle, 2010).The saturated hydraulic conductivity of preferential flow paths is significantly larger than that of the soil matrix (Beven and Germann, 1982;Köhne et al., 2009).A significant portion of subsurface stormflow (Uchida et al., 2004;Zhang et al., 2006;Beven, 1981) is transmitted via preferential flow paths (Nieber and Sidle, 2010).Preferential flow through macropores, fractures, and other local highpermeability zones is extremely rapid, and contributes instantly to high porewater pressures in deep soils (Jarvis, 2007).
Quantification of landslide triggering mechanisms is an essential step in landslide forecasting.Some studies have shown that preferential flow is one of the major mechanisms affecting the timing and location of landslides (Sharma and Nakagawa, 2010;Uchida, 2004;Verachtert et al., 2013).In hillslopes, preferential flow paths, such as macropores, soil pipes, and fissures, have been associated with slope stability (Hencher, 2010;McDonnell, 1990;Uchida et al., 2001;Krzeminska et al., 2012;Debieche et al., 2012).Besides the fact that internal erosion in preferential flow paths deteriorates the slope mass and reduces the soil shear strength, the occurrence of preferential flow can give rapid access to the deeper soil and groundwater system, reduce soil shear strength (due to pore pressure changes), and influence the timing and frequency of landslides (Wienhöfer et al., 2011;Uchida, 2004;Verachtert et al., 2013).
Preferential flow and solute transport have been simulated at various scales including the scales of pores, soil columns, hillslopes, and catchments (Šimůnek et al., 2003;Gerke, 2006;Köhne et al., 2009) using increasingly sophisticated models such as the dual-porosity/dual-permeability model (Gerke and van Genuchten, 1993a;Jarvis et al., 1991;Larsbo and Jarvis, 2003), the multi-permeability model (Wu et al., 2004;Greco, 2002;Gwo et al., 1995), and conceptual models (Armstrong et al., 2000;Weiler, 2005;Mulungu et al., 2005).The dual-permeability model is widely used because of its clear physical concept and powerful simulating ability (Roulier and Jarvis, 2003;Kodešová et al., 2005;Gerke and Köhne, 2004;Köhne et al., 2006;Christiansen et al., 2004;Weiler, 2005;Therrien and Sudicky, 2005;Vogel et al., 2000).The dual-permeability model assumes that the soil consists of two interacting, overlapping pore domains.The matrix domain with relatively low permeability represents the soil micropores where flow is governed by Richards' equation.The preferential flow domain represents the highly permeable preferential flow paths, such as macropores, fractures, cracks, or large pores between soil aggregate.Preferential flow is described by Richards' equation (Šimůnek et al., 2008;Gerke and van Genuchten, 1993a) or the gravity-driven kinematic wave equation (Larsbo and Jarvis, 2003;Jarvis et al., 1991;Greco, 2002).The water exchange between the two domains is driven by the pressure head difference between the two domains (Pirastru and Niedda, 2010;Gerke and van Genuchten, 1993b).Dual-permeability models have proven to be effective for preferential flow simulation, but have not been incorporated into slope stability models.
The objective of this study is to quantify the temporal and the spatial effect of preferential flow on slope stability, and to analyse its underlying hydrological mechanisms using numerical experiments of rainfall-induced shallow landslides.This paper is organised as follows.The subsurface dualpermeability hydrological model is described.The subsurface hydrological model is sequentially coupled with a soil mechanics model and a stress-field-based local factor of safety slope stability method (Sect.2.2).The numerical experiments and parameterisation are discussed in Sect.3. The hydrological and geotechnical results are given in Sect. 4. The influence of preferential flow on subsurface hydrological processes and consequent slope stability is discussed in Sect. 5 by comparing the results of single-and dualpermeability models.

Subsurface flow model
The single-permeability model is described by one Richards' equation to represent flow in a homogenous soil.The dual-permeability model divides the flow domain into two overlapping and interacting continua, where two coupled Richards' equations are used to describe the matrix flow and preferential flow (Gerke and van Genuchten, 1993a): where the subscript f indicates the preferential flow domain and the subscript "m" indicates the matrix domain.C is the differential water capacity (dθ/dh) (L −1 ), is the effective saturation (-), h is the pressure head (L), t is time (T), z is the vertical coordinate (positive upward), K is the isotropic hydraulic conductivity (L T −1 ), S s is the specific storage (L −1 ), w is the volumetric ratio of the preferential flow domain or the matrix domain over the total soil volume (-), and w is the water exchange term (T −1 ) between the two domains.The Brooks-Corey function is used to describe the hydraulic properties of both the matrix and preferential flow domains (Brooks and Corey, 1964) where θ is the water content (L 3 L −3 ), subscripts s and r denote saturation and residual state, K s is the saturated hydraulic conductivity (L T −1 ), and α BC , l BC , n BC , are fitting parameters.
w is the water exchange rate between the two domains (Ray et al., 1997): where α w (L −2 ) is the effective water transfer coefficient, and the relative hydraulic conductivity K a (L T −1 ) is calculated by averaging the hydraulic conductivities of the two pore domains (Arora et al., 2011;Laine-Kaulio et al., 2014): The volumetric ratio of the preferential flow domain and matrix domain sum up to 1: The total water content of the soil is the weighted average of the water contents of the two domains: The same holds for the total saturated hydraulic conductivity of the soil: where K sf and K sm are saturated hydraulic conductivity of preferential flow domain and matrix domain, respectively.Boundary conditions may be specified for pressure head, specified flux, or mixed (Chui and Freyberg, 2009).In the case of a dual-permeability model, specified flux i (infiltration from rainfall) is divided between the matrix and preferential flow domains: where i m and i f are boundary fluxes to the matrix and the preferential flow domains (L T −1 ), respectively.The two domains have an equal opportunity to receive rainfall and are initially equal to rainfall intensity R (Dusek et al., 2008): As the matrix domain has a larger volumetric ratio (w m > w f ), the infiltration process is initially dominated by the matrix domain.The infiltration capacity of each domain is the product of saturated hydraulic conductivity and the pressure head gradient of each domain.Once the specified flux into the matrix is larger than its infiltration capacity, the boundary condition changes to specified pressure head and the specified flux for the preferential flow domain is increased to Once the specified flux into the preferential flow domain is also larger than the infiltration capacity, the boundary conditions of both the matrix and the preferential flow domain are changed to a specified pressure head of zero and overland flow occurs.

Slope stability analysis method
The slope stability analysis is based on the local factor of safety approach (Lu et al., 2012).The plane-strain linear elasticity model is used to calculate the 2-D stress field (Abramson, 2002), which is governed by a momentum balance equation where σ is a stress tensor (M L −1 T −2 ) with three independent stress variables in two-dimensional space, γ is the bulk unit weight of the slope material (ML −2 T −2 ), and b is the unit vector of body forces with two components.The effective stress equation can be formulated as where σ is the effective stress, p w (M L −1 T −2 ) is the porewater pressure, and χ (-) is the matrix suction coefficient, which is usually approximated by the effective saturation (Lu et al., 2010).
The local factor of safety F LFS is defined as the "ratio of the Coulomb stress at the current state of stress to the Coulomb stress of the potential failure state under the Mohr-Coulomb criterion" (Lu et al., 2012): where τ * is the limit Coulomb stress and τ is the actual shear stress (M L −1 T −2 ).Application of the Mohr-Coulomb failure criterion gives the local factor of safety at every point in the hillslope: where c is the effective cohesion (M L −1 T −2 ), φ is the friction angle, and σ 1 and σ 3 are the spatially varying first and third principal effective stress for the variably saturated soil The influence of hydrology on slope stability is manifested in two ways.First, the unit weight function depends on the water content (Eq.9).Second, the effective stress depends on the porewater pressure.In the dual-permeability model, the porewater pressure of the preferential flow domain is used in the computation of the effective stress.
Figure 1 summarises the structure of coupled dualpermeability and slope stability model.Two Richards' equations are coupled by the water exchange function.The hydrological results are sequentially coupled with a soil mechanics model without considering possible feedback of soil deformation on soil properties and the hydrological process.
In this study, COMSOL Multiphysics ® is used to develop a coupled hydrological and slope stability model for both single-permeability and dual-permeability subsurface flow (Shao et al., 2014).For subsurface flow computations, COMSOL Multiphysics facilitates variable mesh resolution to optimise between coarse and fine mesh density to balance 3 Set-up of the numerical experiments

Slope geometry
Consider a slope of 23 • consisting of fine-grained lithology such as clay shales with a more permeable weathered top soil layer (Bogaard, 2002;Berti and Simoni, 2012;Picarelli et al., 2006); this is a typical slope that is vulnerable to failure.The slope is 6 m high and 15 m long and consists of two layers with a 2 m thick homogeneous upper soil layer (see Fig. 2).
The model domain is 42 m by 25 m to reduce the effects of specified boundary conditions along the model boundary.The computational mesh and the boundary conditions are shown in Fig. 2. The boundary conditions of the subsurface flow model are atmospheric at the surface: the left-hand and bottom sides are no-flux boundaries; the right-hand side consists of a seepage boundary condition for the upper soil layer and a specified pressure head to mimic a constant groundwater table for the lower layer.For the soil mechanics model, the surface is a free boundary, the bottom boundary (only horizontal displacements) and the left-and right-hand sides (only vertical displacements) are all roller boundaries.
Since the pressure head in the surface area can change drastically during rainfall, a very dense mesh was used near the surface to accurately model the transient hydrological conditions.The mesh density of the upper layer is approximately 0.25 m (vertical) by 0.5 m (horizontal).A coarser mesh was defined in the lower part of the slope as a less dynamic condition will occur here.

Parameterisation
The volumetric ratio of the preferential flow domain w f is 0.1; a typical range is ∼ 0.025-0.2(Köhne et al., 2002) pore-size distribution of the preferential flow domain allows unsaturated infiltration before the matrix domain is saturated (Dusek et al., 2008).A comparison is made between the hydrological results of the single-permeability and the dualpermeability models.The total weighted saturated hydraulic conductivity of the dual-permeability model is equal to the saturated hydraulic conductivity of the single-permeability model.The water exchange between the matrix and preferential flow domains depends on the hydraulic conductivity between the two domains K a and the water exchange coefficient α w (Eq. 6).Equilibrium between the preferential flow and matrix domains is reached quicker for smaller values (closer to 1) of K sf /K sm and larger values of α w .Moderate values are used for K sf /K sm (100 in the upper layer and 5 in the lower layer) and for α w (0.2 m −2 ).The soil hydraulic parameters are presented in Table 1.Preferential flow plays an important role in the upper soil layer where there is an abundance of macropores, but less so in the lower soil layer where macropores are almost nonexistent (e.g.Bogner et al., 2013).In other words, the volumetric percentage of preferential flow domain is still the same, but in the lower layer the saturated hydraulic conductivity of macropores is more similar to the pores of the matrix.The more permeable top layer is sandy loam and the fine-grained lower layer is clay; the soil hydraulic parameters are taken from the unsaturated soil hydraulic database (UNSODA) (Nemes et al., 2001;Leij, 1996).
Parameterisation of soil hydraulic properties for dualpermeability models is difficult as the two conceptualised domains cannot be experimentally separated.In this paper, the pragmatic approach of Vogel et al. (2000) is adopted: the same hydraulic parameters are adopted for both domains, except for the saturated hydraulic conductivities.
The parameters of the soil mechanics model are also shown in Table 1.The slope stability is a function of two parameters: the friction angle and the effective cohesion.Decreasing the friction angle or increasing the effective cohesion leads to a higher factor of safety.The friction angle is physically constrained in a narrow range, and in this study is fixed.In numerical modelling, effective cohesion c is scale dependent, and is usually defined as a linear function of the slope height to obtain identical values of the safety factor when applying it to different slope sizes (Griffiths and Lane, 1999;Lu et al., 2012).Very low values of effective cohesion result in initially unstable slopes, and very high values of effective cohesion results in slopes that are unconditionally stable.In this study, two sets of cohesion values were selected in the range that may result in rainfall-triggered landslides: a homogeneous case where the effective cohesion of both layers is 5 kPa and a case where the effective cohesion of the upper layer is smaller (c 1 = 3 kPa) than the lower layer (c 2 = 6 kPa).
Intermittent and variable rainfall may significantly influence the pressure response in the shallow layer soil, while the pressure response in the deeper layer soil is driven by percolation and less dynamic to intermittent rainfall.Two rainfall events are modelled: a low-intensity rainfall of 2 mm h −1 for 150 h and a high-intensity rainfall of 20 mm h −1 for 15 h.The two rainfall events are expected to result in distinctly different slope stabilities, even though the total rainfall of both events is equal.The initial condition is the steady porewater pressure distribution obtained from running the model with a daily rainfall of 1.64 mm day −1 (600 mm year −1 ) for 10 years.

Subsurface flow
A schematic diagram of the subsurface flow components in the study area is shown in Fig. 3.Note that the study area is a small part of the model domain (Fig. 2).The main fluxes are the infiltration from rainfall (blue), the inflow/outflow along the left side and bottom (black), the seepage outflow along the surface (red), and the outflow along the right boundary (green).
Hydrological results for the single-and dual-permeability models are shown in Figs. 4 and 5, respectively.The graphs on the left are results for the long-duration, low-intensity rainfall case while the graphs on the right are results for the short-duration, high-intensity rainfall.Integrated fluxes, as shown in Fig. 3, are reported in m 2 h −1 .
For both models, all the rainfall infiltrates into the slope during the beginning of the rain event and when the soil becomes saturated infiltration decreases and saturation excess overland flow occurs.For the single-permeability model and low-intensity rainfall overland flow starts after 95 h (or 190 mm of rainfall) while for the high-intensity rainfall overland flow starts after 8.5 h (or 170 mm of rainfall) (Fig. 4a  and b).
In the dual-permeability model (Fig. 5), the rainfall infiltration is divided over the two domains and additional rainfall infiltrates into the preferential flow domain when the matrix domain reaches infiltration capacity.Recall that the matrix domain is 90 % of the domain, and the preferential flow domain is 10 % of the domain.A smaller fraction of rainfall infiltrates into the preferential flow domain for the case of low-intensity rainfall (10-15 %) than for the case of highintensity rainfall (50-85 %).Overland flow starts after 80 h (or 160 mm of rainfall) for the low-intensity case and after 60 h (or 120 mm of rainfall) for the high-intensity case.
The seepage outflow increases along the left, right, and bottom boundaries during the rainfall event (Figs.4c, d,  and 5c, d) and is smaller than the infiltration rate (storage is increasing).In the dual-permeability model and the lowintensity rain, outflow along the surface boundary starts after 115 h (or 230 mm of rainfall) while for the high-intensity rain outflow starts after 9 h (or 180 mm rainfall).The outflow rate along the surface boundary depends on the groundwater level in the upper layer.In the dual-permeability model, the outflow along the right boundary is approximately 10 times larger for the preferential flow domain than for the matrix domain, which is consistent with their volumetric ratio and their saturated hydraulic conductivity ratio.The water exchange between the two domains in the dual-permeability model is shown in Fig. 5e and f.For the low-intensity rainfall case, the water exchange from the preferential flow domain to the matrix domain increases during the first 100 h and then decreases, while the water exchange from the matrix domain to the preferential flow domain is almost always increasing (more negative).For the high-intensity rainfall case, the water exchange from the matrix to the preferential flow domain is negligible, while the water exchange from the preferential flow domain to the matrix domain reaches more than 0.3 m 2 h −1 , which is similar to the infiltration into the preferential flow domain.After 5 h, approximately 75 % of infiltration into the matrix domain is water exchange from the preferential flow domain (Fig. 5f) and 25 % infiltration from the surface boundary (Fig. 5b).

Water balance
The integrated rainfall and water storage for the study area are shown for both models in Fig. 6.The water balance is obtained by integrating all flow components along the boundaries of the study area.The numerical water balance errors are between 2 and 3 %.
For all cases, the storage increase flattens out when the inflow decreases (Figs. 4 and 5).For the high-intensity rainfall, the dual-permeability model stores 8 % less water than the single-permeability model.The total storage after 150 h of low-intensity rainfall is less than after 15 h of high-intensity rainfall, probably caused by the longer time that water can drain from the study area under low-intensity rain.
For the dual-permeability model, the water exchange has a significant influence on the storage change in each domain.For the low-intensity rainfall, the storage in the preferential flow domain does not increase much after 6 h (Fig. 6).For the high-intensity rainfall, the storage in the preferential flow domain increases rapidly over the first 3 h as very little water infiltrates into the matrix domain due to the low infiltration capacity of the matrix.After 3 h, the preferential flow domain has almost reached full saturation and the large pressure difference between the preferential flow domain and matrix domain causes extensive water exchange (Fig. 5f).

Water content
The water content distribution in the study area is shown in Fig. 7 for both the single-permeability model (left-hand panels) and the dual-permeability model (centre and righthand panels).The water exchange rate between the matrix and preferential flow domains of the dual-permeability  model is shown in Fig. 8.The infiltration process of the dual-permeability model differs significantly from that of the single-permeability model.The initial water content distribution in the matrix, as well as preferential flow domains, is similar for both models.During the rainfall events, the wetting front in the singlepermeability model develops parallel to the surface and propagates downward.This holds for both low and high rainfall intensities (Fig. 7 left-hand column).The wetting front generally reaches the groundwater table at the toe of the slope first, after which the infiltrated water continuously enlarges the saturated area.
In the dual-permeability model, the combined effects of the preferential flow and the matrix flow show a more complicated response.For the low-intensity rainfall, infiltration is dominated by matrix flow, as 90 % of the subsurface consists of the matrix.Because the rainfall intensity is lower than the saturated conductivity of the matrix domain, on the surface boundary, 90 % of the rainfall infiltrates into the matrix domain and 10 % of the rainfall infiltrates into the preferential flow domain.The pressure difference between the two domains drives the water exchange at the matrix wetting front (Figs.5e and 8a).Inside the flow domain of the slope, 10 % of matrix flow transferred to preferential flow due to water  exchange in the soil.At first, water quickly reaches the soil layer interface by preferential flow where it transmits to the matrix, although this exchange flux is very small (Figs.5e  and 8a).After sufficient time (70 h), a much stronger matrix flow (taking about 80 % of the infiltrated rainfall) reaches the soil layer interface and generally reverses the water exchange direction (Fig. 5e).Overall, water exchange during low-intensity rainfall in the study area is dominated by flow from the matrix to the preferential flow domain (Fig. 8a and  b).
For the high-intensity rainfall, the rainfall intensity is 8.4 times the matrix saturated hydraulic conductivity.The percentage of infiltration into the matrix domain decreases from 90 to 50 % within the first half hour, and continues to decrease to less than 20 % after 1.5 h.In contrast, the percentage of rainfall that infiltrates into the preferential flow domain increases from 10 to over 80 % after 2 h.Water in the preferential flow domain quickly reaches the deeper soil layer and forms a perched groundwater table (Fig. 7), where a significant amount of water infiltrates into the matrix (Fig. 5f).

Slope stability
The local factor of safety is computed based on the computed water-pressure distribution (Fig. 7).The distribution of the local factor of safety is shown in Fig. 9 for the initial condition and after 150 h (low-intensity rainfall) and 15 h (highintensity rainfall) for both the single-permeability model and the dual-permeability model and for the case with different cohesion values for the upper and lower layers.The case with equal cohesion values is not shown because the potential failure areas are very small.
A local factor of safety below 1 indicates a potential failure area.The area with a F LFS below 1 was determined every time interval (5 h in case 1, and 0.5 h in case 2) and is shown by the black line in Fig. 9. Slope stability is related to both the specific weight of the wet soil and the porewater pressure in the soil.The specific weight changes due to changes in water storage are relatively small, but changes in water pressure have a significant effect on slope stability, especially in the area of the perched water table .The size of the potential failure area is plotted vs. the cumulative rainfall in Fig. 10 for the two different rainfall events and two sets of cohesion values.The results for the same cohesion values (c 1 = c 2 = 5 kPa) are shown in Fig. 10a.For the low-intensity rainfall, the failure area is very small and is approximately the same for both permeability models.For the high-intensity rainfall, the failure area in the single-permeability model is larger than for the low-intensity rainfall, but the trend is similar.The failure area in the dualpermeability model is significantly larger.Failure starts after 60 mm of rainfall, and the failure area continues to grow during the rainfall infiltration process.
The results for different cohesion values (c 1 = 3 kPa, c 2 = 6 kPa) are shown in Fig. 10b.For the low-intensity rainfall, the failure area is 0.7 m 2 in the single-permeability model after 20 mm of cumulative rainfall.The size of this area shows almost no increase until approximately 220 mm of cumulative rainfall, when the groundwater table starts to rise (Fig. 7).The failure area of the dual-permeability model is 40 % smaller than that of the single-permeability model as the preferential flow domain drains more water into the ma- trix domain.For the high-intensity rainfall, the failure area of the dual-permeability model is larger than of the singlepermeability model, as is the case with equal cohesion values.The failure areas of both models increase fairly quickly to 2 m 2 , or 5 % of the upper layer in the study area.The failure area increases to 5 m 2 in the dual-permeability model and to 3 m 2 in the single-permeability model.
The slope stability result are directly related with subsurface hydrological results.For the low-intensity rainfall, the failure area for the single-permeability model is very similar in size and location to the dual-permeability model as the location of the water table is very similar in both models (Fig. 7).The initial condition of the dual-permeability model is slightly more stable than that of the single-permeability model, since the preferential flow domain has a higher drainage capacity and, consequently, a lower porewater pressure.In the case of low-intensity rainfall, the matrix flow dominates the groundwater recharge and, consequently, the slope instability.Furthermore, the porewater pressure in the preferential flow domain is very low due to its strong drainage capacity.As a result, the failure area calculated by the dual-permeability model under low-intensity rainfall is slightly smaller than that calculated by the singlepermeability model (Fig. 10a).The location of the failure area is similar in the single-and the dual-permeability domain (Fig. 9).
For the high-intensity rainfall, the failure area is significantly larger for the dual-permeability model than for the single-permeability model as the perched water table in the preferential flow domain is much more extensive in the dualpermeability model as compared to the single-permeability model (Fig. 7).The regular wetting front of the singlepermeability model does not reach the interface between soil layers, and the failure area is limited to the toe of the slope.For the dual-permeability model, the high-intensity rainfall results in a rapid infiltration through preferential flow, which quickly reaches the interface between soil layers, and increases the degree of saturation and pressure head of the deeper soil.Positive porewater pressure occurs in the preferential flow domain before the entire slope is fully saturated, and produces a larger failure area than in the equivalent single-permeability model.

Discussion
The role of preferential flow in hydrology focuses mainly on the rapid vertical infiltration of water and contaminant (Christiansen et al., 2004;Kodešová et al., 2005;Laine-Kaulio et al., 2014), or the rapid discharge in hillslope and catchment hydrological studies on discharge generation (Zhang et al., 2006;Mulungu et al., 2005).Preferential flow can cause rapid water-pressure build-up and may trigger landslides, as suggested by many researchers, but it is rarely quantified.Prior to this study, a detailed evaluation of the influence of preferential flow on slope stability using a coupled dual-permeability and slope stability model had not been carried out.In this study, a physically based numerical model was used to investigate the hydrologic response under predefined conditions.Our starting point was an idealised, but representative, 2-D slope with fixed parameters and boundary conditions (except rainfall forcing).A linear-elastic model was used to calculate the stress field to quantify the failure area and timing, but plastic deformation after failure was not considered; the model is not able to quantify landslide displacements.In this section, the underlying approximations of the numerical model are explored and the influence of the chosen parameter sets on the model outcome is discussed.The numerical experimental results are compared with field studies and other published numerical experiments.

Continuum model
Soil heterogeneity is one of the most difficult problems in both hydrology and soil mechanics studies.As an alternative to the continuum approach used here, preferential flow may be simulated by explicitly including fissures, pipes, or fracture networks in discrete (or discontinuous) models.Several field studies (Hencher, 2010;Uchida, 2004;Verachtert et al., 2013) and numerical experiments (Tsutsumi and Fujita, 2008;Chang et al., 2014) have focused on the investigation and simulation of pipe flow (in soil) and fracture flow (in rock).In order to accurately describe the geometry of the preferential flow paths, the high-resolution macropore image reconstruction approach (Hu et al., 2014) or the statistical approach (Köhne et al., 2009) may be applied.Numerical simulations of these natural macropore networks require large amounts of geometry information (Nieber and Sidle, 2010) and computational time and are consequently limited to small-scale studies with a limited number of pipes (Tsutsumi and Fujita, 2008) or cracks (Moonen et al., 2008).
The dual-permeability model is a useful tool to simulate subsurface stormflow and solute transport in a hillslope when the parameterisation is able to capture the hydraulic characteristics of each domain (Laine-Kaulio, 2011;Laine-Kaulio et al., 2014).As the dual-permeability model describes the subsurface as a continuum of two linked domains, it is suitable for heterogeneous slopes with a high density of preferential flow paths and not for slopes with only a few large fissures or cracks (van der Spek et al., 2013).
In this paper, flow in both domains is described with the Darcy-Richards equation, which is valid when the macropores have a relatively small size, and the macropore flow is still viscous (Köhne and Mohanty, 2005;Laine-Kaulio et al., 2014).When fluid velocities are high and flow becomes turbulent, Darcy's equation is not valid (Beven and Germann, 2013) as may be the case in large cracks or fissures under near-saturated or ponded infiltration (Beven and Germann, 1982).The existence of pore necks and dead ends in pref- erential flow paths reduce the occurrence of turbulent flow (Jarvis, 2007).

Coupling term in the dual-permeability model
In the dual-permeability model, the two domains are in general not at equilibrium.The parameterisation of the water exchange between the two domains governs whether a dualpermeability model behaves as a dual-permeability model or mimics the behaviour of a single-permeability model.The water exchange is governed by the pressure difference between the two domains and two parameters: the water exchange coefficient and the average hydraulic conductivity between the two domains (Eq. 6).The average hydraulic conductivity in turn is a function of the hydraulic conductivities of the two domains, which are a function of the pressure head.The larger the product, the quicker the two domains equilibrate.Specifically, the sensitivity analysis of Ray et al. (1997) shows that a very large water exchange coefficient results in almost instantaneous equilibrium between the two domains.Consequently, the dual-permeability model behaves like a single-permeability model, which means that the preferential flow has no influence on the hydrological and soil mechanics results.This behaviour does not agree with experimental studies of natural soils, which often show a non-equilibrium phenomenon between the two domains (Šimůnek et al., 2003;Jarvis, 2007).The parameterisation of the water exchange used in this study was compared with previous studies (see Table 2).Estimation of the water exchange coefficient from physical measurements is very difficult.The most widely used equation is (Gerke and van Genuchten, 1993b) where β is a scaling factor, d is half the representative distance between two macropores, and γ w is a geometrydependent shape factor that equals 3 for rectangular slabs and 15 for spheres (Ray et al., 1997).Parameter values for the water exchange term used in several studies are summarised in Table 2. Vogel et al. (2000) and Gerke and Köhne (2004) conceptualise the preferential flow domain as rectangular matrix blocks arranged as parallel slabs.A reduction factor of 0.01 or 0.001 was used to significantly reduce the water exchange between the two domains, because the hydraulic conductivity at the matrix/fracture interface was conceptualised to be controlled by relatively impermeable coatings that are composed of minerals and organic matter (Ray et al., 2004;Gerke and Köhne, 2002).Köhne and Mohanty (2005) conceptualise the dual domain as a hollow cylindrical matrix that is filled with coarse sand in the middle to mimic the macropore domain.Arora et al. (2011) based their parameters on a high density of macropore columns, and they calculated K a by averaging the hydraulic conductivities of the two pore domains (as adopted in this paper; see Eq. 7).Arora et al. (2011) and Köhne and Mohanty (2005) did not consider the influence of coatings on the permeability, nor was this done in this study.
It may be seen from Table 2 that the magnitude of the product α w K sa is similar for all studies, even though some of the other values (notably the ratios K sa /K sm and the values of α w ) differ by several orders of magnitude.As such, the water exchange between all these models is likely similar.

Computation of effective stress
In the dual-permeability model, the porewater pressure of the matrix and the preferential flow domains are different and water flows from the domain with a higher pressure to the domain with a lower pressure.Van der Spek et al. (2013) showed that in the case of varved clays with a low hydraulic conductivity of the soil matrix and a low density of fissures, the time delay between water entering the fissure network and an increase in pressure in the matrix is relatively large.This study concerns a system with a very high density of macropores and consequently the numerical simulations show only a small time delay for the pressure propagation from the preferential flow domain to the matrix domain.The porewater pressure of the preferential flow domain is used for the effective stress calculation in the slope stability analysis, but failure time and area are only slightly different when the matrix porewater pressure is used for the slope stability analysis.Field evidence (Uchida et al., 2001;Hencher, 2010;Wienhöfer et al., 2011) and numerical experiments (Nieber and Sidle, 2010;Laine-Kaulio et al., 2014) suggest that individual preferential flow networks are hydraulically connected, and that the high porewater pressure build-up in the preferential flow paths is directly correlated with slope stability.

Implications of preferential flow for slope stability
This study is not the first to address the influence of preferential flow on subsurface flow and slope stability.Preferential flow has an effect on infiltration and drainage fluxes and as such influences the triggering factors for rainfallinduced landslides.Moreover, storage capacity relates to the pore distribution in a soil and controls the antecedent condition or the cause of landslide occurrence (Hamdhan and Schweiger, 2011).The complexity hides in the combination of rainfall characteristics and soil hydraulic properties, together with the physiographic properties like slope, soil thickness, bedrock topography and so on, which determine the resultant porewater pressure response.The model runs and analyses show that rainfall intensity needs to be related to both the soil infiltration rate of the matrix domain and the preferential flow domain.Natural hillslopes show a different response from the matrix as compared to the preferential flow domain depending on the rainfall intensity.The dualpermeability model can simulate the preferential flow which will rapidly reach deep soil and increase porewater pressure.However, a single-permeability model cannot reproduce the phenomenon of preferential flow dominated fast pressure response in deep soil.Parameterisation of a dual-permeability model is difficult in practice (Laine-Kaulio et al., 2014).Therefore, the use of single-permeability models with effective soil hydraulic parameters prevails in regional hazard assessment (Hamdhan and Schweiger, 2011;Zhou et al., 2014).Rainfall-intensity duration plots for regional hazard assessment are well established and abundantly used but do not include soil and hydrological information (Guzzetti et al., 2007(Guzzetti et al., , 2008)).They empirically relate precipitation intensity and duration to observed landslides.The inclusion of more detailed hydrometeorological information in these analyses is ongoing.Recently, von Ruette et al. (2014) showed the importance of spatially and temporally heterogeneous rainfall on the initiation of landslides.In a synthetic study they showed that spatially distributed rainfall resulted in an increase of the number of shallow landslides as compared to uniform or intermittent rainfall (short periods with higher rainfall intensities but spatially homogeneous).They concluded that "low-rainfall intensities (below soil infiltration capacity) and long durations resulted in more infiltration, lower stream discharge, and more saturations and thus failure".This is in full agreement with the results for low rainfall intensities in this study.
Generally speaking, this holds for every case where infiltration capacity of the matrix remains higher than the rainfall intensity even in the presence of preferential flow paths.For low-intensity rainfall, the water-pressure increase simulated with a single-permeability model is generally larger than with a dual-permeability model as drainage by the preferential flow paths is underestimated.Soil drainage is a typical threshold process of the soil to get rid of its high porewater pressure and in this way stabilises the slope.Consequently, the stability is slightly underestimated with a singlepermeability model for low-intensity rainfall.
The reverse is true, however, for high rainfall intensities, when the matrix reaches infiltration capacity early on.In these cases the preferential flow system dominates because water that cannot infiltrate into the matrix domain infiltrates into the preferential flow domain instead, resulting in a large pressure increase with a negative effect on slope stability.A much smaller pressure increase is simulated with a single-permeability model for the same highintensity rainfall.Consequently, the stability is overestimated with a single-permeability model even when equivalent parameters are used.

Conclusions
A coupled dual-permeability and slope stability model was developed to simulate the influence of preferential flow on subsurface hydrology and consequent slope failure area.The dual-permeability model is able to simulate both preferential flow and matrix flow.The slope failure area was determined with a local factor of safety analysis.Numerical experiments were carried out to study the effect of rainfall events on slope stability with both a single-permeability (no preferential flow) and a dual-permeability model.A 23 • slope consisting of two soil layers was used in the study.The upper layer is sandy loam and the bottom layer is clay.Both the case where the cohesion of the two layers are equal, and the case where the cohesion of the upper layer is smaller than the lower layer were simulated.Two types of rainfall events were considered: low-intensity, long-duration rainfall and high-intensity short-duration rainfall.The total amount of water of both rainfall events was equal.The effect of preferential flow on slope stability was studied by comparing the failure area obtained with a single-permeability model and a dual-permeability model for the same rainfall event.
For low-intensity rainfall, the failure area of both models is similar when the cohesion of the upper and lower layers is equal, but the failure area is significantly larger in the singlepermeability model as compared to the dual-permeability model when the cohesion of the upper layer is lower than the cohesion of the lower layer.During low-intensity rainfall, preferential flow has a positive effect on slope stability as it drains water from the matrix domain and decreases the water pressure.For high-intensity rainfall, the failure area of the dualpermeability model is significantly larger than the singlepermeability model whether the cohesion values of the two layers are equal or not.During high-intensity rainfall, the rainfall intensity is larger than the infiltration capacity of the matrix domain so that most of the rainfall infiltrates into the preferential flow domain.As a result, the water pressure increases very quickly in the preferential flow domain resulting in a much larger failure area than is the case for the singlepermeability model.
In summary, two different effects of preferential flow on slope stability were identified with a coupled dualpermeability and soil mechanics model.Preferential flow has a positive effect on slope stability during low-intensity rainfall and a negative effect on slope stability during highintensity rainfall.The magnitude of the effect is a function of the soil hydraulic properties and soil mechanical properties of a specific slope.Identification of parameter ranges for which this behaviour is significant requires further investigation.

Figure 1 .
Figure 1.Structure of coupled dual-permeability model and soil mechanics model.

Figure 2 .
Figure 2. Computational mesh and boundary conditions.

Figure 3 .
Figure 3. Flow component and water balance of study area.

Figure 4 .
Figure 4. Integrated fluxes for single-permeability model and 2 mm h −1 (left panels) and 20 mm h −1 (right panels) rainfall.Rainfall and infiltration (a, b), and outflow at the right, outflow at the left and bottom, and outflow at the surface (c, d).

Figure 5 .
Figure 5. Integrated fluxes for dual-permeability model and 2 mm h −1 (left panels) and 20 mm h −1 (right panels) rainfall.Rainfall and infiltration (a, b); outflow at the right, outflow at the left and bottom, and outflow at the surface (c, d); exchange between matrix domain (MT) and preferential flow domain (PF) (e, f), positive for flow from PF to MT and negative for flow from MT to PF.

Figure 6 .
Figure 6.Storage increase of single-permeability model and dual-permeability model.

Figure 8 .
Figure 8. Water exchange rate distribution.Positive values (red): mean water exchange from preferential flow domain to matrix; negative values (blue): mean water exchange from matrix to preferential flow domain.

Figure 10 .
Figure 10.Development of the failure area under different rainfall intensities and soil cohesion.

Table 1 .
Summary of parameters.

Table 2 .
Parameters setting of water exchange coefficients in different literature.