A steady-state saturation model to determine the subsurface travel time ( STT ) in complex hillslopes

The travel time of subsurface flow in complex hillslopes (hillslopes with different plan shape and profile curvature) is an important parameter in predicting the subsurface flow in catchments. This time depends on the hillslopes geometry (plan shape and profile curvature), soil properties and climate conditions. The saturation capacity of hillslopes affect the travel time of subsurface flow. The saturation capacity, and subsurface travel time of compound hillslopes depend on parameters such as soil depth, porosity, soil hydraulic conductivity, plan shape (convergent, parallel or divergent), hillslope length, profile curvature(concave, straight or convex) and recharge rate to the groundwater table. An equation for calculating subsurface travel time for all complex hillslopes was presented. This equation is a function of the saturation zone length (SZL) on the surface. Saturation zone length of the complex hillslopes was calculated numerically by using the hillslope-storage kinematic wave equation for subsurface flow, so an analytical equation was presented for calculating the saturation zone length of the straight hillslopes and all plan shapes geometries. Based on our results, the convergent hillslopes become saturated very soon and they showed longer SZL with shorter travel time compared to the parallel and divergent ones. The subsurface average flow rate in convergent hillslopes is much less than the divergent ones in the steady state conditions. Concerning to subsurface travel time , convex hillslopes have more travel time in comparison to straight and concave hillslopes. The convex hillslopes exhibit more average flow rate than concave hillslopes and their saturation capacity is very low. Finally, the effects of recharge rate variations, average bedrock slope and soil depth on saturation zone extension were investigated.


