Natural Hazards and Earth System Sciences On the application of kinematic models to simulate the diffusive processes of debris flows

Debris flows generally propagate along steep mountain torrents with dynamics primarily governed by gravitational and frictional forces. Thus, debris flows modelling can be successfully performed through the application of kinematic models, which consider only the effects of slope and friction and neglect the remaining terms of the momentum equation. However, the diffusion processes that can be observed in the field, such as the spreading of the debris flow wave as it flows downstream, can not be theoretically predicted by kinematic models, since diffusion is a second-order process neglected in the kinematic approximation. In this paper, this issue is discussed and an application for both a generalized diffusion wave model and a kinematic model is proposed of a debris flow which occurred in an Italian instrumented torrent to identify, in a real case scenario, the effective value of the neglected terms in the kinematic approximation.


Introduction
Debris flows usually occur in mountain torrents characterised by steep channel slopes.Debris flows may significantly differ as far as the number of surges following the main wave (or superimposed to it), maximum discharge, front height and magnitude are concerned.These general characteristics may vary from case to case, depending on triggering conditions, channel geometry and the rheology of the mixture (Takahashi, 1991).Nevertheless, debris flows share some similarities.In particular, debris flows show two diffusive characteristics as they travel downstream: 1. spreading of the hydrograph; Correspondence to: M. Arattano (massimo.arattano@irpi.to.cnr.it) 2. decay of front height (Pierson, 1986;Johnson, 1970;Genevois et al., 2000;Arattano and Marchi, 2008;Iverson, 1997;Hürlimann et a., 2003).
Similar characteristics can also be observed in flume experiments performed to simulate debris flow mechanics (Pudasaini et al., 2005;Lanzoni and Tubino, 1993).The reduction of the front height and the spreading of the debris flow hydrograph may also be observed for lahars (Pierson et al., 1990).
Mathematical models are commonly employed to simulate the propagation of debris flows along the torrent and predict the front height evolution, the hydrograph deformation and then the deposition on the fan.The simulation of the propagation phase of the flow within the channel boundaries is normally performed through 1-D models, using the so called Saint-Venant equations, while the deposition is simulated through 2-D models.
A flow that takes place in a steep channel is governed primarily by gravitational and frictional forces and, to a much smaller extent, by the local and convective inertial forces, which originate from the local variation of velocity and by the flow depth gradient, respectively.For this reason, kinematic models have been proposed in literature to simulate the propagation on steep slopes of either water waves (Lighthill and Whitham, 1955) and debris flows (Arattano and Savage, 1994).Kinematic models, in fact, neglect the local variation of velocity and the flow depth gradient in the Saint-Venant momentum conservation equation and take into account only the gravitational and frictional forces.Kinematic models apply when the flow takes place along steep channels, as it actually occurs for the mountain torrents along which debris flows propagate.However, from a theoretical viewpoint, kinematic waves cannot diffuse because the diffusion is given by those very two terms mentioned above that are neglected in the equation of motion (the momentum equation).Kinematic waves can convect downstream and transport mass in the process, but they cannot dissipate their discharge or flow stage.Therefore, the effective capability of kinematic models to simulate real debris flows and their observable diffusive characteristics as they travel downstream, as in the case of the Arattano and Savage (1994) model, should be better understood and explained.Moreover, the effect of the diffusive terms should be quantified to verify its amount and a discussion is needed about the possibility to adequately simulate debris-flow processes through kinematic models.

