Table Of Content

Shock-Capturing and Front-Tracking Methods for Granular Avalanches 5 Y.C. Tai1, S. Noelle2, J.M.N.T. Gray3 & K. Hutter1 1 1Institutfu¨rMechanik,TechnischeUniversitätDarmstadt,64289Darmstadt,Germany. 0 [email protected],[email protected] 2 2Institutfu¨rAngewandteMathematik,UniversitätBonn,53115Bonn,Germany. n [email protected] a 3DepartmentofMathematics,UniversityofManchester,ManchesterM139PL,UK. J [email protected] 0 2 Shockformationsareobservedingranularavalancheswhensupercritical ] A flow merges into a region of subcritical flow. In this paper we employ a shock-capturing numerical scheme for the one-dimensional Savage-Hutter N theoryofgranularflowtodescribethisphenomenon. ALagrangianmoving . h mesh scheme applied to the non-conservative form of the equations repro- t a duces smooth solutions of these free boundary problems very well, but m fails when shocks are formed. A non-oscillatory central difference (NOC) [ schemewithTVDlimiterorWENOcellreconstructionfortheconservative 1 equations is therefore introduced. For the avalanche free boundary prob- v lemsitmustbecombinedwithafront-trackingmethod,developedhere,to 6 properly describe the margin evolution. It is found that this NOCscheme 5 combined with the front-tracking module reproduces both the shock wave 7 4 and the smooth solution accurately. A piecewise quadratic WENO recon- 0 structionimprovesthesmoothnessofthesolutionnearlocal extrema. The . 1 schemesarecheckedagainstexactsolutionsfor(1)anupwardmovingshock 0 wave,(2)themotionofaparaboliccapdownaninclinedplaneand(3)the 5 motion of a parabolic cap down a curved slope ending in a flat run-out 1 : region, where a shock is formed as theavalanche comes to a halt. v i X r Key Words: Granular avalanche, shock-capturing, non-oscillatory central scheme, free a moving boundary,front-tracking. AMS Subject Classification code: 65M06, 35L67, 35R35, 86A60. 1. INTRODUCTION Snow avalanches,landslides, rock falls and debris flows are extremely dangerous anddestructivenaturalphenomenaofwhichtheoccurancehasincreasedduringthe past few decades. Their human impact has become so significant that the United Nations declared 1990–2000 International Decade for Natural Disasters Reduction (IDNDR). Research on the protection of habitants from floods, debris flows and 1 2 Y.C.TAI,S.NOELLE,J.M.N.T.GRAY&K.HUTTER avalanches is under way worldwide, and many institutions focus on the numerical prediction of such flows under ideal as well as realistic conditions. One of the models that has become popular in recent yearsis the Savage-Hutter (SH) avalanche theory for granular materials [33, 34]. In the past decade numerical techniques were developed to solve the SH-governing differential equations for typical moving boundary value problems [33, 34, 13, 9, 10, 19, 16, 14, 8, 41, 6, 7]. These techniques are based on a Lagrangian moving mesh finite-difference scheme in which the granular material is divided into quadrilateral cells (2D) or trian- gular prisms with flat tops (3D). Exact similarity solutions of the SH-equations were constructed in spatially one-dimensional chute flows [33, 35, 31] and for two- dimensional unconfined flows [15, 12]. In the case of chute flows it was shown that the solutions obtained by the Lagrangian integration procedure approximate the exactparabolicsimilaritysolutionveryaccurately,andthesetheoreticalandnumer- icalresultsareingoodagreementwithexperimentalavalanchedata. Similaragree- ment between theoretical, numerical and experimental data was also obtained for thetwo-dimensionalflowconfigurations(cf. abovereferences). IntheseLagrangian schemes explicit artificial numerical diffusion was incorporated to maintain stability. Indoingsothequalityofresolutiondeteriorates. Infact,theadequacyofthese numericalsolutionscanbechallengedbecauseofuncontrolledspreadingduetothis diffusion. ItwasalsoobservedthattheLagrangianschemesloosetheirstability(or else unjustified artificial diffusion must be applied) whenever internal shocks are formed. This appears to occur whenever the avalanche moves from an extending to a contracting flow configuration. These shocks are travelling waves which form bumps with steep gradients on the free surface, which is thicker on the down slope side. Itisthereforenaturaltodevelopconservativehigh-resolutionshock-capturing numerical techniques that are able to resolve the steep surface gradients and iden- tify the shocks often observed in experiments but not captured by the Lagrangian finite difference scheme. The development of high-resolution shock-capturing schemes has a long history which we cannoteven sketchhere (see e.g. the classicalreferences [3, 40, 11, 42] or the recenttextbooks [25,4, 20, 39]. The mostcommonapproachis to first develop a one-dimensional TVD (total-variation-diminuishing) upwind scheme for a scalar conservationlaw andthen apply it to systems using one-dimensionalcharacteristic decompositions or approximate Riemann solvers. Upwind schemes have been used very successfully for gas dynamical calculations, where the Riemann problem can be solved exactly and many approximate Riemann solvers are available. For more complicatedsystemslikethe granularflowmodelconsideredherecharacteristicde- compositions are often not available, and the Riemann problem cannot be solved analytically. Therefore we have chosen an alternative approach to high-resolution shock-capturing,namely the recentnon-oscillatorycentral(NOC) schemes first introduced by Nessyahu and Tadmor [30]. While upwind schemes are higher order extensions of the classical Godunov scheme, central schemes build upon the (also classical) Lax-Friedrichs scheme [23]. This scheme avoids characteristic decompositions and Riemann solvers by the use of a staggered grid. When used together with piecewise constant spatial reconstructions, the Lax-Friedrichs scheme is more diffusive than Godunov’s scheme. However, when one combines the scheme with TVD-type piecewise linear reconstructions, it becomes competitive with the up- SHOCK-CAPTURINGANDFRONT-TRACKINGFORGRANULARAVALANCHES 3 wind schemes. Recently, central schemes have been extended in many directions, seeeg. [1,2,18,27]formultidimensionalextensions,[32]for anadaptivestaggered scheme,[28,26]forthird-andhigher-orderschemesand[22,21]forcentralschemes onnon-staggeredgrids,whicharepreciselyatthe borderlineofcentralandupwind schemes. HereweadaptthesecondorderNOCschemeofNessyahuandTadmortoinclude anearthpressurecoefficientwhichhasajumpdiscontinuityastheflowtravelsfrom anexpandingintoacontractingregion,andtotreatthesourcetermwhichisdueto the spatiallyvarying topographyandthe gravitationalforce. The resultingscheme works well both in smooth regions and at shocks, which are captured within two mesh cells and without any oscillations. Besides the formation of shock fronts in the interior,avalanchesmay also have a vacuum front at their margins. Similarly as for the equations of gas dynamics, the hyperbolicsystemdegeneratesatthe vacuumstate. Manyshockcapturingupwind schemes produce negative heights at these points and subsequently break down or become completely unstable. While our NOC scheme is remarkably stable at the margins, it does not capture the vacuum front as well as the Lagrangian moving mesh scheme. To overcome this imperfection, we augment the NOC scheme with an algorithm that tracks the vacuum front. The combined front-tracking non- oscillatorycentralschemeisaccurateandrobustbothatshocksandatthemargins of the granular avalanche. The ensuing analysis commences in 2 with the presentation of the governing § SH-equations in conservativeand non-conservativeform; then the jump conditions ofmassandmomentumatsingularsurfaceswillbestatedandthesolutiontoasin- gleshockwave(ahydraulicjump)willbepresented; 2closeswiththeconstruction § of exact similarity solutions of a parabolic heap moving down a rough incline. 3 § introduces the numerical techniques; at first the Lagrangian integration technique is described; itis followedbythe presentationofthe non-oscillatorycentral(NOC) scheme. In 4 we augment the NOC scheme (which uses a fixed Eulerian grid) § with a Lagrangiantype front-tracking method in the marginal cells. 5 elaborates § on numerical results. The travelling shock wave cannot be handled by the La- grangianmethod, but the NOC scheme can do so with very little diffusion accross the shock. On the other hand, the parabolic similarity solution is well produced by the Lagrangian integration technique, but much less accurately by the NOC schemes unless the Lagrangian front-tracking is introduced for the marginal cells. It is alsoshownthat the NOC scheme with piecewise linearspatial reconstructions applying standard TVD type slope limiters exhibits some oscillations near smooth localmaxima. We removethese oscillationsbyincorporatingapiecewise quadratic WENO (weighted essentially non-oscillatory)reconstruction into our scheme. Our final numerical experiment combines all the difficulties treated in the paper: an avalanche with a vacuum front at the margins expands as it flows downhill and contractsasithits theflatrunout(sotheearthpressurecoefficientchangesdiscon- tinuouslyatthe transitionpoint). As the avalanchecomes toa haltatthe bottom, a shock wave develops and propagates upslope. Our NOC front tracking scheme handles this challenging flow very satisfactorily. 6 presents conclusions and gives § an outlook to further work. 4 Y.C.TAI,S.NOELLE,J.M.N.T.GRAY&K.HUTTER 2. GOVERNING EQUATIONS AdetailedderivationoftheSavage-Huttertheoryhasbeengivenin[33,34]. Here weconfineourselvestoabriefdescription. Althoughcohesionlessgranularmaterials exhibit dilatancy effects numerous experiments have confirmed that during rapid dense flow it is reasonable to assume that the avalanche is incompressible with constant uniform density ρ .” During flow the body behaves as a Mohr-Coulomb 0 plastic material at yield. As the avalanche slides over the rigid basal topography a Coulomb dry friction force resists the motion. The basal shear stress is therefore equal to the normal basal pressure multiplied by a coefficient of friction tanδ, where δ is termed the basal friction angle [19]. Scaling analysis isolates the physically significant terms in the governing equations and identifies those terms that can be neglected. Plane flow configurations are our focus in this paper, so depth integrationreducesthe theoryto onespatialdimension. The leadingorder,dimen- sionless, depth integrated equations for the local thickness of the avalanche h and the momentum hu (u is the downslope velocity) reduce to ∂h ∂ + (hu)=0, (1) ∂t ∂x ∂(hu) ∂ + (hu2+β h2/2)=hs (2) x x ∂t ∂x with net driving force ∂zb s =sinζ sgn(u)tanδ cosζ+λκu2 εcosζ , (3) x − − ∂x (cid:0) (cid:1) wherexisthearclengthmeasuredalongtheavalanchetrack,zb denotestheheight of the basal topography relative to the track (usually zb = 0 in one spatial di- mension) and ζ and λκ are the local slope inclination angle and curvature of the track, respectively. The term sgn(u) selects the orientation of the dry Coulomb drag friction, and ε 1 is the aspect ratio of a typical thickness and length of the ≪ avalanche. Note that equations (1) and (2) are written in conservative form [8], while in the original SH theory the smoothness assumption allows the momentum balance equation to transform to an evolution equation for the velocity, viz., du ∂h 1 ∂β x =s β h . (4) x x dt − ∂x − 2 ∂x The factor β is defined as β =εcosζK and the earth pressure coefficient K is x x x x given by the ad hoc assumption K for ∂u/∂x>0, xact K = (5) x  K for ∂u/∂x<0,  xpass with  K =2 1 1 cos2φ/cos2δ sec2φ 1, (6) xact/pass ∓ − − (cid:18) (cid:19) p andφistheinternalfrictionangleofthegranularmaterial. Notethatthevaluesof the earth pressure coefficient K are based on the postulation of a Mohr-Coulomb x SHOCK-CAPTURINGANDFRONT-TRACKINGFORGRANULARAVALANCHES 5 plasticbehaviourforthe cohessionlessyieldonthe basalslidingsurface,seeSavage & Hutter [33, 34] for details. In this theory the earth pressure coefficient K is x assumed to be function of the velocity gradient, i.e. K =K (∂u/∂x). x x Thegoverningequationslookliketheshallow-waterequations,butbecauseofthe jump in the earth pressure coefficients K , the source term s and the free xact/pass x boundary at the front and rearmargins, it becomes much more complicated to de- velopanappropriatenumericalschemetodescribethe flow. TheoriginalLagrange finite-difference scheme [33] is implemented for the equation system (1) and (4) in Lagrangian form, with primitive variables h and u. The shock capturing scheme developed here is applied to the system in conservative form (1) and (2), where the conserved quantities are the avalanche thickness h and the depth integrated momentum m=hu. In vector notation, equations (1) and (2) take the form w +f =s, (7) t x where h m 0 w = , f = and s= . (8) m m2/h+β h2/2 hs (cid:18) (cid:19) (cid:18) x (cid:19) (cid:18) x (cid:19) This form is more convenient for mathematical analysis than (1) and (2). 2.1. Jump Condition and Travelling wave The Savage-Hutter theory can be used to model the upslope propagating travelling shock wave observed in experiments [5, 7] by introducing the jump conditions (see Fig.1) of the balance equations (1) and (2) for mass and momentum [[h(u V )]]=0, (9) n − [[hu(u V )+ 1β h2]]=0, (10) − n 2 x where V is the normalspeedof the singularsurface. Letus suppose that [[β ]]=0 n x (for example, this is always satisfied if φ = δ, i.e., K = K ). Substituting xact xpass (9) into (10) (i.e. eliminating V ) yields the following relation between the depth n ratio, H :=h−/h+, and the velocity difference H +1 H 1 2 (u+ u−)2 =β h− − . (11) x − 2 H (cid:18) (cid:19) For an upslope travelling shock wave with travelling wave speed V and corre- n sponding depth ratio H, the factor β is a function of material and topographic x parameters, φ, δ and ζ, which are given by the selected material and topography. Provided that the depths before and after the shock, h+ and h−, are known (they can be determined by experiment) and the downslope velocity is also given (it is normallyequalto zero),then the upslope velocitycanbe determined byusing (11) 1/2 H 1 H +1 u+ =u− − β h− . (12) x ± H 2 (cid:20) (cid:21) 6 Y.C.TAI,S.NOELLE,J.M.N.T.GRAY&K.HUTTER σ nσ −Vn u− h− h+ u+ ζ FIG. 1. The plane travelling shock wave can be interpreted as a jump in thickness and velocityseparatingthebodyoftheavalancheintotwopartsonaplanewithinclinedangleζ. h+ andh− arethethicknesses ofbothsidesandu+ andu− arethevelocities,respectively, whereas thisjumptravelswithvelocityVn up-slopes. Note that the term under the square root is positive for all positive H. If H = 1 then u+ = u−, which indicates that no shock wave (discontinuity) takes place. Thus, velocity jumps and depth jumps occur together. Byinspectionofthe massbalanceequation(9),thevelocityofthe shockis given by Hu− u+ H +1 1/2 − − V = − = u β h . (13) n H 1 ∓ x 2H2 − (cid:20) (cid:21) Note that as h+ tends to h=h−, then u+ tends to u=u− and V u [β h]1/2, (14) n x → ∓ sowehaverecoveredthe characteristicspeeds ofthe shallowwaterequations. Now we apply Lax’ shock inequalities [24] to single out the physically relevant branches of the shock curves: for the first family, with characteristic speed u √β h, we x − require that 1/2 u+ β h+ 1/2 > V = u− β h−H +1 > u− β h− 1/2, − x n − x 2H2 − x (cid:20) (cid:21) (cid:2) (cid:3) (cid:2) (cid:3) whichimpliesH >1(recallthattheupslopestate“+”liestotheleftoftheshock). Analogously, for the second family, with characteristic speed u+√β h, we obtain x H <1. Forexample,anupwardjump(h+ <h−)canonlybe carriedbya shockof the first family, and in this case u+ >u− >V , so particles which cross the shock n are condensed and slow down. 2.2. Similarity Solution Considerthemotionofafinitemassofgranularmaterialalongaflatplane,i.e. ζ isconstantandλκ=0in(3). In[33]oneparticularsimilaritysolutiontoamoving boundary problem of finite mass was derived; this solution is now generalised (see [36]. To this end we introduce a moving coordinate system with velocity t u (t)=u (0)+ (sinζ tanδ cosζ)dt (15) 0 0 − Z0 SHOCK-CAPTURINGANDFRONT-TRACKINGFORGRANULARAVALANCHES 7 on an plane with inclination angle ζ. This velocity is due to the net driving force s in (3), where we assume that the velocity is positive for positive times, i.e. x sgn(u)= 1. The relative velocity u˘ in the moving coordinate system is then given by u˘=u u (t) . (16) 0 − Asymmetricbulkisconsideredandtheoriginofthemovingcoordinatesystemis selectedtolie atthe centre wherethe surfacegradient,∂h/∂x,is zero. Tokeepthe symmetric depth profile during the motion the relative velocity is further assumed to be skew–symmetric, u˘(ξ, t)= u˘( ξ, t), where − − t ′ ′ ξ =x u (t)dt (17) 0 − Z0 indicates the distance from the origin in the moving coordinates. Provided that g(t)is the distance fromthe coordinateoriginto the marginattime t,the physical domain occupied by the granular mass can be mapped from [ g(t), g(t)] to the − fixed domain [ 1, 1] by − 1 t ′ ′ η = x u (t)dt , where η [ 1,1]. (18) 0 g(t) − ∈ − (cid:26) Z0 (cid:27) With this coordinate mapping, (x, t) (η, τ), the model equations (1) and (2) → reduce to ′ ∂h g ∂h 1 ∂ η + (hu˘)=0, (19) ∂t − g ∂η g∂η ′ ∂u˘ g ∂u˘ 1 ∂u˘ ∂h η + u˘ +β =0, (20) x ∂t − g ∂η g ∂η ∂η (cid:18) (cid:19) ′ where the τ is again replaced by t and we have used g =dg/dt= u /η. 0 − Now we assume that u˘(η,t) varies linearly in η. Since the margins move with ′ ′ relativespeeds g (t),thisyieldsu˘(η,t)=ηg (t). Nowtheevolutionequations(19) ± and (20) reduce to ′ ∂h g + h=0, (21) ∂t g β ∂h ′′ x ηg + =0, (22) g ∂η where g′′ = d2g/dt2. Integrating (22) subject to the boundary condition either h(η =1)=0 or h(η = 1)=0, it follows that the thickness is described by − ′′ g(t)g (t) h(η, t)= (1 η2). (23) 2β − x Thisimpliesthattheavalanchebodykeepsaparabolicthicknessdistributionduring the motion. With the thickness distribution (23) one can easily obtain the total mass M to be ξFt 1 2g′′g2 M = h(ξ, t)dξ = h(η, t)g(t)dη = . (24) Zξf Z−1 3 βx 8 Y.C.TAI,S.NOELLE,J.M.N.T.GRAY&K.HUTTER Since mass is conserved, d 2g ′ ′′ ′′′ 0= M = (2g g +gg ). (25) dt 3β x This relation can also be derived directly from the mass balance equation (21). ′ Changing the independent variable t to g(t) and letting p(t) = g (t), equation (24) can be written as dp K p = , (26) dg g2 where 3β M = 2K. The similarity solution is then obtained by solving (26) with x initial condition, g(0)=g and p(0)=p 0 0 1 1 p2(t)=2K +p2. (27) g − g(t) 0 (cid:18) 0 (cid:19) With the definition α = 2K, β =p2 and G=(α +β )g it follows that g g0 g 0 g g √GG′ =(α +β )3/2. (28) g g √G 2K − We now use the relation d √G √G√G 2K+2Kln(√G+√G 2K) = dG − − √G 2K h i − and integrate equation (28) to yield √G√G 2K+2Kln(√G+√G 2K) − − (29) √G√G 2K+2Kln(√G+√G 2K) =(α +β )3/2t. g g − − − t=0 h i With g =1, p =0 we obtain the Savage-Hutter solution [33] 0 0 √g g 1+ln √g+ g 1 =√2Kt, (30) − − p (cid:16) p (cid:17) for which g(t) > 1. Both (29) and (30) are implicit evolution equations for g(t). ′ Once g(t) is deduced, with the presumption u˘(η,t)=ηg (t), the complete solution is then given by (23) and (27), 1/2 1 1 3M u˘(η, t)=η 2K +p2 , h(η, t)= (1 η2), (31) g − g(t) 0 4g(t) − (cid:26) (cid:18) 0 (cid:19) (cid:27) where η is defined in (18). In the present similarity solution it is presumed that u/u = 1, which means that u > 0 for all t 0. From (16) and the presumption | | ′ ≥ u˘(η,t)=ηg (t) it follows that ′ u(t)=u (t)+u˘(t)>0 g (t)<u (t), for all t 0. (32) 0 0 ⇒ ≥ It is very important to verify that the velocity is consistent with condition (32) to keep the parabolic similarity solution valid. The generalisation (29) of (30) SHOCK-CAPTURINGANDFRONT-TRACKINGFORGRANULARAVALANCHES 9 wasneededto haveexactsolutionswithnon-vanishinginitialvelocities(forfurther details see [36]). 3. NUMERICAL SCHEME The numerical schemes employed in this paper are designed to explicitly solve the system of equations in 1D and we here introduce a Lagrangian algorithm and an Eulerian shock-capturing NOC (Non-Oscillatory Central) scheme. In the Lagrangian technique [33, 34] the avalanche body is divided into several cells. The purpose is to find the velocity of the cell boundaries in order to deter- mine the cell boundary locations for eachtime step; so it is a moving-gridmethod, whereas, the NOC scheme is built on a stationary uniform grid and gives a high resolution of the shock solutions without any spurious oscillations near a discontinuity. In the Lagrangian method the value of the depth hn is defined as the volume j average within the jth cell for time tn, which is bounded by bj−1(t) and bj(t), and the boundary b (t) moves with the velocity u . Whilst, in the NOC scheme the j j valueofthediscretisedvariableUn, U =h, misdefinedonthemeshasthevolume j average within the jth mesh cell centred at position x for time tn, where the jth j cell is bounded by xj+1/2 and xj−1/2. 3.1. Lagrangian method Inthe Lagrangianmethod[33,34]the avalanchebody isdividedintoN material cells, wherex=bj−1(t) andx=bj(t) denote the boundariesofthe cell j attime t, see Fig.2. These boundaries move with the avalanche velocity, i.e. d b (t)=u (t)=u(b (t),t). j j j dt Integrating the mass balance equation (1) over the cell yields bj bj ∂h ∂ d d + (hu) dx= hdx=0 V =0, (33) ∂t ∂x dt ⇒ dt cellj Z (cid:26) (cid:27) Z bj−1 bj−1 and implies that the volume (mass) of the cell is conserved during the motion. Because of this, the mean height of the jth cell can be determined by V hn = cellj . (34) j bn bn j − j−1 The computations proceed as follows. It is assumed that bn, hn and un+1/2 j j j are given as initial values and the new location of the cell boundary bn+1 after an j elapsed time t is given by △ bn+1 =bn+ tun+1/2. (35) j j △ j 10 Y.C.TAI,S.NOELLE,J.M.N.T.GRAY&K.HUTTER h h 3 N-2 h h 2 N-1 h h 1 N b0 c1 b1 c2 b2 c3 b3 bN-3 cN-2 bN-2 cN-1 bN-1 cN bN FIG. 2. The avalanche body isdivided into N elements with average depth hj, where cj isthecentreofthejth element. Notethatherethevelocityu indicatestheboundaryvelocityofb . Themomentum j j balance(4)allowsthevelocityofthecellboundaryattimetn+1/2 tobedetermined, un+1/2 =un−1/2+ t sn εcosζ (K )n ∂h n εhnj+1/2 ∂(cosζKx) n . j j △ ( j − j x j(cid:18)∂x(cid:19)j − 2 (cid:18) ∂x (cid:19)j) (36) The net driving accelerationsn as given by (3) is j sn =sinζ sgn un−1/2 tanδ cosζ +λκ un−1/2 2 εcosζ ∂zb , (37) j j − j j j j − j ∂x (cid:16) (cid:17) (cid:26) (cid:16) (cid:17) (cid:27) (cid:18) (cid:19)j where ζ represents the local inclination angle, κ is the local curvature, and zb j j denotes the local basal topography. Note that the last term at the right-hand side of (36) contains the gradientof the earth pressurecoefficient, which is neglected in the numerical scheme of Savage and Hutter [33, 34]. The earth pressure coefficient K is determined by the ad hoc definition x K , for un−1/2 un−1/2, (K )n = xact j+1 ≥ j (38) x j  K , for un−1/2 < un−1/2  xpass j+1 j in [33, 34]. The surface (depth) gradients in (36) are determined by the depths of the adjacent elements ∂h n (hn hn) 2(hn hn) = j+1− j = j+1− j , (39) ∂x cn cn bn bn (cid:18) (cid:19)j j+1− j j+1− j−1 where cn represents the centre of the jth cell, cn = (bn+bn )/2, at time t = tn, j j j j−1 see Fig.2. The height at the cell boundary, h , is given by their mean values j+1/2 in adjacent cells, h = 1(h +h ), and the gradient of the earth pressure j+1/2 2 j j+1 coefficient is ∂(cosζKx) n = cosζj+1(Kx)nj+1−cosζj(Kx)nj . (40) ∂x cn cn (cid:18) (cid:19)j j+1− j However, while this method is excellent for classical smooth solutions, it looses numerical stability if shocks develop. Shocks are initiated when the avalanche velocity is faster than its characteristic speed and the avalanche front reaches the