Introduction
Subsurface flow is percolating water that encounters an impending horizon in shallow soil, where the water is diverted horizontally and reaches the stream channel.Due to the high permeability of topsoil and generally greater potential gradients in these upper sloping horizons, water following a topsoil path reaches the stream channel much quicker than the groundwater flow does.Some of this water arrives at the channel soon enough to contribute to the storm hydrograph and is classified as subsurface storm flow.The dynamic interaction between the saturated-unsaturated subsurface flow and surface flow has been examined by many researchers (Freeze and Harlan, 1969;Freeze, 1971Freeze, , 1972aFreeze, , 1972b;;Beven, 1982) through numerical simulations.The dynamics of water in a catchment and particularly at the surface/subsurface interface is still poorly understood.For simulating surface and subsurface flow in catchments, the changing of the saturated and unsaturated area by spatial and temporal rainfall distributions during storms is very important.
In the past, the Geomorphological Instantaneous Unit Hydrograph (GIUH) was used for simulating surface runoff (Rodriguez-Iturbe, 1979;Gupta et al., 1980;Rodriguez-Iturbe et al., 1982;Chutha and Dooge, 1990;Lee andYen, 1997, 2005;Olivera and Maidment, 1999).Recently, the GIUH model has been applied to consider both the surface and subsurface flow processes (Lee and Chang, 2005).This method is based on travel time probability distributions for runoff in surface flow and subsurface flow regions and channels.Travel time is defined as the average time required for water particles to travel from the top of the hillslope via the subsurface hillslope layers to the outlet.Henderson and Wooding (1964) simulated the surface and subsurface flow by using kinematic-wave approximation.The Henderson and Wooding equations showed that the travel time of the subsurface flow is proportional to the soil porosity and inversely proportional to the slope and hydraulic conductivity.Yet, their method cannot describe the effects of recharge rate, plan shape of hillslope(divergent, parallel, convergent), profile curvature (convex, planer and concave) and soil depth on subsurface travel time.The effect of the mentioned parameters on the surface travel time has been proved in past researches (Henderson, 1966;Eagleson, 1970;Overton and Meadows, 1976;Singh, 1882;Agiralioglu, 1985 andAkan, 1993).
In a simple hillslope experiencing a uniform net recharge, the analytical derivation of the response time behavior involves solving the one or two dimensional transient flow partial differential equation for hillslopes, popularly known as the Boussinesq equation (Boussinesq, 1877).
A general analytical solution to this non-linear equation has never been achieved.A number of researchers have solved simplified forms of this equation analytically, mostly for steady-state and for various special cases.Verhoest and Troch (2000), Troch et al. (2002) and Troch et al. (2004) developed analytical solutions for the Boussinesq equation using linearization and the method of characteristics, respectively.Huyck et al. (2005) developed an analytical solution to the linearized Boussinesq equation for realistic aquifer shapes and temporally variable recharge rates.Troch et al. (2003) and Hilberts et al. (2004) demonstrated that (numerical) solutions of the 1D hillslope-storage Boussinesq (hsB) equation account explicitly for plan shape (by means of the hillslope width function) and profile curvature (local bedrock slope angle and hillslope soil depth function) of the hillslope.To investigate the key role of geometric characteristics of hillslopes (plan shape and profile curvature) on shallow landslides, Talebi et al. (2008a) presented a steadystate analytical hillslope stability model based on kinematic wave subsurface storage dynamics.Comparison between the hillslope-storage Boussinesq and Richards' equation models for various scenarios and hillslope configurations shows that the hsB model is able to capture the general features of the storage and outflow responses of complex hillslopes (Paniconi et al., 2003;Hilberts et al., 2004).Berne et al. (2005) used the hsB model for the similarity analysis of subsurface flow response of hillslopes with complex geometry.He linearized the hsB equations by exponential width functions and introduced the hillslope Pe'clet number, an efficient similarity parameter for describing the hillslope subsurface flow response.Aryal (2005) and O'Loughlin (2005) have shown that the hillslope travel time in subsurface flow is dependent on hillslope length, hydraulic conductivity, plan shape, profile curvature and recharge rate.They demonstrated equations of saturation zone boundary for hillslope in steady state and introduced three equations for calculating complex hillslopes travel time based on Zaslavsky and Rogowski (1969) geometry equations.The objectives of this paper are: (i) introduce an equation for subsurface travel time of all complex hillslopes with regard to parameters such as the saturation zone length , total length, soil porosity, profile curvature, soil hydraulic conductivity , and average bedrock slope, (ii) calculate the saturation zone length of nine basic hillslopes in steady-state conditions, (iii) explore the effects of different factors such as the soil depth, the recharge rate, bedrock slope angle on travel time and saturation zone length, (iv) present analytical expressions for calculating saturation zone length in straight hillslopes for different shape functions (convergent, parallel, divergent) and finally, (v) compare the drainage capacity of all complex hillslopes based on their average discharge rates.
2 Model formulation 2.1 Hillslope geometry Evans (1980) characterized hillslopes by the combined curvature in the gradient direction (profile curvature) and the direction perpendicular to the gradient (contour or plan curvature).The surface of an individual hillslope is represented by the following function: where z is the elevation (m), x is horizontal distance measured in the downstream length (m) direction of the surface, y is the horizontal distance (m) from the slope centre in the direction perpendicular to the length direction (the width direction), E is the minimum elevation (m) of the surface above an arbitrary datum, H is the maximum elevation (m) difference defined by the surface, L is the total horizontal length of hillslope (m), n is a profile curvature parameter, and ω is a plan curvature parameter.Figure 1 shows a hillslope with a three-dimensional view of a convergent hillslope on top of an impermeable layer and a straight bedrock.Figure 2 illustrates nine basic hillslope types that are formed by combining three plan and three profile curvatures.The geometrical parameters for the nine characterized  hillslopes are listed in Table 1.The values of the hydrological parameters have been listed in Table 2.The assumptions applied to modeling subsurface government equations are: The saturated hydraulic conductivity is assumed to be uniform with depth, the hydraulic gradient is equal to the local surface slope, soil depth is uniform and recharge rate is constant (the steady state conditions).