The kinematic model
Modelling a water-sediment flow on a steep slope is generally made using the Saint-Venant equations, which are obtained by assuming a hydrostatic pressure distribution over the flow depth, no erosions and depositions and slowly varying cross-sections.The Saint-Venant equations consist of the continuity equation, which is (assuming a rectangular crosssection shape): and the momentum equation: In Eq. ( 2), the terms on the left side are the bed slope and the pressure gradient and the terms on the right side are the local and convective inertia and the energy gradient, respectively.A kinematic wave can be defined as a shallow water wave with the assumption that the gravitational and frictional forces prevail on the inertial forces, so that the momentum equation reduces to S 0 =S.Consequently, in kinematic waves a direct relationship between the flow velocity and the flow depth holds: The coefficients k s ,n, in Eq. (3) depend on the rheological behaviour of the flowing mixture and may vary in a wide range (Arattano et al., 2006).According to Ponce (1992), a kinematic wave can be defined in different ways, first and foremost as a wave that transports mass, as expressed in Eq. (1).As pointed out by Ponce (1991Ponce ( , 1992)), Eq. (3) plays a key role in the model.In fact, neglecting to take into account the terms for the local variation of velocity and the flow depth gradient, it imposes the physical and mathematical constraint that the wave cannot diffuse.In fact, diffusion is a second-order process that is described precisely by those neglected terms (local inertia, convective inertia and pressure gradient).Therefore, a kinematic wave cannot show any dissipation, in particular, it should not be able to account for any spreading in space and time of both the discharge and flow stage (Ponce, 1991(Ponce, , 1992)).
If some sort of dissipation was observed in a kinematic wave numerical simulation, it would only be due to the conversion of the partial differential Eq. (1) into a finite difference equation to run the simulation (Ponce, 1991(Ponce, , 1992)).The amount of this numerical dissipation would depend on the chosen space and time grids ( x, t) and would propagate in the numerical computation, thus, affecting the final results.Consequently, the comparison between recorded data and the results obtained by the mathematical simulation can be misleading, as numerically induced errors could not be negligible.
While all this applies to numerical modelling, analytical versions of the kinematic model should not show any dissipation at all, because in the analytical models no numerical integration schemes are used (e.g.finite difference or finite elements schemes) and, therefore, no numerical diffusion is possible.Nevertheless, in 1994 Arattano and Savage proposed an analytically deduced kinematic wave model for debris flows that is apparently (and paradoxically) capable of simulating diffusive characteristics.It consists of two parametric equations, obtained by assuming a triangular shape of the wave at t=0 as one of the boundary conditions.
The two parametric equations of the Arattano and Savage (1994) kinematic model give the position of the front: and the flow depth behind the front: The model, which has been successfully applied to field data (Arattano and Savage, 1994;Arattano et al., 2006), analytically predicts a progressive reduction of the front depth as the wave flows downstream and a spreading of the debris flow hydrograph that corresponds to a continuous increase of the wave length.These are two typical diffusive manifestations that a kinematic model should not be able to simulate.Thus, this apparent paradox needs a brief discussion.
One of the assumptions made by the authors is that the total volume of the debris flow remains constant during the entire propagation process, that is no erosions or depositions take place or erosions and depositions balance out (Arattano and Savage, 1994).At time t=0, the shape of the debris flow is assumed to be triangular, with a total length L and front depth H (Fig. 1), according to: The assumption of a constant volume is expressed as follows Nat.Hazards Earth Syst.Sci., 10, 1689-1695, 2010 www.nat-hazards-earth-syst-sci.net/10/1689/2010/(assuming a rectangular shape of the channel): Integration of Eq. ( 8) starts at x=0 and stops at x=x f where the integral is equal to the total volume.Since the starting point of the integration remains fixed (x=0), the integration is implicitly carried out assuming that, for any given ε>0 arbitrarily chosen: As a consequence of the assumption of a constant volume and of the integration over the interval x=[0;x f ], the front height must progressively decrease.This is a result of the particular choice of the boundary conditions needed to solve the system equations given by Eq. ( 1) and Eq. ( 3).This can be easily understood if we consider that the portion of the wave behind the peak, that is the decreasing limb of the hydrograph, follows the peak with a slower velocity.This is due to the direct dependency of velocity on flow height, as expressed by Eq. ( 3), and to the progressive decrease of flow height behind the main wave front due to the assumption of a triangular form for the initial debris mass.Thus, the whole debris flow wave elongates and the continuity equation then imposes the peak height to decrease.It must be noticed that a different choice of the boundary conditions would have lead to completely different results.As a matter of fact, very different boundary conditions can be imposed for a kinematic model and these can then greatly influence the final results.As an extreme example, it could even be possible to impose that the integration starts from the front that flows downstream with a celerity given by: that ranges between the front itself x f and the position of the tail.In this way, the integral (8) becomes: where x(h=H ) is the abscissa where the debris flow height is h=H , and x t is the tail of the debris flow, where the flow depth is not necessarily equal to zero.
Even though this would lead to unrealistic results, it could nevertheless be mathematically imposed and the equation could be consequently solved to obtain: These equations predict that the tail thickens downstream and that the head of the debris flow remains unchanged (Fig. 2).Note that the latter is a consequence of the particular choice of the integration interval and it is not imposed a-priori.
As mentioned earlier, these results are unrealistic but they are presented here to show the strong dependency of the solution on the specific boundary conditions that are chosen.The capability of the Arattano and Savage (1994) model to simulate debris flow propagation and the diffusive effects that take place in nature only derives from a peculiar choice of the boundary conditions, while the model still remains incapable of simulating the natural diffusions taking place along the debris flow wave.
Since we have shown that there is no real capability to simulate any real diffusion through a kinematic wave, it remains to be established if the amount of effective diffusion actually taking place in debris flow waves is small enough to allow the use of such an approximation to satisfactorily simulate debris flow propagation.In the following, a nonlinear model proposed by Todini (2007) will be applied to investigate this issue.The model will be used to simulate a real debris flow that took place in Rio Moscardo (North-Western Italy, Fig. 3) on 23 July 2004, and to compare the respective roles played by the diffusion and by the kinematic terms.A detailed description of the debris flow that has been simulated can be found in Arattano et al. (2006).In the following figure, a sketch of the instrumented watershed basin is shown, with the location of the ultrasonic gauges.
The debris flow stage was measured by two ultrasonic gauges, at two different cross-sections.The lymnographs allowed us to investigate on the rheology of the recorded debris flow (Arattano at al., 2006) by means of a mathematical model.In the latter, the upstream hydrograph (Fig. 4) was used to set the upstream boundary conditions.
The simulated results have been compared to the flow heights recorded at the downstream cross-section by the second ultrasonic sensor.