The hillslope-storage kinematic wave equation
The soil moisture storage capacity (S c ) has been defined by Fan and Bras (1998) (Troch et al., 2003;Talebi et al., 2008a)  as: where f is the drainable porosity, w(x) is the width of the hillslope (m) at a distance x and D(x) is the average soil depth (m) at x (see Fig. 1a).S c defines the pore space along the hillslope and accounts for both plan shape, through the width function, and the profile curvature, through the soil depth function.
Similarly, the soil moisture storage S(x,t) has been defined by Troch et al. (2002) as: where h(x,t) is the average height over the width of the groundwater table at x and t.Introducing the integrated discharge over the width of the hillslope Q(x,t), the continuity equation becomes (Troch et al.2002): where N (t) is the recharge to the saturated layer (m/s).The subsurface flow rates can be described with a kinematic wave approximation of Darcy's law as (Troch et al., 2002): where z is the elevation of the bedrock above a given datum, k is the soil hydraulic conductivity (m/s).In the context of subsurface flow, it is reasonable to assume the following initial and boundary conditions: S(0,t) = 0 ∀t where g(x) represents the initial soil moisture storage along the hillslope.Troch et al. (2002) solved Eq. (4) analytically using the method of characteristics.The solution is given by: where A(x) is the upstream drainage area at location x (integral from 0 to x of w(x)).This equation expresses the storage profile along the hillslope in the steady-state condition.Analytical solutions to Boussinesq's equation are very useful to understand the dynamics of subsurface flow processes along a hillslope.

Prediction of the saturation zone in complex hillslopes
In this study, we used the steady state analytical solution of hillslope-storage kinematic wave equation that was presented by Troch et al. (2002) for predicting and extending the saturation zone in compound hillslopes.
Figure 3 shows a convergent hillslopes under recharge conditions.As can be seen, many parameters like recharge rate (N ), soil depth (D), hillslopes length (L), soil hydraulic conductivity (k), average slope (S) , profile curvature parameter (n), and plan curvature parameter (ω) affect the hillslope saturation zone extension.
According to Fig. 3 any point of the hillslope which the storage equals the storage capacity (S(x) = S c (x)), belongs to the saturation zone.If we call the ratio of actual storage to storage capacity as "Relative Saturation (σ )", one can say that any point of the hillslope where the relative saturation reaches one (σ ≥ 1), would be a saturation point.
The steady-state relative saturation function is now given by Talebi et al. (2008a): where T=kD is soil transmissivity (m 2 /s) and a(x)=A(x)/w(x) is drainage area per unit hillslope width (m).The variable σ (x)describes the steady-state wetness of the soil and is conceptually similar to the topographic index ln( a tanβ ) of Beven and Kirkby (1979), wetness index (W ) derived by O'Loughlin (1986) and Montgomery and Dietrich (1994).The location of the saturation zone boundary can be determined by inserting σ (x) = 1inEq.(8)orS(x) = S c (x) Relative saturated storage along the nine basic hillslopes for different recharge rates (solid line: N=30 mm/day; dashed line: N=20 mm/day; dotted line: N=10 mm/day).(k=0.0001m/s,D=2m,f=0.34,β=15 • ) and using the storage function from Eq. ( 7) and storage capacity from Eq. ( 2), we then obtain: By solving the Eq. ( 9) numerically, the location of the saturation zone boundary (x sat ) could be determined.The xcoordinate x sat , is where the mean groundwater table height is maximum.The saturation capacity beyond the saturation zone boundary (x sat ≤ x ≤ L) depends on the relative saturation at those points.

Calculation of the saturation zone length in the nine basic complex hillslopes
In general, relative saturation in the hillslope is a determiner of the soil saturation capacity.This parameter was also used by Troch et al. (2002) and Talebi et al. (2008a) in their research.Figure 4 shows variations of the relative saturation along the nine basic hillslopes of Tables 1 and 2 for different recharge rates.According to Fig. 4 all hillslopes react to the recharge variations differently.The saturation zone occurs in a certain recharge rate corresponding to the geometric attributes and the soil characteristics of the hillslopes which is called "Saturation Recharge Rate (SRR)".The recharge rate that causes the occurrence of the saturation zone in every hillslopes was calculated for all slopes.Figure 5 shows the SRR for nine basic of hillslopes.
The concave and convergent hillslopes are saturated very soon.As can be seen, the SSR in divergent hillslopes is averagely seven times more than the convergent slopes and the SSR in convex slopes is averagely nine times more than the concave slopes.For example, the SRR for convergent-concave hillslope is 2 mm/day (minimum rate) and it is 129 mm/day for divergent-convex hillslopes (maximum rate).
According to the studies on concave and straight hillslopes, the saturation zone in these hillslopes after saturation, occurs at the lower reaches of the hillslope between the edge of the saturated boundary and the hillslope outlet (ridge) completely, so the saturation zone length is obtained from the relation: SZL=L − x sat .In the case of a convex hillslope, the saturation zone occurs in the distance between the edge of saturation boundary and the hillslopes outlet; and close to the ridge, the storage is less than the storage capacity, as seen in Fig. 4. The relative saturated storage profile in concave and straight slopes is linear to parabolic and in convex slopes is semi-ellipse.
Figure 6 depicts the SZL of all complex hillslopes for various recharge rates (10 mm/day-30 mm/day).The recharge rate is a very effective factor in the saturation rate, for instance, the convergent-concave hillslopes with recharge under 20 mm/day show more reaction to the saturation rate in comparison to straight and convergent-convex hillslopes.In recharge rates over 20 mm/day the maximum of the SZL corresponds to the convex-convergent hillslopes.
In all recharge rates, the convergent hillslopes tend to saturate much more than the parallel and divergent ones.SZL in the convergent hillslopes are greater than parallel and divergent hillslopes.Greater SZL corresponds to concave hill- slopes compared to the convex ones.The Convex-divergent hillslopes show minimum response to saturation; therefore, in a recharge rate below 129 mm/day no saturation zone is created.In general, occurrence of the saturation zone in hillslopes causes an increase in pore pressure followed by a decrease in the stability of hillslopes.Talebi et al. (2008a) proved that the stability of the convergent hillslope is less than the divergent ones and the same is true about the concave hillslopes compared with the convex ones.
Soil depth is also an important factor affecting relative saturation.Figure 7 shows the change of relative saturation and SZL for the soil depths from 0.5 m-2 m.The less soil depth is created the more saturation zone, because the storage capacity is decreased and the soil will become saturated faster.
An increase in the soil hydraulic conductivity as well as the average bed rock slope yields a decrease in the saturation zone in complex hillslopes.Changes in the bedrock slope cause changes in the plan shape coefficient ω = ±(H /L 2 ), then changes the relative saturation of the hillslopes.Figure 8 depicts the effect of the average bedrock slope angle on the extension of the saturation zone .

Analytical solution for the saturation zone length in the straight hillslopes
In straight hillslopes there is no slope variation (n = 1) and the width function is as follows (Talebi et al., 2008a): of x becomes (Talebi et al., 2008a): Note that in convex/concave hillslopes, the drainage area in each point should be determined numerically and we cannot product an analytical equation for calculating the saturation zone in these slopes.The steady-state saturated storage profile for straight hillslopes based on (Talebi et al., 2008a) can be calculated as (Appendix A): Dividing Eq. ( 12) by Eq. ( 2), Talebi et al. (2008a) obtained the relative saturation function for straight hillslopes: belongs to the saturation zone.σ (x) ≥ 1 each point of the hillslopes with equating Eq. ( 13) to one, we derive the saturation zone length as follows: Where S is the average slope (=H /L).Equation ( 14) expresses the SZL in straight hillslopes.The SZL depends on the recharge rate, the plan shape, the soil hydraulic conductivity, the soil depth and the total hillslope length.If solution of the Eq. ( 14) is negative or a complex number (no valid value), it emphasizes that the saturation zone does not exist, so SZL is zero; otherwise, if the solution is positive.The plan shape parameter(ω) of straight-parallel hillslopes is zero, by limiting ω toward zero; the Eq. ( 14) is changed to: Positive solutions of the Eq. ( 15) present the saturation zone length in straight-parallel hillslopes after saturation.