Kinematic-diffusive model
As previously mentioned, Eq. (3), which is a simplified form of the Saint-Venant momentum equation, implies that no physical diffusion is taken into account in the simulations, since diffusion is a second-order effect.However, physical diffusion cannot be neglected a priori.In fact, along the increasing limb of the hydrograph in particular, the depth gra- dient is no more negligible and, as a consequence, Eq. ( 3) does not hold any more.The discharge is higher than that calculated with Eq. ( 3) along the rising limb of the hydrograph, and it is smaller along its falling limb, as it is expressed by the so-called looped stage-discharge curves (Chow, 1959).Lighthill and Witham (1955), in order to take into account a small amount of diffusion, tried to extend the kinematic model by considering the terms that depend on the flow gradient and neglecting the momentum source terms and the terms that depend on local inertia.
By considering the continuity equation and the momentum equation together (where the local inertial and convective inertial terms are neglected), it is obtained (Ponce, 1991(Ponce, , 1992;;Todini, 2007): with c and υ that are the kinematic wave celerity and hydraulic diffusivity, respectively, defined as: Equation ( 14) can be solved either analytically or numerically.
In the present work, a nonlinear model recently proposed by Todini (2007), for the computation of Q (c and υ vary with time) will be used to compare the roles played by the diffusion terms and by the kinematic terms.
The computation of h, can be made by solving the continuity equation, Eq. ( 1), by means of a first order approximation: where Q n i and h n i represent the numerical approximation: Nat. Hazards Earth Syst.Sci., 10, 1689-1695, 2010 www.nat-hazards-earth-syst-sci.net/10/1689/2010/Q(x i ,t n );h(x i ,t n ), with x i = x 0 + i x; t n = t 0 + n t, being i = 1,2,...,M; n = 1,2,...,N, and α, β numerical parameters of integration (Abbot and Cunge, 1982).The model gives Q = Q(c,υ), where both celerity c and diffusivity υ depend on the Courant number Cn, defined as follows (Todini, 2007): wave.celerity computation.celerity(18) The model allows the computation of the total height of the debris flow in space and time, it is mass conservative, and allows the comparison of the role played by the diffusion terms in Eq. ( 14) and that played by kinematic terms.
4 Application of the kinematic-diffusive model Todini's (2007) diffusive model was applied to simulate the propagation of a debris flow which occurred on 23 July 2004 in the Moscardo torrent, an instrumented basin located on the northeastern Italian Alps (Fig. 3) (Arattano et al., 1997).
The debris flow depths were measured at two different cross-sections, located 75 m apart, by means of two ultrasonic sensors.The availability of two hydrographs (Figs. 4  and 5) has allowed the calibration of the rheological parameters of Eq. ( 3) (Arattano et al., 2006).The upstream hydrograph has been used as an upstream boundary condition, while the downstream hydrograph has been used to compare the simulation results and the recorded data.
The values of the rheological parameters of Eq. ( 3) have been assumed to be the same as those obtained by Arattano et al. (2006), employing the complete Saint-Venant equations, namely k s =4 m −0.8 /s and n=1.2.The results (Fig. 5) show a good match between computed and recorded hydrographs.
Note that the passage of the debris flow has increased the bed elevation at the downstream station.This could be due to the filling of some erosion under the downstream ultrasonic sensor that took place before the arrival of the debris flow.Thus, the bed elevation surveyed after the passage of the debris flow has been used in the simulation and this explains the significantly different flow stage that the simulation and the recordings show in Fig. 5 before the arrival of the debris flow front.
The goal of evaluating the amount of physical diffusion in the recorded debris flow has been pursued calculating the ratio: This ratio compares the sum of the first and second terms of Eq. ( 14) to the third one.It is noticeable that if in Eq. ( 14) one assumes υ=0, the model expressed only by the sum of the first and second term, that is:  is the same as the kinematic model.Consequently the higher the rate given by Eq. ( 19), the higher the influence of kinematic terms is, with respect to the diffusive term (second partial derivative).In Fig. 6, a plot with the value of this rate for the downstream hydrograph is shown.Figure 6 shows that the diffusion term can be of the same order of magnitude as the remaining terms.This means that diffusive effects cannot be neglected a priori.Their main effects are exerted in the front of the debris flow, where the curvature of Q is also not negligible.Therefore, the application prefers a mathematical model capable of taking into account the diffusion terms.
However, the uncertainties regarding the determination of the rheological parameters, together with their influence on the simulation results and with the influence of channel geometry, can be so great to significantly hide, in some cases, the diffusion effects.This latter is probably the main reason why the Arattano and Savage (1994) model has been capable of such a good representation of the field observations, even though it cannot represent any actual diffusive effect.Since Todini's (2007) model is applied through a numerical simulation, it has inevitably introduced some amount of numerical diffusion.In order to compare the effects due to numerical and physical diffusion, it is possible to use the cell Reynolds number, D, introduced by Todini (2007).
D is given by the ratio between the physical and numerical diffusivity, that is: The larger D is, the higher is the role of the physical diffusivity.
In the simulated debris flow (Fig. 7) D is always larger than 1, with a maximum of 5 in the proximity of the front of the debris flow.This confirms that, assuming x=1 m and t=1 s, the physical diffusivity prevails on the numerical one.
Thus, a numerical model can be safely employed without greatly affecting the results because of the numerical dissipation.

Conclusions
The diffusive processes that take place in debris flow can be mathematically simulated only if second order terms are taken into account in the Saint-Venant equation.Consequently kinematic models, though widely used in the literature, cannot fairly reproduce the dissipation that take place in real debris flows.However, due to a particular choice of boundary conditions, it has been shown that a kinematic model expressed in analytical form allows the simulating of two large-scale diffusion effects: the downstream spreading of the debris-flow wave and the progressive decrease of its main front height.Then a mathematical convective-diffusive model has been used to simulate a field debris flow which occurred in 2004 in an Italian torrent, whose rheological behaviour has been described by means of a Chézy-like for-mula.The effects of local diffusion have been determined and compared with the effects due to other terms.Physical diffusive processes and processes due to numerical diffusion have also been quantified and compared.The results showed that physical diffusion processes are not negligible with respect to kinematic (mass transfer) processes.Therefore, the application of a mathematical model capable of taking into account the diffusion terms should be preferable and it has been shown that a numerical model can be safely employed without greatly affecting the results because of the numerical dissipation that it inevitably introduces.

Fig. 5 .
Fig. 5. Comparison of recorded and simulated hydrograph (simulation performed with x=1 m; t=1 s) at the downstream reach.
the flow c = kinematic celerity Cn = Courant number (dimensionless) D = cell Reynolds number (dimensionless) h = flow depth H = initial height of the debris flow front h f = height of the front of the debris flow h tail = height of the tail of the debris flow k s ,n = rheological parameters L = initial length of the debris flow M = total number of space steps of integration N = total number of time steps of integration Q = total discharge (water and sediment) S = energy gradient (dimensionless) S 0 = channel bottom slope (dimensionless) t 0 = time at which h f is at the x f progressive U = mean velocity of the flow x, t = independent variables (time and space) x f = abscissa of the front of the debris flow x t = abscissa of the tail of the debris flow x, t = finite differences in x and t α, β = numerical parameters (dimensionless) υ = hydraulic diffusivity