The complex hillslope travel time in the subsurface flow
The time of concentration has been used by some authors to define response times of hillslopes.Ben-Zvi (1984) defined the time of concentration as the time taken from the initiation of rainfall to the time when the catchment discharge attains (nearly) 0.8 of the equilibrium discharge.Beven (1982) defined the time of concentration as the time at which a steadystate flow profile is developed over the entire hillslope, assuming a constant input rate for a sufficient length of time.In this paper, the time of concentration is defined as the average time required for water particles to travel from the top of the hillslope, via the subsurface hillslope layers, to the outlet.Aryal et al. (2005) used the Zaslavsky and Rogowski (1969) geometry equations and derived three equations for travel time in hillslopes with concave , convex and straight profiles and all plan form geometries, however, in this paper, since we have used the Evans (1980) equation for modeling of slopes geometry, one equation is presented for all complex hillslopes.In a soil profile over an impermeable layer, the interstitial velocity in soil pores according to Darcy's law is (Aryal2005): where s * is the local slope.The profile curvature affects slope changes and the velocity of water in soil.The local slope for the compound hillslopes is derived from Eq. ( 1): Putting v = dx/dt and substituting the value of s * [from Eq. ( 17) ] in Eq. ( 16) and integrating with bounds t = 0 to T , and x = 0 to x sat for travel time of the unsaturation zone, we obtain: Note that when the saturation occurs over part of the lower hillslope, the total travel time diminishes.In this case, the overall travel time for subsurface flow is reduced by the Hydrol.Earth Syst.Sci., 14, 891-900, 2010 www.hydrol-earth-syst-sci.net/14/891/2010/ travel time required to traverse the saturation zone length (see Marani et al., 2001;Aryal et al., 2002).The analytical Eq. ( 18) represents the subsurface travel time in all complex hillslopes.The coordination of the saturated zone boundary in each hillslope is a key parameter in calculating travel time.All parameters which affect the development of the saturation zone of the hillslopes also change the travel time.Equation (18) expresses that the travel time of the nine hillslopes is a function of the saturation zone length (SZL = L − x sat ), the total length (L), the effective porosity (f ), the profile curvature parameter (n), the soil hydraulic conductivity (k) , and the average bedrock slope (S).
The subsurface travel time in the steady-state conditions involves the storage rate in the system, and the outlet discharge.Inserting Eq. ( 17) into Eq.( 5) gives the ratio of the storage to the outlet discharge in this case: This ratio for hillslopes with constant profile curvature remains the same, but varies along the hillslope.Combining Eq. ( 18) and Eq. ( 19) one can write (Talebi et al., 2008c): The value of the storage in the system is estimated from Eq. ( 7) and the outlet discharge from Eq. ( 5).By using Eq. ( 20), we derive the travel time from the area under the graph of S/Q along the unsaturation zone length.Both Eqs. ( 18) and (20) describe the subsurface travel time in complex hillslopes but Eq. ( 18) is simpler and avoids any calculations for S(x) and Q(x).In slopes with fixed profile curvature, the ratio of the storage to the flow rate in steady-state conditions remains constant; hence this is the unsaturation zone length (the effective length) which influences the travel time in these hillslopes.

The subsurface travel time in straight hillslopes
In straight hillslopes (n=1) Eq. ( 18) becomes: This equation presents the subsurface travel time of the straight hillslopes after saturation.Replacing Eq. ( 14) into Eq.( 21) yields (Appendix B): As seen in Eq. ( 22), the subsurface travel time of the straight hillslopes after saturation is a function depending on the recharge rate, the plan shape, the soil hydraulic conductivity, the soil depth, the hillslope length, the soil porosity and the bedrock slope angle.Equation ( 22) does not relate to SZL and profile curvature.This equation is an extended equation of the Henderson (1964) equation.The runoff travel time for the subsurface flow region that was presented by Henderson and Wooding (1964) is: Most researchers have used Eq. ( 23) for calculating subsurface travel time in overlands in order to predicting the subsurface flow hydrograph.In Eq. ( 23) the effect of recharge rate, geometry and soil depth has been ignored.This equation is only valid for the straight hillslopes before the saturation conditions.The Eq. ( 22) is simplified for straight-parallel (ω = 0) hillslopes to: It can be stated that the travel time in straight-parallel hillslopes is a function of f L/kS before saturation and is a function of Df/N after saturation develops.This concept has also been proved by Aryal et al. (2005).They proved that the travel time in straight-parallel hillslopes is a function only of smd/∇q after saturation occurs, where ∇q is the net change in flux and smd is the soil moisture deficit.The relationship between the initial soil moisture and the soil moisture deficit is:

Calculation of the subsurface travel time in the nine basic hillslopes
The travel time and SZL of the all complex hillslopes according to the attributes in Table 1 and Table 2 are presented in Table 3. Figure 9 represents also the variations of the SZL as the recharge reaches 50 mm/day.Figure 10 shows the subsurface travel time of the nine basic complex hillslope for N=50 mm/day.The histograms showed in Fig. 10 and Table 3 illustrate that the convergent hillslopes exhibit less travel time than parallel and divergent ones and it is also more for the convex hillslopes compared to the straight and concave hillslopes.The travel time in divergent hillslopes is approximately double those of convergent hillslopes.
As verified by our studies, the least travel time corresponds to the concave-convergent hillslopes and the greatest to the convex-divergent ones.When the saturation zone length is increased, the length of unsaturation zone is decreased and the effect of the surface flow is more important than the subsurface flow.In this situation the travel time of the subsurface flow will be reduced.
The steady-state outflow at each point of the hillslope is equal to: The average outflow along hillslopes would be: The average subsurface flow can be obtained from Eq. ( 27) for nine basic hillslopes (see Fig. 11).Figure 11 illustrates that the average flow of the convergent hillslopes is less than the divergent ones.Also the convex hillslopes tend to show much more flow than the concave ones.The least flow relates to the concave-convergent hillslope, with 0.04 m 3 /hr while the highest corresponds to the convex-parallel, with 3 m 3 /hr .2005) also measured the water tables and outflow rates from a drainage experiment in a laboratory setup by two sets of linear convergent hillslopes and linear divergent hillslopes and our results are consistence with their results as they showed that the convergent hillslopes drain more slowly than the divergent ones.

Conclusions
In this paper, we proved the hillslope-storage kinematic wave model is suitable for investigating the response of the complex hillslopes and some of our results are similar to the Aryal et al. (2005) results in the steady state condition but the hsB model can be also extended for unsteady state condition.Troch et al. (2003) and Talebi et al. (2008b) have used the unsteady-state hsB model in their researches.The main aim of the present study is to benefit from its results for the modeling saturation zone extension in unsteady-state condition based on temporal distributions of rainfall during storms in hillslopes in future studies.The convergent hillslopes possess less flow rate discharge in comparison to the parallel and divergent ones, a fact admitted by the experimental results obtained by Troch et al. (2003) andHilbert et al. (2004) in the laboratory.In convergent hillslopes the ground water table is higher than the divergent ones, leading to more saturation with larger saturation zone length in contrast to parallel and divergent hillslopes.
Since the travel time is determined along the unsaturation zone, the saturation zone length reduces the subsurface travel time.In hillslopes with fixed profile curvature (convergent, parallel, divergent), the ratio of the storage to the flow rate in steady-state conditions remains constant, hence, this is the effective length which influences the travel time in these hillslopes.The maximum saturation zone length in convergent hillslopes, it is inferred that they have the minimum effective length, with shorter subsurface travel time relative to the divergent and parallel hillslopes.Any alteration in the plan shape makes change(s) in the effective length, eventually resulting in changes in the travel time.So the travel time is a function of the profile curvature, the plan shape, and the characteristics of the soil and recharge.
Based on our results, the convex hillslopes show smaller saturation zone than the concave hillslopes and take greater travel time than the straight and concave ones; whereas in convex-convergent hillslopes due to their convergent, there is much inclination towards saturation.
The least travel time corresponds to the concaveconvergent hillslopes, and the greatest to the convexdivergent ones.

Fig. 1 .
Fig. 1.(a) A three dimensional view of a convergent hillslope overlying a straight bedrock profile, (b) a definition sketch of the cross section of a one-dimensional hillslope aquifer overlying a bedrock with a constant bedrock slope angle (after Talebi et al., 2008a).

Fig. 2 .
Fig. 2. A three-dimensional view (top) and a two-dimensional plot of the contour lines and slope divides (bottom) of the nine hillslopes considered in this study (after Hilbert et al., 2004).

Table 3 .
SZL and STT in the nine basic hillslopes (N=30 mm/day)